xsec_tot_uncert_optimization.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 
37 
39 
42 //#include "NDAna/my_nuecc_inc/NueCCIncBins.h"
43 //#include "CAFAna/XSec/FluxMultiverseSyst.h"
45 
46 #include "TH1.h"
47 #include "TH2.h"
48 #include "THStack.h"
49 #include "TCanvas.h"
50 #include "TLegend.h"
51 #include "TROOT.h"
52 #include <iostream>
53 #include <vector>
54 #include <string>
55 
56 #include "CAFAna/Cuts/SpillCuts.h"
57 #include "CAFAna/Core/Utilities.h"
59 #include "CAFAna/Core/Spectrum.h"
60 #include "CAFAna/Cuts/TimingCuts.h"
61 #include "CAFAna/Analysis/Plots.h"
62 #include "CAFAna/Analysis/Style.h"
64 #include "CAFAna/Core/Ratio.h"
65 
66 #include "TH1D.h"
67 #include "TCanvas.h"
68 #include "TGraph.h"
69 #include "TMultiGraph.h"
70 #include "TLine.h"
71 #include "TLegend.h"
72 
73 using namespace ana;
74 
75 const double pot = 8.09e20; // Data POT for ND after 3A
76 const double mc_pot = 3.54206e+21; //MC POT for ND after 3A
77 //const double mc_syst_pot = mc_pot
78 
79 const bool up_down_combined = true; // If true, do up-down for Genie, Calib and Light syst.
80  // Otherwised, treat all up and down shifts individually
81 
82 //-----------------------------
83 TGraphAsymmErrors* MyPlotWithSystErrorBand(TH1F*& nom,
84  std::map<std::string,TH1F*>& ups,
85  std::map<std::string,TH1F*>& dns,
86  int col, int errCol,
87  float headroom, bool newaxis)
88 {
89  if(col == -1){
90  col = kTotalMCColor;
91  errCol = kTotalMCErrorBandColor;
92  }
93  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
94 
95  TCanvas *c = new TCanvas("c","c");
96 
97  nom->SetLineColor(col);
98  nom->GetXaxis()->CenterTitle();
99  nom->GetYaxis()->CenterTitle();
100  if(newaxis) nom->Draw("hist ]["); // Set the axes up
101 
102  double yMax = nom->GetBinContent(nom->GetMaximumBin());
103 
104  TGraphAsymmErrors* g = new TGraphAsymmErrors;
105  return g; // TODO: handle shape and cherenkov, update maps and vectors
106  /*
107  for(int binIdx = 1; binIdx <= nom->GetNbinsX(); ++binIdx){
108  const double y = nom->GetBinContent(binIdx);
109  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
110 
111  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
112 
113  double errUp = 0;
114  double errDn = 0;
115 
116  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
117  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
118  double lo = y-dns[systIdx]->GetBinContent(binIdx);
119 
120  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
121 
122  errUp += hi*hi;
123  errDn += lo*lo;
124 
125  // TODO: what happens if they're both high or both low?
126  } // end for systIdx
127 
128 
129  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
130  } // end for i
131 
132  g->SetFillColor(errCol);
133  g->Draw("e2 same");
134  g->GetYaxis()->SetRangeUser(0.00001, headroom*yMax);
135  nom->GetYaxis()->SetRangeUser(0.0001, headroom*yMax);
136  //c->SaveAs("OptimizationPlots/Combined_Errors.pdf");
137 
138  TCanvas *c1 = new TCanvas("c1","c1");
139  nom->Draw("hist ][ same");
140  //c1->SaveAs("OptimizationPlots/combined_error_overlay.pdf");
141 
142  return g;
143  */
144 }
145 
146 //--------------------------------------
147 
148 TH1F* GetSelectedStatisticalUncertainty(Spectrum spectrum_nominal, double pot_function)
149 {
150  TH1F* nominal = (TH1F*)spectrum_nominal.ToTH1(pot_function);
151 
152  TH1F* nominal_cut_value = new TH1F("nominal_cut_value",
153  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
154  nominal->GetXaxis()->GetNbins(),
155  nominal->GetXaxis()->GetXmin(),
156  nominal->GetXaxis()->GetXmax());
157 
158 
159  //Integrate To Get Total Number of Events selected for a particular cut value
160  for(int ibin = 1; ibin <= nominal->GetXaxis()->GetNbins(); ibin++){
161  float fNumSelected = nominal->Integral(ibin,nominal->GetXaxis()->GetNbins());
162  nominal_cut_value->SetBinContent(ibin,fNumSelected);
163  }
164 
165  TCanvas *c1 = new TCanvas("c1","c1");
166  nominal_cut_value->Draw();
167  std::string outname = "OptimizationPlots/debug/selected_vs_lower_cut.pdf";
168  c1->SaveAs(outname.c_str());
169  std::string outname_root = "OptimizationPlots/root/selected_vs_lower_cut.root";
170  c1->SaveAs(outname_root.c_str());
171  c1->Close();
172  delete c1;
173 
174  //Now get the statistical error
175  TH1F* stat_err_cut_value = new TH1F("error_cut_value",
176  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
177  nominal->GetXaxis()->GetNbins(),
178  nominal->GetXaxis()->GetXmin(),
179  nominal->GetXaxis()->GetXmax());
180 
181  for(int ibin = 1; ibin <= nominal->GetXaxis()->GetNbins(); ibin++){
182  float fSelectedErr = nominal_cut_value->GetBinContent(ibin);
183  stat_err_cut_value->SetBinContent(ibin,sqrt(fSelectedErr));
184  }
185 
186  TCanvas *c2 = new TCanvas("c2","c2");
187  stat_err_cut_value->Draw();
188  outname = "OptimizationPlots/debug/statistical_error_selected_vs_lower_cut.pdf";
189  c2->SaveAs(outname.c_str());
190  outname_root = "OptimizationPlots/root/statistical_error_selected_vs_lower_cut.root";
191  c2->SaveAs(outname_root.c_str());
192  c2->Close();
193  delete c2;
194 
195  TH1F* fractional_uncertainty = (TH1F*)stat_err_cut_value->Clone("h1");
196  fractional_uncertainty->Divide(nominal_cut_value);
197  TCanvas *c3 = new TCanvas("c3","c3");
198  fractional_uncertainty->Draw();
199  outname = "OptimizationPlots/debug/relative_statistical_error_selected_vs_lower_cut.pdf";
200  c3->SaveAs(outname.c_str());
201  outname_root = "OptimizationPlots/root/relative_statistical_error_selected_vs_lower_cut.root";
202  c3->SaveAs(outname_root.c_str());
203  c3->Close();
204  delete c3;
205 
206  return stat_err_cut_value;
207 
208 }
209 
210 TH1F* GetDenominator( Spectrum spectrum_selected,
211  Spectrum spectrum_background,
212  double pot_function)
213 {
214  TH1F* selected = (TH1F*)spectrum_selected.ToTH1(pot_function);
215  TH1F* background = (TH1F*)spectrum_background.ToTH1(pot_function);
216 
217  TH1F* denominator_cut_value = new TH1F("",
218  TString(";Cut on ") + selected->GetXaxis()->GetTitle(),
219  selected->GetXaxis()->GetNbins(),
220  selected->GetXaxis()->GetXmin(),
221  selected->GetXaxis()->GetXmax());
222 
223 
224  //Integrate To Get Total Number of Events selected for a particular cut value
225  for(int ibin = 1; ibin <= selected->GetXaxis()->GetNbins(); ibin++){
226  float fNumSelected = selected->Integral(ibin,selected->GetXaxis()->GetNbins());
227  float fNumBackground = background->Integral(ibin,background->GetXaxis()->GetNbins());
228  denominator_cut_value->SetBinContent(ibin,fNumSelected - fNumBackground);
229  }
230 
231  TCanvas *c1 = new TCanvas("c1","c1");
232  denominator_cut_value->GetYaxis()->SetTitle("N_{sel}-N_{bkg}");
233  denominator_cut_value->GetXaxis()->SetTitle(TString("Cut on ") + selected->GetXaxis()->GetTitle());
234  denominator_cut_value->Draw();
235  std::string outname = "OptimizationPlots/debug/denominator_vs_lower_cut.pdf";
236  c1->SaveAs(outname.c_str());
237  std::string outname_root = "OptimizationPlots/root/denominator_vs_lower_cut.root";
238  c1->SaveAs(outname_root.c_str());
239  c1->Close();
240  delete c1;
241 
242  return denominator_cut_value;
243 }
244 
245 TH1F* GetBackgroundStatisticalUncertainty(Spectrum spectrum_nominal, double pot_function)
246 {
247  TH1F* nominal = (TH1F*)spectrum_nominal.ToTH1(pot_function);
248 
249  TH1F* nominal_cut_value = new TH1F("nominal_cut_value_bkg",
250  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
251  nominal->GetXaxis()->GetNbins(),
252  nominal->GetXaxis()->GetXmin(),
253  nominal->GetXaxis()->GetXmax());
254 
255 
256  //Integrate To Get Total Number of Events selected for a particular cut value
257  for(int ibin = 1; ibin <= nominal->GetXaxis()->GetNbins(); ibin++){
258  float fNumSelected = nominal->Integral(ibin,nominal->GetXaxis()->GetNbins());
259  nominal_cut_value->SetBinContent(ibin,fNumSelected*pot/mc_pot);
260  }
261 
262  TCanvas *c1 = new TCanvas("c1","c1");
263  nominal_cut_value->GetYaxis()->SetTitle("N_{bkg}");
264  nominal_cut_value->Draw();
265  std::string outname = "OptimizationPlots/debug/background_vs_lower_cut.pdf";
266  c1->SaveAs(outname.c_str());
267  std::string outname_root = "OptimizationPlots/root/background_vs_lower_cut.root";
268  c1->SaveAs(outname_root.c_str());
269  c1->Close();
270  delete c1;
271 
272  //Now get the statistical error
273  TH1F* stat_err_cut_value = new TH1F("error_cut_value_bkg",
274  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
275  nominal->GetXaxis()->GetNbins(),
276  nominal->GetXaxis()->GetXmin(),
277  nominal->GetXaxis()->GetXmax())
278  ;
279  for(int ibin = 1; ibin <= nominal->GetXaxis()->GetNbins(); ibin++){
280  float fSelectedErr = nominal_cut_value->GetBinContent(ibin);
281  stat_err_cut_value->SetBinContent(ibin,sqrt(fSelectedErr*pot/mc_pot));
282  }
283 
284  TCanvas *c2 = new TCanvas("c2","c2");
285  stat_err_cut_value->GetYaxis()->SetTitle("#delta N_{bkg}^{stat}");
286  stat_err_cut_value->Draw();
287  outname = "OptimizationPlots/debug/statistical_error_background_vs_lower_cut.pdf";
288  c2->SaveAs(outname.c_str());
289  outname_root = "OptimizationPlots/root/statistical_error_background_vs_lower_cut.root";
290  c2->SaveAs(outname_root.c_str());
291  c2->Close();
292  delete c2;
293 
294  TH1F* fractional_uncertainty = (TH1F*)stat_err_cut_value->Clone("h1_bkg");
295  fractional_uncertainty->Divide(nominal_cut_value);
296  TCanvas *c3 = new TCanvas("c3","c3");
297  fractional_uncertainty->GetYaxis()->SetTitle("#frac{#delta N_{bkg}^{stat}}{N_{bkg}}");
298  fractional_uncertainty->Draw();
299  outname = "OptimizationPlots/relative_statistical_error_background_vs_lower_cut.pdf";
300  c3->SaveAs(outname.c_str());
301  outname_root = "OptimizationPlots/root/relative_statistical_error_background_vs_lower_cut.root";
302  c3->SaveAs(outname_root.c_str());
303  c3->Close();
304  delete c3;
305 
306  return stat_err_cut_value;
307 }
308 
310  TH1F* denominator,
311  std::map<std::string,Spectrum> ShiftUp,
312  std::map<std::string,Spectrum> ShiftDown,
313  double pot_function)
314 {
315  TH1F* nominal = (TH1F*)spectrum_nominal.ToTH1(pot_function);
316  std::map<std::string,TH1F*> hist_shift_up;
317  std::map<std::string,TH1F*> hist_shift_down;
318 
319  for(auto it = ShiftUp.begin(); it != ShiftUp.end(); ++it)
320  hist_shift_up.insert(std::make_pair(it->first,(TH1F*)ShiftUp.at(it->first).ToTH1(pot_function)));
321 
322  for(auto it = ShiftDown.begin(); it != ShiftDown.end(); ++it)
323  hist_shift_down.insert(std::make_pair(it->first,(TH1F*)ShiftDown.at(it->first).ToTH1(pot_function)));
324 
325  //Need histograms to hold all intermediate values
326  //currently there are 5 systematics shifts in both up and down vectors
327  TH1F* nominal_cut_value = new TH1F("nominal_cut_value_syst",
328  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
329  nominal->GetXaxis()->GetNbins(),
330  nominal->GetXaxis()->GetXmin(),
331  nominal->GetXaxis()->GetXmax());
332  TH1F* genie_up_cut_value = new TH1F("genie_up_cut_value",
333  TString(";Cut on ") + nominal->GetXaxis()->GetTitle() + ";N_{bkg} above cut",
334  nominal->GetXaxis()->GetNbins(),
335  nominal->GetXaxis()->GetXmin(),
336  nominal->GetXaxis()->GetXmax());
337  TH1F* genie_down_cut_value = new TH1F("genie_down_cut_value",
338  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
339  nominal->GetXaxis()->GetNbins(),
340  nominal->GetXaxis()->GetXmin(),
341  nominal->GetXaxis()->GetXmax());
342  TH1F* cal_up_cut_value = new TH1F("cal_up_cut_value",
343  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
344  nominal->GetXaxis()->GetNbins(),
345  nominal->GetXaxis()->GetXmin(),
346  nominal->GetXaxis()->GetXmax());
347  TH1F* cal_down_cut_value = new TH1F("cal_down_cut_value",
348  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
349  nominal->GetXaxis()->GetNbins(),
350  nominal->GetXaxis()->GetXmin(),
351  nominal->GetXaxis()->GetXmax());
352  TH1F* light_up_cut_value = new TH1F("light_up_cut_value",
353  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
354  nominal->GetXaxis()->GetNbins(),
355  nominal->GetXaxis()->GetXmin(),
356  nominal->GetXaxis()->GetXmax());
357  TH1F* light_down_cut_value = new TH1F("light_down_cut_value",
358  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
359  nominal->GetXaxis()->GetNbins(),
360  nominal->GetXaxis()->GetXmin(),
361  nominal->GetXaxis()->GetXmax());
362  TH1F* shape_cut_value = new TH1F("shape_cut_value",
363  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
364  nominal->GetXaxis()->GetNbins(),
365  nominal->GetXaxis()->GetXmin(),
366  nominal->GetXaxis()->GetXmax());
367  TH1F* cheren_cut_value = new TH1F("cheren_cut_value",
368  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
369  nominal->GetXaxis()->GetNbins(),
370  nominal->GetXaxis()->GetXmin(),
371  nominal->GetXaxis()->GetXmax());
372 
373  //Integrate To Get Total Number of Events selected for a particular cut value
374  for(int ibin = 1; ibin <= nominal->GetXaxis()->GetNbins(); ibin++){
375  float fNumSelected = nominal->Integral(ibin,nominal->GetXaxis()->GetNbins());
376  nominal_cut_value->SetBinContent(ibin,fNumSelected);
377 
378  float fGenieUp = hist_shift_up.at("genie")->Integral(ibin,nominal->GetXaxis()->GetNbins());
379  float fGenieDown = hist_shift_down.at("genie")->Integral(ibin,nominal->GetXaxis()->GetNbins());
380  genie_up_cut_value->SetBinContent(ibin, fGenieUp);
381  genie_down_cut_value->SetBinContent(ibin,fGenieDown);
382 
383  float fCalUp = hist_shift_up.at("calib")->Integral(ibin,nominal->GetXaxis()->GetNbins());
384  float fCalDown = hist_shift_down.at("calib")->Integral(ibin,nominal->GetXaxis()->GetNbins());
385  cal_up_cut_value->SetBinContent(ibin, fCalUp);
386  cal_down_cut_value->SetBinContent(ibin,fCalDown);
387 
388  float fLightUp = hist_shift_up.at("light")->Integral(ibin,nominal->GetXaxis()->GetNbins());
389  float fLightDown = hist_shift_down.at("light")->Integral(ibin,nominal->GetXaxis()->GetNbins());
390  light_up_cut_value->SetBinContent(ibin, fLightUp);
391  light_down_cut_value->SetBinContent(ibin,fLightDown);
392 
393  float fShape = hist_shift_up.at("shape")->Integral(ibin,nominal->GetXaxis()->GetNbins());
394  shape_cut_value->SetBinContent(ibin, fShape);
395 
396  float fCheren = hist_shift_up.at("cherenkov")->Integral(ibin,nominal->GetXaxis()->GetNbins());
397  cheren_cut_value->SetBinContent(ibin, fCheren);
398  }
399 
400  // TCanvas *c_syst = new TCanvas("c_syst","c_syst");
401  // nominal->SetLineColor(kBlack);
402 
403  // hist_shift_up.at("genie")->SetLineColor(kRed);
404  // hist_shift_up.at("calib")->SetLineColor(kGreen);
405  // hist_shift_up.at("light")->SetLineColor(kOrange);
406  // hist_shift_up.at("shape")->SetLineColor(kMagenta);
407  // hist_shift_up.at("cherenkov")->SetLineColor(kBlue);
408 
409  // hist_shift_down.at("genie")->SetLineColor(kRed+1);
410  // hist_shift_down.at("calib")->SetLineColor(kGreen+1);
411  // hist_shift_down.at("light")->SetLineColor(kOrange+1);
412 
413  // TLegend *leg = new TLegend(0.65,0.7,0.9,0.9);
414  // leg->AddEntry(nominal,"Nominal");
415  // leg->AddEntry(hist_shift_up.at("genie"),"Genie up");
416  // leg->AddEntry(hist_shift_down.at("genie"),"Genie down");
417  // leg->AddEntry(hist_shift_up.at("calib"),"Cal up");
418  // leg->AddEntry(hist_shift_down.at("calib"),"Cal down");
419  // leg->AddEntry(hist_shift_up.at("light"),"Light up");
420  // leg->AddEntry(hist_shift_down.at("light"),"Light down");
421  // leg->AddEntry(hist_shift_up.at("shape"),"Shape");
422  // leg->AddEntry(hist_shift_up.at("cherenkov"),"Cherenkov");
423 
424  // hist_shift_up.at("genie")->Draw("hist");
425  // hist_shift_down.at("genie")->Draw("histsame");
426  // hist_shift_up.at("calib")->Draw("histsame");
427  // hist_shift_down.at("calib")->Draw("histsame");
428  // hist_shift_up.at("light")->Draw("histsame");
429  // hist_shift_down.at("light")->Draw("histsame");
430  // hist_shift_up.at("shape")->Draw("histsame");
431  // hist_shift_up.at("cherenkov")->Draw("histsame");
432  // nominal->Draw("histsame");
433 
434  // leg->Draw("same");
435  // c_syst->SaveAs("OptimizationPlots/bkg_syst_comparison.pdf");
436  // c_syst->SaveAs("OptimizationPlots/root/bkg_syst_comparison.root");
437  // c_syst->Close();
438  // delete c_syst;
439 
440  TCanvas *c_syst_cut = new TCanvas("c_syst_cut","c_syst_cut");
441  nominal_cut_value->SetLineColor(kBlack);
442 
443  genie_up_cut_value->SetLineColor(kRed);
444  cal_up_cut_value->SetLineColor(kGreen);
445  light_up_cut_value->SetLineColor(kOrange);
446  shape_cut_value->SetLineColor(kMagenta);
447  cheren_cut_value->SetLineColor(kBlue);
448 
449  genie_down_cut_value->SetLineColor(kRed+1);
450  cal_down_cut_value->SetLineColor(kGreen+1);
451  light_down_cut_value->SetLineColor(kOrange+1);
452 
453  TLegend *leg_cut = new TLegend(0.65,0.7,0.9,0.9);
454  leg_cut->AddEntry(nominal_cut_value,"Nominal");
455  leg_cut->AddEntry(genie_up_cut_value,"Genie up");
456  leg_cut->AddEntry(genie_down_cut_value,"Genie down");
457  leg_cut->AddEntry(cal_up_cut_value,"Cal up");
458  leg_cut->AddEntry(cal_down_cut_value,"Cal down");
459  leg_cut->AddEntry(light_up_cut_value,"Light up");
460  leg_cut->AddEntry(light_down_cut_value,"Light down");
461  leg_cut->AddEntry(shape_cut_value,"Shape");
462  leg_cut->AddEntry(cheren_cut_value,"Cherenkov");
463 
464  genie_up_cut_value->Draw("");
465  genie_down_cut_value->Draw("same");
466  cal_up_cut_value->Draw("same");
467  cal_down_cut_value->Draw("same");
468  light_up_cut_value->Draw("same");
469  light_down_cut_value->Draw("same");
470  shape_cut_value->Draw("same");
471  cheren_cut_value->Draw("same");
472  nominal_cut_value->Draw("same");
473  leg_cut->Draw("same");
474  c_syst_cut->SaveAs("OptimizationPlots/bkg_syst_cut_comparison.pdf");
475  c_syst_cut->SaveAs("OptimizationPlots/root/bkg_syst_cut_comparison.root");
476  c_syst_cut->Close();
477  delete c_syst_cut;
478 
479  std::map<std::string,TH1F*> shifts_up;
480  shifts_up.insert(std::make_pair("genie",genie_up_cut_value));
481  shifts_up.insert(std::make_pair("calib",cal_up_cut_value));
482  shifts_up.insert(std::make_pair("light",light_up_cut_value));
483  shifts_up.insert(std::make_pair("shape",shape_cut_value));
484  shifts_up.insert(std::make_pair("cherenkov",cheren_cut_value));
485 
486  std::map<std::string,TH1F*> shifts_down;
487  shifts_down.insert(std::make_pair("genie",genie_down_cut_value));
488  shifts_down.insert(std::make_pair("calib",cal_down_cut_value));
489  shifts_down.insert(std::make_pair("light",light_down_cut_value));
490 
491  new TCanvas();
492  MyPlotWithSystErrorBand(nominal_cut_value, shifts_up, shifts_down, -1,
493  -1, 1.3, true);
494  std::string outname = "OptimizationPlots/debug/background_with_syst_errors_vs_lower_cut.pdf";
495  gPad->SaveAs(outname.c_str());
496  std::string outname_root = "OptimizationPlots/root/background_with_syst_errors_vs_lower_cut.root";
497  gPad->SaveAs(outname_root.c_str());
498 
499  TH1F* syst_err_cut_value_tot = new TH1F("syst_err_cut_value_tot",
500  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
501  nominal->GetXaxis()->GetNbins(),
502  nominal->GetXaxis()->GetXmin(),
503  nominal->GetXaxis()->GetXmax());
504  TH1F* syst_err_cut_value_genie = new TH1F("syst_err_cut_value_genie",
505  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
506  nominal->GetXaxis()->GetNbins(),
507  nominal->GetXaxis()->GetXmin(),
508  nominal->GetXaxis()->GetXmax());
509  TH1F* syst_err_cut_value_calib = new TH1F("syst_err_cut_value_calib",
510  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
511  nominal->GetXaxis()->GetNbins(),
512  nominal->GetXaxis()->GetXmin(),
513  nominal->GetXaxis()->GetXmax());
514  TH1F* syst_err_cut_value_light = new TH1F("syst_err_cut_value_light",
515  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
516  nominal->GetXaxis()->GetNbins(),
517  nominal->GetXaxis()->GetXmin(),
518  nominal->GetXaxis()->GetXmax());
519  TH1F* syst_err_cut_value_shape = new TH1F("syst_err_cut_value_shape",
520  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
521  nominal->GetXaxis()->GetNbins(),
522  nominal->GetXaxis()->GetXmin(),
523  nominal->GetXaxis()->GetXmax());
524  TH1F* syst_err_cut_value_cheren = new TH1F("syst_err_cut_value_cheren",
525  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
526  nominal->GetXaxis()->GetNbins(),
527  nominal->GetXaxis()->GetXmin(),
528  nominal->GetXaxis()->GetXmax());
529  TH1F* syst_err_cut_value_genie_up = new TH1F("syst_err_cut_value_genie_up",
530  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
531  nominal->GetXaxis()->GetNbins(),
532  nominal->GetXaxis()->GetXmin(),
533  nominal->GetXaxis()->GetXmax());
534  TH1F* syst_err_cut_value_calib_up = new TH1F("syst_err_cut_value_calib_up",
535  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
536  nominal->GetXaxis()->GetNbins(),
537  nominal->GetXaxis()->GetXmin(),
538  nominal->GetXaxis()->GetXmax());
539  TH1F* syst_err_cut_value_light_up = new TH1F("syst_err_cut_value_light_up",
540  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
541  nominal->GetXaxis()->GetNbins(),
542  nominal->GetXaxis()->GetXmin(),
543  nominal->GetXaxis()->GetXmax());
544  TH1F* syst_err_cut_value_genie_down = new TH1F("syst_err_cut_value_genie_down",
545  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
546  nominal->GetXaxis()->GetNbins(),
547  nominal->GetXaxis()->GetXmin(),
548  nominal->GetXaxis()->GetXmax());
549  TH1F* syst_err_cut_value_calib_down = new TH1F("syst_err_cut_value_calib_down",
550  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
551  nominal->GetXaxis()->GetNbins(),
552  nominal->GetXaxis()->GetXmin(),
553  nominal->GetXaxis()->GetXmax());
554  TH1F* syst_err_cut_value_light_down = new TH1F("syst_err_cut_value_light_down",
555  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
556  nominal->GetXaxis()->GetNbins(),
557  nominal->GetXaxis()->GetXmin(),
558  nominal->GetXaxis()->GetXmax());
559 
560  if(up_down_combined) {
561  for(int ibin = 1; ibin <= nominal->GetXaxis()->GetNbins(); ibin++){
562 
563  const double i_nominal = nominal_cut_value->GetBinContent(ibin);
564 
565  double i_bkg_syst_uncert = 0;
566  for(auto it_syst = shifts_up.begin(); it_syst != shifts_up.end(); ++it_syst){
567 
568  double isyst_bkg_syst_uncert=0;
569  if(it_syst->first == "shape" || it_syst->first == "cherenkov") // 3==shape or 4==cherenkov
570  isyst_bkg_syst_uncert = shifts_up.at(it_syst->first)->GetBinContent(ibin) - i_nominal;
571  else
572  isyst_bkg_syst_uncert = (shifts_up.at(it_syst->first)->GetBinContent(ibin) -
573  shifts_down.at(it_syst->first)->GetBinContent(ibin)) / 2;
574 
575  i_bkg_syst_uncert += isyst_bkg_syst_uncert*isyst_bkg_syst_uncert;
576  //Debug syst contributions
577  if(it_syst->first == "genie")
578  syst_err_cut_value_genie->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
579  if(it_syst->first == "calib")
580  syst_err_cut_value_calib->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
581  if(it_syst->first == "light")
582  syst_err_cut_value_light->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
583  if(it_syst->first == "shape")
584  syst_err_cut_value_shape->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
585  if(it_syst->first == "cherenkov")
586  syst_err_cut_value_cheren->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
587  } // end for systIdx
588 
589  syst_err_cut_value_tot->SetBinContent(ibin, sqrt(std::abs(i_bkg_syst_uncert)));
590  }
591 
592  TCanvas *c2 = new TCanvas("c2","c2");
593  syst_err_cut_value_tot->Draw();
594  outname = "OptimizationPlots/systematic_error_background_vs_lower_cut.pdf";
595  c2->SaveAs(outname.c_str());
596  outname_root = "OptimizationPlots/root/systematic_error_background_vs_lower_cut.root";
597  c2->SaveAs(outname_root.c_str());
598  c2->Close();
599  delete c2;
600 
601  TCanvas *c_bkg_syst_decomp = new TCanvas("c_bkg_syst_decomp","c_bkg_syst_decomp");
602  syst_err_cut_value_tot->GetYaxis()->SetTitle("#delta N_{bkg}");
603  syst_err_cut_value_tot->SetLineColor(kBlack);
604  syst_err_cut_value_genie->SetLineColor(kRed);
605  syst_err_cut_value_calib->SetLineColor(kGreen);
606  syst_err_cut_value_light->SetLineColor(kOrange);
607  syst_err_cut_value_shape->SetLineColor(kMagenta);
608  syst_err_cut_value_cheren->SetLineColor(kBlue);
609 
610  TLegend *leg_eff_decomp = new TLegend(0.65,0.7,0.9,0.9);
611  leg_eff_decomp->AddEntry(syst_err_cut_value_tot,"Total");
612  leg_eff_decomp->AddEntry(syst_err_cut_value_genie,"Genie");
613  leg_eff_decomp->AddEntry(syst_err_cut_value_calib,"Calib");
614  leg_eff_decomp->AddEntry(syst_err_cut_value_light,"Light");
615  leg_eff_decomp->AddEntry(syst_err_cut_value_shape,"Shape");
616  leg_eff_decomp->AddEntry(syst_err_cut_value_cheren,"Cherenkov");
617 
618  syst_err_cut_value_tot->Draw("hist");
619  syst_err_cut_value_genie->Draw("histsame");
620  syst_err_cut_value_calib->Draw("histsame");
621  syst_err_cut_value_light->Draw("histsame");
622  syst_err_cut_value_shape->Draw("histsame");
623  syst_err_cut_value_cheren->Draw("histsame");
624  leg_eff_decomp->Draw("same");
625  c_bkg_syst_decomp->SaveAs("OptimizationPlots/delta_bkg_syst_decomp.pdf");
626  c_bkg_syst_decomp->SaveAs("OptimizationPlots/root/delta_bkg_syst_decomp.root");
627 
628  TH1* copy_syst_err_cut_value_tot = (TH1F*)syst_err_cut_value_tot->Clone("copy_syst_err_cut_value_tot");
629  c_bkg_syst_decomp->SetLeftMargin(0.15);
630  copy_syst_err_cut_value_tot->GetYaxis()->SetTitle("#frac{#delta N_{bkg}}{N_{sel}-N_{bkg}}");
631  copy_syst_err_cut_value_tot->GetYaxis()->SetRangeUser(0,1);
632  copy_syst_err_cut_value_tot->Divide(denominator);
633  syst_err_cut_value_genie->Divide(denominator);
634  syst_err_cut_value_calib->Divide(denominator);
635  syst_err_cut_value_light->Divide(denominator);
636  syst_err_cut_value_shape->Divide(denominator);
637  syst_err_cut_value_cheren->Divide(denominator);
638 
639  copy_syst_err_cut_value_tot->Draw("hist");
640  syst_err_cut_value_genie->Draw("histsame");
641  syst_err_cut_value_calib->Draw("histsame");
642  syst_err_cut_value_light->Draw("histsame");
643  syst_err_cut_value_shape->Draw("histsame");
644  syst_err_cut_value_cheren->Draw("histsame");
645  leg_eff_decomp->Draw("same");
646  c_bkg_syst_decomp->SaveAs("OptimizationPlots/bkg_syst_decomp.pdf");
647  c_bkg_syst_decomp->SaveAs("OptimizationPlots/root/bkg_syst_decomp.root");
648 
649  c_bkg_syst_decomp->Close();
650  delete c_bkg_syst_decomp;
651  }
652  else {
653  for(int ibin = 1; ibin <= nominal->GetXaxis()->GetNbins(); ibin++){
654 
655  const double i_nominal = nominal_cut_value->GetBinContent(ibin);
656 
657  double i_bkg_syst_uncert = 0;
658 
659  for(auto it_syst = shifts_up.begin(); it_syst != shifts_up.end(); ++it_syst){
660 
661  double isyst_bkg_syst_uncert=0;
662 
663  isyst_bkg_syst_uncert = shifts_up.at(it_syst->first)->GetBinContent(ibin) - i_nominal;
664 
665  i_bkg_syst_uncert += isyst_bkg_syst_uncert*isyst_bkg_syst_uncert;
666 
667  //Debug syst contributions
668  if(it_syst->first == "genie")
669  syst_err_cut_value_genie_up->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
670  if(it_syst->first == "calib")
671  syst_err_cut_value_calib_up->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
672  if(it_syst->first == "light")
673  syst_err_cut_value_light_up->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
674  if(it_syst->first == "shape")
675  syst_err_cut_value_shape->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
676  if(it_syst->first == "cherenkov")
677  syst_err_cut_value_cheren->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
678  } // end for systIdx
679 
680  for(auto it_syst = shifts_down.begin(); it_syst != shifts_down.end(); ++it_syst){
681 
682  double isyst_bkg_syst_uncert=0;
683 
684  isyst_bkg_syst_uncert = shifts_down.at(it_syst->first)->GetBinContent(ibin) - i_nominal;
685 
686  i_bkg_syst_uncert += isyst_bkg_syst_uncert*isyst_bkg_syst_uncert;
687 
688  //Debug syst contributions
689  if(it_syst->first == "genie")
690  syst_err_cut_value_genie_down->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
691  if(it_syst->first == "calib")
692  syst_err_cut_value_calib_down->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
693  if(it_syst->first == "light")
694  syst_err_cut_value_light_down->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
695  } // end for systIdx
696 
697  syst_err_cut_value_tot->SetBinContent(ibin, sqrt(std::abs(i_bkg_syst_uncert)));
698  }
699 
700  TCanvas *c2 = new TCanvas("c2","c2");
701  syst_err_cut_value_tot->Draw();
702  outname = "OptimizationPlots/systematic_error_background_vs_lower_cut.pdf";
703  c2->SaveAs(outname.c_str());
704  outname_root = "OptimizationPlots/root/systematic_error_background_vs_lower_cut.root";
705  c2->SaveAs(outname_root.c_str());
706  c2->Close();
707  delete c2;
708 
709  TCanvas *c_bkg_syst_decomp = new TCanvas("c_bkg_syst_decomp","c_bkg_syst_decomp");
710  syst_err_cut_value_tot->GetYaxis()->SetTitle("#delta N_{bkg}");
711  syst_err_cut_value_tot->SetLineColor(kBlack);
712  syst_err_cut_value_genie_up->SetLineColor(kRed);
713  syst_err_cut_value_calib_up->SetLineColor(kGreen);
714  syst_err_cut_value_light_up->SetLineColor(kOrange);
715  syst_err_cut_value_genie_down->SetLineColor(kRed+1);
716  syst_err_cut_value_calib_down->SetLineColor(kGreen+1);
717  syst_err_cut_value_light_down->SetLineColor(kOrange+1);
718  syst_err_cut_value_shape->SetLineColor(kMagenta);
719  syst_err_cut_value_cheren->SetLineColor(kBlue);
720 
721  TLegend *leg_eff_decomp = new TLegend(0.65,0.7,0.9,0.9);
722  leg_eff_decomp->AddEntry(syst_err_cut_value_tot,"Total");
723  leg_eff_decomp->AddEntry(syst_err_cut_value_genie_up,"Genie up");
724  leg_eff_decomp->AddEntry(syst_err_cut_value_calib_up,"Cal up");
725  leg_eff_decomp->AddEntry(syst_err_cut_value_light_up,"Light up");
726  leg_eff_decomp->AddEntry(syst_err_cut_value_genie_down,"Genie down");
727  leg_eff_decomp->AddEntry(syst_err_cut_value_calib_down,"Cal down");
728  leg_eff_decomp->AddEntry(syst_err_cut_value_light_down,"Light down");
729  leg_eff_decomp->AddEntry(syst_err_cut_value_shape,"Shape");
730  leg_eff_decomp->AddEntry(syst_err_cut_value_cheren,"Cherenkov");
731 
732  syst_err_cut_value_tot->Draw("hist");
733  syst_err_cut_value_genie_up->Draw("histsame");
734  syst_err_cut_value_calib_up->Draw("histsame");
735  syst_err_cut_value_light_up->Draw("histsame");
736  syst_err_cut_value_genie_down->Draw("histsame");
737  syst_err_cut_value_calib_down->Draw("histsame");
738  syst_err_cut_value_light_down->Draw("histsame");
739  syst_err_cut_value_shape->Draw("histsame");
740  syst_err_cut_value_cheren->Draw("histsame");
741  leg_eff_decomp->Draw("same");
742  c_bkg_syst_decomp->SaveAs("OptimizationPlots/delta_bkg_syst_decomp.pdf");
743  c_bkg_syst_decomp->SaveAs("OptimizationPlots/root/delta_bkg_syst_decomp.root");
744 
745  TH1* copy_syst_err_cut_value_tot = (TH1F*)syst_err_cut_value_tot->Clone("copy_syst_err_cut_value_tot");
746  c_bkg_syst_decomp->SetLeftMargin(0.15);
747  copy_syst_err_cut_value_tot->GetYaxis()->SetTitle("#frac{#delta N_{bkg}}{N_{sel}-N_{bkg}}");
748  copy_syst_err_cut_value_tot->GetYaxis()->SetRangeUser(0,1);
749  copy_syst_err_cut_value_tot->Divide(denominator);
750  syst_err_cut_value_genie_up->Divide(denominator);
751  syst_err_cut_value_calib_up->Divide(denominator);
752  syst_err_cut_value_light_up->Divide(denominator);
753  syst_err_cut_value_genie_down->Divide(denominator);
754  syst_err_cut_value_calib_down->Divide(denominator);
755  syst_err_cut_value_light_down->Divide(denominator);
756  syst_err_cut_value_shape->Divide(denominator);
757  syst_err_cut_value_cheren->Divide(denominator);
758 
759  copy_syst_err_cut_value_tot->Draw("hist");
760  syst_err_cut_value_genie_up->Draw("histsame");
761  syst_err_cut_value_calib_up->Draw("histsame");
762  syst_err_cut_value_light_up->Draw("histsame");
763  syst_err_cut_value_genie_down->Draw("histsame");
764  syst_err_cut_value_calib_down->Draw("histsame");
765  syst_err_cut_value_light_down->Draw("histsame");
766  syst_err_cut_value_shape->Draw("histsame");
767  syst_err_cut_value_cheren->Draw("histsame");
768  leg_eff_decomp->Draw("same");
769  c_bkg_syst_decomp->SaveAs("OptimizationPlots/bkg_syst_decomp.pdf");
770  c_bkg_syst_decomp->SaveAs("OptimizationPlots/root/bkg_syst_decomp.root");
771 
772  c_bkg_syst_decomp->Close();
773  delete c_bkg_syst_decomp;
774 
775  }
776  return syst_err_cut_value_tot;
777 }
778 
780  Spectrum spec_truth_nominal,
781  std::map<std::string,Spectrum> ShiftUp,
782  std::map<std::string,Spectrum> ShiftDown,
783  std::map<std::string,Spectrum> TruthUp,
784  std::map<std::string,Spectrum> TruthDown,
785  double pot_function)
786 {
787  TH1F* signal_nominal = (TH1F*)spec_signal_nominal.ToTH1(pot_function);
788  TH1F* truth_nominal = (TH1F*)spec_truth_nominal.ToTH1(pot_function);
789 
790  std::map<std::string,TH1F*> hist_shift_up;
791  std::map<std::string,TH1F*> hist_shift_down;
792  std::map<std::string,TH1F*> truth_shift_up;
793  std::map<std::string,TH1F*> truth_shift_down;
794 
795  for(auto it = ShiftUp.begin(); it != ShiftUp.end(); ++it){
796  hist_shift_up.insert(std::make_pair(it->first,(TH1F*)ShiftUp.at(it->first).ToTH1(pot_function)));
797  truth_shift_up.insert(std::make_pair(it->first,(TH1F*)TruthUp.at(it->first).ToTH1(pot_function)));
798  }
799  for(auto it = ShiftDown.begin(); it != ShiftDown.end(); ++it) {
800  hist_shift_down.insert(std::make_pair(it->first,(TH1F*)ShiftDown.at(it->first).ToTH1(pot_function)));
801  truth_shift_down.insert(std::make_pair(it->first,(TH1F*)TruthDown.at(it->first).ToTH1(pot_function)));
802  }
803 
804  //Need histograms to hold all intermediate values
805  //currently there are 5 systematics shifts in both up and down vectors
806  TH1F* nominal_cut_value_eff = new TH1F("nominal_cut_value_syst_eff",
807  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
808  signal_nominal->GetXaxis()->GetNbins(),
809  signal_nominal->GetXaxis()->GetXmin(),
810  signal_nominal->GetXaxis()->GetXmax());
811  TH1F* genie_up_cut_value_eff = new TH1F("genie_up_cut_value_eff",
812  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle() + ";#epsilon",
813  signal_nominal->GetXaxis()->GetNbins(),
814  signal_nominal->GetXaxis()->GetXmin(),
815  signal_nominal->GetXaxis()->GetXmax());
816  TH1F* genie_down_cut_value_eff = new TH1F("genie_down_cut_value_eff",
817  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
818  signal_nominal->GetXaxis()->GetNbins(),
819  signal_nominal->GetXaxis()->GetXmin(),
820  signal_nominal->GetXaxis()->GetXmax());
821  TH1F* cal_up_cut_value_eff = new TH1F("cal_up_cut_value_eff",
822  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
823  signal_nominal->GetXaxis()->GetNbins(),
824  signal_nominal->GetXaxis()->GetXmin(),
825  signal_nominal->GetXaxis()->GetXmax());
826  TH1F* cal_down_cut_value_eff = new TH1F("cal_down_cut_value_eff",
827  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
828  signal_nominal->GetXaxis()->GetNbins(),
829  signal_nominal->GetXaxis()->GetXmin(),
830  signal_nominal->GetXaxis()->GetXmax());
831  TH1F* light_up_cut_value_eff = new TH1F("light_up_cut_value_eff",
832  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
833  signal_nominal->GetXaxis()->GetNbins(),
834  signal_nominal->GetXaxis()->GetXmin(),
835  signal_nominal->GetXaxis()->GetXmax());
836  TH1F* light_down_cut_value_eff = new TH1F("light_down_cut_value_eff",
837  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
838  signal_nominal->GetXaxis()->GetNbins(),
839  signal_nominal->GetXaxis()->GetXmin(),
840  signal_nominal->GetXaxis()->GetXmax());
841  TH1F* shape_cut_value_eff = new TH1F("shape_cut_value_eff",
842  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
843  signal_nominal->GetXaxis()->GetNbins(),
844  signal_nominal->GetXaxis()->GetXmin(),
845  signal_nominal->GetXaxis()->GetXmax());
846  TH1F* cheren_cut_value_eff = new TH1F("cheren_cut_value_eff",
847  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
848  signal_nominal->GetXaxis()->GetNbins(),
849  signal_nominal->GetXaxis()->GetXmin(),
850  signal_nominal->GetXaxis()->GetXmax());
851 
852  //Integrate To Get Total Number of Events selected for a particular cut value
853  //Divide By Each samples total integral to convert to efficiency
854 
855  double denominator_nominal = truth_nominal->Integral(1,truth_nominal->GetXaxis()->GetNbins());
856  double denominator_genie_up = truth_shift_up.at("genie")->Integral(1,truth_shift_up.at("genie")->GetXaxis()->GetNbins());
857  double denominator_genie_down = truth_shift_down.at("genie")->Integral(1,truth_shift_down.at("genie")->GetXaxis()->GetNbins());
858  double denominator_calib_up = truth_shift_up.at("calib")->Integral(1,truth_shift_up.at("calib")->GetXaxis()->GetNbins());
859  double denominator_calib_down = truth_shift_down.at("calib")->Integral(1,truth_shift_down.at("calib")->GetXaxis()->GetNbins());
860  double denominator_light_up = truth_shift_up.at("light")->Integral(1,truth_shift_up.at("light")->GetXaxis()->GetNbins());
861  double denominator_light_down = truth_shift_down.at("light")->Integral(1,truth_shift_down.at("light")->GetXaxis()->GetNbins());
862  double denominator_shape = truth_shift_up.at("shape")->Integral(1,truth_shift_up.at("shape")->GetXaxis()->GetNbins());
863  double denominator_cherenkov = truth_shift_up.at("cherenkov")->Integral(1,truth_shift_up.at("cherenkov")->GetXaxis()->GetNbins());
864 
865  for(int ibin = 1; ibin <= signal_nominal->GetXaxis()->GetNbins(); ibin++){
866  nominal_cut_value_eff->SetBinContent(ibin,
867  signal_nominal->Integral(ibin,signal_nominal->GetXaxis()->GetNbins())/denominator_nominal);
868 
869  genie_up_cut_value_eff->SetBinContent(ibin,
870  hist_shift_up.at("genie")->Integral(ibin,hist_shift_up.at("genie")->GetXaxis()->GetNbins())/denominator_genie_up);
871  genie_down_cut_value_eff->SetBinContent(ibin,
872  hist_shift_down.at("genie")->Integral(ibin,hist_shift_down.at("genie")->GetXaxis()->GetNbins())/denominator_genie_down);
873  cal_up_cut_value_eff->SetBinContent(ibin,
874  hist_shift_up.at("calib")->Integral(ibin,hist_shift_up.at("calib")->GetXaxis()->GetNbins())/denominator_calib_up);
875  cal_down_cut_value_eff->SetBinContent(ibin,
876  hist_shift_down.at("calib")->Integral(ibin,hist_shift_down.at("calib")->GetXaxis()->GetNbins())/denominator_calib_down);
877  light_up_cut_value_eff->SetBinContent(ibin,
878  hist_shift_up.at("light")->Integral(ibin,hist_shift_up.at("light")->GetXaxis()->GetNbins())/denominator_light_up);
879  light_down_cut_value_eff->SetBinContent(ibin,
880  hist_shift_down.at("light")->Integral(ibin,hist_shift_down.at("light")->GetXaxis()->GetNbins())/denominator_light_down);
881  shape_cut_value_eff->SetBinContent(ibin,
882  hist_shift_up.at("shape")->Integral(ibin,hist_shift_up.at("shape")->GetXaxis()->GetNbins())/denominator_shape);
883  cheren_cut_value_eff->SetBinContent(ibin,
884  hist_shift_up.at("cherenkov")->Integral(ibin,hist_shift_up.at("cherenkov")->GetXaxis()->GetNbins())/denominator_cherenkov);
885  }
886 
887  TCanvas *c_eff = new TCanvas("c_eff","c_eff");
888 
889  nominal_cut_value_eff->SetLineColor(kBlack);
890 
891  genie_up_cut_value_eff->SetLineColor(kRed);
892  cal_up_cut_value_eff->SetLineColor(kGreen);
893  light_up_cut_value_eff->SetLineColor(kOrange);
894  shape_cut_value_eff->SetLineColor(kMagenta);
895  cheren_cut_value_eff->SetLineColor(kBlue);
896 
897  genie_down_cut_value_eff->SetLineColor(kRed+1);
898  cal_down_cut_value_eff->SetLineColor(kGreen+1);
899  light_down_cut_value_eff->SetLineColor(kOrange+1);
900 
901  TLegend *leg = new TLegend(0.65,0.7,0.9,0.9);
902  leg->AddEntry(nominal_cut_value_eff,"Nominal");
903  leg->AddEntry(genie_up_cut_value_eff,"Genie up");
904  leg->AddEntry(genie_down_cut_value_eff,"Genie down");
905  leg->AddEntry(cal_up_cut_value_eff,"Cal up");
906  leg->AddEntry(cal_down_cut_value_eff,"Cal down");
907  leg->AddEntry(light_up_cut_value_eff,"Light up");
908  leg->AddEntry(light_down_cut_value_eff,"Light down");
909  leg->AddEntry(shape_cut_value_eff,"Shape up");
910  leg->AddEntry(cheren_cut_value_eff,"Cherenkov up");
911 
912  nominal_cut_value_eff->GetYaxis()->SetRangeUser(0,0.08);
913  nominal_cut_value_eff->GetYaxis()->SetTitle("#epsilon");
914  genie_up_cut_value_eff->Draw("");
915  genie_down_cut_value_eff->Draw("same");
916  cal_up_cut_value_eff->Draw("same");
917  cal_down_cut_value_eff->Draw("same");
918  light_up_cut_value_eff->Draw("same");
919  light_down_cut_value_eff->Draw("same");
920  shape_cut_value_eff->Draw("same");
921  cheren_cut_value_eff->Draw("same");
922  nominal_cut_value_eff->Draw("same");
923  leg->Draw("same");
924  c_eff->SaveAs("OptimizationPlots/eff_syst_comparison.pdf");
925  c_eff->SaveAs("OptimizationPlots/root/eff_syst_comparison.root");
926  c_eff->Close();
927  delete c_eff;
928 
929  std::map<std::string,TH1F*> shifts_up;
930  shifts_up.insert(std::make_pair("genie",genie_up_cut_value_eff));
931  shifts_up.insert(std::make_pair("calib",cal_up_cut_value_eff));
932  shifts_up.insert(std::make_pair("light",light_up_cut_value_eff));
933  shifts_up.insert(std::make_pair("shape",shape_cut_value_eff));
934  shifts_up.insert(std::make_pair("cherenkov",cheren_cut_value_eff));
935 
936  std::map<std::string,TH1F*> shifts_down;
937  shifts_down.insert(std::make_pair("genie",genie_down_cut_value_eff));
938  shifts_down.insert(std::make_pair("calib",cal_down_cut_value_eff));
939  shifts_down.insert(std::make_pair("light",light_down_cut_value_eff));
940 
941  new TCanvas();
942  MyPlotWithSystErrorBand(nominal_cut_value_eff, shifts_up,
943  shifts_down, -1,
944  -1, 1.3, true);
945  std::string outname = "OptimizationPlots/debug/efficiency_with_syst_errors_vs_lower_cut.pdf";
946  gPad->SaveAs(outname.c_str());
947  std::string outname_root = "OptimizationPlots/root/efficiency_with_syst_errors_vs_lower_cut.root";
948  gPad->SaveAs(outname_root.c_str());
949 
950  TH1F* syst_err_cut_value_eff_tot = new TH1F("syst_err_cut_value_eff_tot",
951  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
952  signal_nominal->GetXaxis()->GetNbins(),
953  signal_nominal->GetXaxis()->GetXmin(),
954  signal_nominal->GetXaxis()->GetXmax());
955 
956  TH1F* syst_err_cut_value_eff_truth = new TH1F("syst_err_cut_value_eff_truth",
957  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
958  signal_nominal->GetXaxis()->GetNbins(),
959  signal_nominal->GetXaxis()->GetXmin(),
960  signal_nominal->GetXaxis()->GetXmax());
961 
962  TH1F* syst_err_cut_value_eff_genie = new TH1F("syst_err_cut_value_eff_genie",
963  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
964  signal_nominal->GetXaxis()->GetNbins(),
965  signal_nominal->GetXaxis()->GetXmin(),
966  signal_nominal->GetXaxis()->GetXmax());
967  TH1F* syst_err_cut_value_eff_calib = new TH1F("syst_err_cut_value_eff_calib",
968  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
969  signal_nominal->GetXaxis()->GetNbins(),
970  signal_nominal->GetXaxis()->GetXmin(),
971  signal_nominal->GetXaxis()->GetXmax());
972  TH1F* syst_err_cut_value_eff_light = new TH1F("syst_err_cut_value_eff_light",
973  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
974  signal_nominal->GetXaxis()->GetNbins(),
975  signal_nominal->GetXaxis()->GetXmin(),
976  signal_nominal->GetXaxis()->GetXmax());
977  TH1F* syst_err_cut_value_eff_shape = new TH1F("syst_err_cut_value_eff_shape",
978  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
979  signal_nominal->GetXaxis()->GetNbins(),
980  signal_nominal->GetXaxis()->GetXmin(),
981  signal_nominal->GetXaxis()->GetXmax());
982  TH1F* syst_err_cut_value_eff_cheren = new TH1F("syst_err_cut_value_eff_cheren",
983  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
984  signal_nominal->GetXaxis()->GetNbins(),
985  signal_nominal->GetXaxis()->GetXmin(),
986  signal_nominal->GetXaxis()->GetXmax());
987 
988  TH1F* syst_err_cut_value_eff_genie_up = new TH1F("syst_err_cut_value_eff_genie_up",
989  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
990  signal_nominal->GetXaxis()->GetNbins(),
991  signal_nominal->GetXaxis()->GetXmin(),
992  signal_nominal->GetXaxis()->GetXmax());
993  TH1F* syst_err_cut_value_eff_calib_up = new TH1F("syst_err_cut_value_eff_calib_up",
994  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
995  signal_nominal->GetXaxis()->GetNbins(),
996  signal_nominal->GetXaxis()->GetXmin(),
997  signal_nominal->GetXaxis()->GetXmax());
998  TH1F* syst_err_cut_value_eff_light_up = new TH1F("syst_err_cut_value_eff_light_up",
999  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
1000  signal_nominal->GetXaxis()->GetNbins(),
1001  signal_nominal->GetXaxis()->GetXmin(),
1002  signal_nominal->GetXaxis()->GetXmax());
1003  TH1F* syst_err_cut_value_eff_genie_down = new TH1F("syst_err_cut_value_eff_genie_down",
1004  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
1005  signal_nominal->GetXaxis()->GetNbins(),
1006  signal_nominal->GetXaxis()->GetXmin(),
1007  signal_nominal->GetXaxis()->GetXmax());
1008  TH1F* syst_err_cut_value_eff_calib_down = new TH1F("syst_err_cut_value_eff_calib_down",
1009  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
1010  signal_nominal->GetXaxis()->GetNbins(),
1011  signal_nominal->GetXaxis()->GetXmin(),
1012  signal_nominal->GetXaxis()->GetXmax());
1013  TH1F* syst_err_cut_value_eff_light_down = new TH1F("syst_err_cut_value_eff_light_down",
1014  TString(";Cut on ") + signal_nominal->GetXaxis()->GetTitle(),
1015  signal_nominal->GetXaxis()->GetNbins(),
1016  signal_nominal->GetXaxis()->GetXmin(),
1017  signal_nominal->GetXaxis()->GetXmax());
1018 
1019  if(up_down_combined) {
1020  for(int ibin = 1; ibin <= signal_nominal->GetXaxis()->GetNbins(); ibin++){
1021 
1022  const double i_nominal = nominal_cut_value_eff->GetBinContent(ibin);
1023 
1024  double i_bkg_syst_uncert = 0;
1025 
1026  for(auto it_syst = shifts_up.begin(); it_syst != shifts_up.end(); ++it_syst){
1027 
1028  double isyst_bkg_syst_uncert=0;
1029 
1030  if(it_syst->first == "shape" || it_syst->first == "cherenkov") // 3==shape or 4==cherenkov
1031  isyst_bkg_syst_uncert = shifts_up.at(it_syst->first)->GetBinContent(ibin) - i_nominal;
1032  else
1033  isyst_bkg_syst_uncert = (shifts_up.at(it_syst->first)->GetBinContent(ibin) -
1034  shifts_down.at(it_syst->first)->GetBinContent(ibin)) / 2;
1035 
1036  i_bkg_syst_uncert += isyst_bkg_syst_uncert*isyst_bkg_syst_uncert;
1037 
1038  //Debug syst contributions
1039  if(it_syst->first == "genie")
1040  syst_err_cut_value_eff_genie->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
1041  if(it_syst->first == "calib")
1042  syst_err_cut_value_eff_calib->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
1043  if(it_syst->first == "light")
1044  syst_err_cut_value_eff_light->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
1045  if(it_syst->first == "shape")
1046  syst_err_cut_value_eff_shape->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
1047  if(it_syst->first == "cherenkov")
1048  syst_err_cut_value_eff_cheren->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
1049  } // end for syst loop
1050 
1051  syst_err_cut_value_eff_tot->SetBinContent(ibin, sqrt(std::abs(i_bkg_syst_uncert)));
1052  }
1053 
1054  TCanvas *c2 = new TCanvas("c2","c2");
1055  syst_err_cut_value_eff_tot->Draw();
1056  outname = "OptimizationPlots/debug/systematic_error_efficiency_vs_lower_cut.pdf";
1057  c2->SaveAs(outname.c_str());
1058  outname_root = "OptimizationPlots/root/systematic_error_efficiency_vs_lower_cut.root";
1059  c2->SaveAs(outname_root.c_str());
1060  c2->Close();
1061  delete c2;
1062 
1063  TCanvas *c_eff_decomp = new TCanvas("c_eff_decomp","c_eff_decomp");
1064  syst_err_cut_value_eff_tot->GetYaxis()->SetTitle("#delta#epsilon");
1065  syst_err_cut_value_eff_tot->SetLineColor(kBlack);
1066  syst_err_cut_value_eff_genie->SetLineColor(kRed);
1067  syst_err_cut_value_eff_calib->SetLineColor(kGreen);
1068  syst_err_cut_value_eff_light->SetLineColor(kOrange);
1069  syst_err_cut_value_eff_shape->SetLineColor(kMagenta);
1070  syst_err_cut_value_eff_cheren->SetLineColor(kBlue);
1071 
1072  TLegend *leg_eff_decomp = new TLegend(0.65,0.7,0.9,0.9);
1073  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_tot,"Total");
1074  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_genie,"Genie");
1075  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_calib,"Calib");
1076  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_light,"Light");
1077  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_shape,"Shape");
1078  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_cheren,"Cherenkov");
1079 
1080  syst_err_cut_value_eff_tot->Draw("hist");
1081  syst_err_cut_value_eff_genie->Draw("histsame");
1082  syst_err_cut_value_eff_calib->Draw("histsame");
1083  syst_err_cut_value_eff_light->Draw("histsame");
1084  syst_err_cut_value_eff_shape->Draw("histsame");
1085  syst_err_cut_value_eff_cheren->Draw("histsame");
1086  leg_eff_decomp->Draw("same");
1087  c_eff_decomp->SaveAs("OptimizationPlots/efficiency_syst_error_numerator.pdf");
1088  c_eff_decomp->SaveAs("OptimizationPlots/root/efficiency_syst_error_numerator.root");
1089 
1090  TH1* copy_syst_err_cut_value_eff_tot = (TH1F*)syst_err_cut_value_eff_tot->Clone("copy_syst_err_cut_value_eff_tot");
1091  copy_syst_err_cut_value_eff_tot->GetYaxis()->SetTitle("#frac{#delta#epsilon}{#epsilon}");
1092  copy_syst_err_cut_value_eff_tot->GetYaxis()->SetRangeUser(0,1);
1093  copy_syst_err_cut_value_eff_tot->Divide(nominal_cut_value_eff);
1094  syst_err_cut_value_eff_genie->Divide(nominal_cut_value_eff);
1095  syst_err_cut_value_eff_calib->Divide(nominal_cut_value_eff);
1096  syst_err_cut_value_eff_light->Divide(nominal_cut_value_eff);
1097  syst_err_cut_value_eff_shape->Divide(nominal_cut_value_eff);
1098  syst_err_cut_value_eff_cheren->Divide(nominal_cut_value_eff);
1099 
1100  copy_syst_err_cut_value_eff_tot->Draw("hist");
1101  syst_err_cut_value_eff_genie->Draw("histsame");
1102  syst_err_cut_value_eff_calib->Draw("histsame");
1103  syst_err_cut_value_eff_light->Draw("histsame");
1104  syst_err_cut_value_eff_shape->Draw("histsame");
1105  syst_err_cut_value_eff_cheren->Draw("histsame");
1106  leg_eff_decomp->Draw("same");
1107  c_eff_decomp->SaveAs("OptimizationPlots/efficiency_syst_error_decomp.pdf");
1108  c_eff_decomp->SaveAs("OptimizationPlots/root/efficiency_syst_error_decomp.root");
1109 
1110  c_eff_decomp->Close();
1111  delete c_eff_decomp;
1112  }
1113  else {
1114  for(int ibin = 1; ibin <= signal_nominal->GetXaxis()->GetNbins(); ibin++){
1115 
1116  const double i_nominal = nominal_cut_value_eff->GetBinContent(ibin);
1117 
1118  double i_bkg_syst_uncert = 0;
1119 
1120  for(auto it_syst = shifts_up.begin(); it_syst != shifts_up.end(); ++it_syst){
1121 
1122  double isyst_bkg_syst_uncert=0;
1123 
1124  isyst_bkg_syst_uncert = shifts_up.at(it_syst->first)->GetBinContent(ibin)-i_nominal;
1125 
1126  i_bkg_syst_uncert += isyst_bkg_syst_uncert*isyst_bkg_syst_uncert;
1127 
1128  //Debug syst contributions
1129  if(it_syst->first == "genie")
1130  syst_err_cut_value_eff_genie_up->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
1131  if(it_syst->first == "calib")
1132  syst_err_cut_value_eff_calib_up->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
1133  if(it_syst->first == "light")
1134  syst_err_cut_value_eff_light_up->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
1135  if(it_syst->first == "shape")
1136  syst_err_cut_value_eff_shape->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
1137  if(it_syst->first == "cherenkov")
1138  syst_err_cut_value_eff_cheren->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
1139  } // end for syst loop
1140 
1141  for(auto it_syst = shifts_down.begin(); it_syst != shifts_down.end(); ++it_syst){
1142 
1143  double isyst_bkg_syst_uncert=0;
1144 
1145  isyst_bkg_syst_uncert = shifts_down.at(it_syst->first)->GetBinContent(ibin)-i_nominal;
1146  std::cout << " down " << it_syst->first << " " << isyst_bkg_syst_uncert << std::endl;
1147  i_bkg_syst_uncert += isyst_bkg_syst_uncert*isyst_bkg_syst_uncert;
1148 
1149  //Debug syst contributions
1150  if(it_syst->first == "genie")
1151  syst_err_cut_value_eff_genie_down->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
1152  if(it_syst->first == "calib")
1153  syst_err_cut_value_eff_calib_down->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
1154  if(it_syst->first == "light")
1155  syst_err_cut_value_eff_light_down->SetBinContent(ibin, std::abs(isyst_bkg_syst_uncert));
1156  } // end for syst loop
1157 
1158  syst_err_cut_value_eff_tot->SetBinContent(ibin, sqrt(std::abs(i_bkg_syst_uncert)));
1159  }
1160  TCanvas *c2 = new TCanvas("c2","c2");
1161  syst_err_cut_value_eff_tot->Draw();
1162  outname = "OptimizationPlots/debug/systematic_error_efficiency_vs_lower_cut.pdf";
1163  c2->SaveAs(outname.c_str());
1164  outname_root = "OptimizationPlots/root/systematic_error_efficiency_vs_lower_cut.root";
1165  c2->SaveAs(outname_root.c_str());
1166  c2->Close();
1167  delete c2;
1168 
1169  TCanvas *c_eff_decomp = new TCanvas("c_eff_decomp","c_eff_decomp");
1170  syst_err_cut_value_eff_tot->GetYaxis()->SetTitle("#delta#epsilon");
1171  syst_err_cut_value_eff_tot->SetLineColor(kBlack);
1172  syst_err_cut_value_eff_genie_up->SetLineColor(kRed);
1173  syst_err_cut_value_eff_calib_up->SetLineColor(kGreen);
1174  syst_err_cut_value_eff_light_up->SetLineColor(kOrange);
1175  syst_err_cut_value_eff_genie_down->SetLineColor(kRed+1);
1176  syst_err_cut_value_eff_calib_down->SetLineColor(kGreen+1);
1177  syst_err_cut_value_eff_light_down->SetLineColor(kOrange+1);
1178  syst_err_cut_value_eff_shape->SetLineColor(kMagenta);
1179  syst_err_cut_value_eff_cheren->SetLineColor(kBlue);
1180 
1181  TLegend *leg_eff_decomp = new TLegend(0.65,0.7,0.9,0.9);
1182  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_tot,"Total");
1183  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_genie_up,"Genie up");
1184  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_calib_up,"Calib up");
1185  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_light_up,"Light up");
1186  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_genie_down,"Genie down");
1187  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_calib_down,"Calib down");
1188  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_light_down,"Light down");
1189  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_shape,"Shape");
1190  leg_eff_decomp->AddEntry(syst_err_cut_value_eff_cheren,"Cherenkov");
1191 
1192  syst_err_cut_value_eff_tot->Draw("hist");
1193  syst_err_cut_value_eff_genie_up->Draw("histsame");
1194  syst_err_cut_value_eff_calib_up->Draw("histsame");
1195  syst_err_cut_value_eff_light_up->Draw("histsame");
1196  syst_err_cut_value_eff_genie_down->Draw("histsame");
1197  syst_err_cut_value_eff_calib_down->Draw("histsame");
1198  syst_err_cut_value_eff_light_down->Draw("histsame");
1199  syst_err_cut_value_eff_shape->Draw("histsame");
1200  syst_err_cut_value_eff_cheren->Draw("histsame");
1201  leg_eff_decomp->Draw("same");
1202  c_eff_decomp->SaveAs("OptimizationPlots/efficiency_syst_error_numerator.pdf");
1203  c_eff_decomp->SaveAs("OptimizationPlots/root/efficiency_syst_error_numerator.root");
1204 
1205  TH1* copy_syst_err_cut_value_eff_tot = (TH1F*)syst_err_cut_value_eff_tot->Clone("copy_syst_err_cut_value_eff_tot");
1206  copy_syst_err_cut_value_eff_tot->GetYaxis()->SetTitle("#frac{#delta#epsilon}{#epsilon}");
1207  copy_syst_err_cut_value_eff_tot->GetYaxis()->SetRangeUser(0,1);
1208  copy_syst_err_cut_value_eff_tot->Divide(nominal_cut_value_eff);
1209  syst_err_cut_value_eff_genie_up->Divide(nominal_cut_value_eff);
1210  syst_err_cut_value_eff_calib_up->Divide(nominal_cut_value_eff);
1211  syst_err_cut_value_eff_light_up->Divide(nominal_cut_value_eff);
1212  syst_err_cut_value_eff_genie_down->Divide(nominal_cut_value_eff);
1213  syst_err_cut_value_eff_calib_down->Divide(nominal_cut_value_eff);
1214  syst_err_cut_value_eff_light_down->Divide(nominal_cut_value_eff);
1215  syst_err_cut_value_eff_shape->Divide(nominal_cut_value_eff);
1216  syst_err_cut_value_eff_cheren->Divide(nominal_cut_value_eff);
1217 
1218  copy_syst_err_cut_value_eff_tot->Draw("hist");
1219  syst_err_cut_value_eff_genie_up->Draw("histsame");
1220  syst_err_cut_value_eff_calib_up->Draw("histsame");
1221  syst_err_cut_value_eff_light_up->Draw("histsame");
1222  syst_err_cut_value_eff_genie_down->Draw("histsame");
1223  syst_err_cut_value_eff_calib_down->Draw("histsame");
1224  syst_err_cut_value_eff_light_down->Draw("histsame");
1225  syst_err_cut_value_eff_shape->Draw("histsame");
1226  syst_err_cut_value_eff_cheren->Draw("histsame");
1227  leg_eff_decomp->Draw("same");
1228  c_eff_decomp->SaveAs("OptimizationPlots/efficiency_syst_error_decomp.pdf");
1229  c_eff_decomp->SaveAs("OptimizationPlots/root/efficiency_syst_error_decomp.root");
1230 
1231  c_eff_decomp->Close();
1232  delete c_eff_decomp;
1233  }
1234  return syst_err_cut_value_eff_tot;
1235 }
1236 
1237 TH1F* GetEfficiencyDenominator(Spectrum spectrum_nominal, Spectrum spectrum_truth,
1238  double pot_function)
1239 {
1240  TH1F* nominal = (TH1F*)spectrum_nominal.ToTH1(pot_function);
1241  TH1F* truth = (TH1F*)spectrum_truth.ToTH1(pot_function);
1242 
1243  //Need histograms to hold all intermediate values
1244  //currently there are 5 systematics shifts in both up and down vectors
1245  TH1F* nominal_cut_value_eff = new TH1F("hEfficiency",
1246  TString(";Cut on ") + nominal->GetXaxis()->GetTitle(),
1247  nominal->GetXaxis()->GetNbins(),
1248  nominal->GetXaxis()->GetXmin(),
1249  nominal->GetXaxis()->GetXmax());
1250 
1251  //Integrate To Get Total Number of Events selected for a particular cut value
1252  //Divide By Each samples total integral to convert to efficiency
1253 
1254  double denominator = truth->Integral(1,truth->GetXaxis()->GetNbins());
1255 
1256  for(int ibin = 1; ibin <= nominal->GetXaxis()->GetNbins(); ibin++){
1257  float fNumSelected = nominal->Integral(ibin,nominal->GetXaxis()->GetNbins());
1258  nominal_cut_value_eff->SetBinContent(ibin,fNumSelected/denominator);
1259  }
1260 
1261  TCanvas *c2 = new TCanvas("c2","c2");
1262  nominal_cut_value_eff->Draw();
1263  std::string outname = "OptimizationPlots/debug/efficiency_only.pdf";
1264  c2->SaveAs(outname.c_str());
1265  std::string outname_root = "OptimizationPlots/root/efficiency_only.root";
1266  c2->SaveAs(outname_root.c_str());
1267  c2->Close();
1268  delete c2;
1269 
1270  return nominal_cut_value_eff;
1271 }
1272 
1273 
1274 /////////////////////////////////////////////////////////
1275 /////////////////////////////////////////////////////////
1276 /////////////////////////////////////////////////////////
1277 /////////////////////////////////////////////////////////
1278 /////////////////////////////////////////////////////////
1279 /////////////////////////////////////////////////////////
1280 
1281 
1283 {
1284  // std::string mkdir_name = "mkdir -p OptimizationPlots_" + param + "/root";
1285  // gSystem->Exec(mkdir_name.c_str());
1286 
1287  gSystem->Exec("rm -rf OptimizationPlots");
1288  gSystem->Exec("mkdir -p OptimizationPlots/root");
1289  gSystem->Exec("mkdir -p OptimizationPlots/debug");
1290 
1291  std::string param_sel;
1292  std::string param_obs_sig;
1293  std::string param_obs_bkg;
1294  std::string param_truth;
1295 
1296  if(param=="energy") {
1297  param_sel="RecoPionEnergySel";
1298  param_obs_sig="RecoPionEnergySig";
1299  param_obs_bkg="RecoPionEnergyBkg";
1300  param_truth="RecoPionEnergyTru";
1301  }
1302  else if(param=="pionbdtscore") {
1303  param_sel="BestEventPionScore";
1304  param_obs_sig="BestEventPionScoreSignal";
1305  param_obs_bkg="BestEventPionScoreBackground";
1306  param_truth="BestEventPionScoreTru";
1307  }
1308  else {
1309  std::cout << "ERROR: Parameter to optimize against is not specified or not supported yet!" << std::endl;
1310  return;
1311  }
1312 
1313  TFile *file1 = new TFile("/nova/app/users/projas/dev_10_24_18/ccpiinc_mc_studies-full-calibdown.root", "read");
1314  TFile *file2 = new TFile("/nova/app/users/projas/dev_10_24_18/ccpiinc_mc_studies-full-calibup.root", "read");
1315  TFile *file3 = new TFile("/nova/app/users/projas/dev_10_24_18/ccpiinc_mc_studies-full-ckvproton.root", "read");
1316  TFile *file4 = new TFile("/nova/app/users/projas/dev_10_24_18/ccpiinc_mc_studies-full-negoffset.root", "read");
1317  TFile *file5 = new TFile("/nova/app/users/projas/dev_10_24_18/ccpiinc_mc_studies-full-posoffset.root", "read");
1318  TFile *file6 = new TFile("/nova/app/users/projas/dev_10_24_18/ccpiinc_mc_studies-full-calibshape.root");
1319  TFile *filerw = new TFile("/nova/app/users/projas/dev_10_24_18/ccpiinc_mc_studies-full-geniemv.root", "read");
1320 
1321 
1322  //O == Preselected 1 == Presel && Signal 2 == Presel & Bkgrd
1323 
1324  std::vector<GenieMultiverseSpectra> multiversespec;
1325  std::vector<Spectrum> calibupspect;
1326  std::vector<Spectrum> calibdownspect;
1327  std::vector<Spectrum> lightupspect;
1328  std::vector<Spectrum> lightdownspect;
1329  std::vector<Spectrum> calibshapespect;
1330  std::vector<Spectrum> cherenspect;
1331 
1332 
1333  std::vector<std::unique_ptr<ana::Spectrum>>genie_spects_sel;
1334  std::vector<std::unique_ptr<ana::Spectrum>>genie_spects_obbkg;
1335  std::vector<std::unique_ptr<ana::Spectrum>>genie_spects_obsig;
1336  std::vector<std::unique_ptr<ana::Spectrum>>genie_spects_tru;
1337 
1338 
1339  if(param=="energy") {
1340  for(int k = 0; k < 100; k++){
1341  char name[50];
1342  sprintf(name, "mc_geniemv_recoenergy_%i_sel",k);
1343  genie_spects_sel.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
1344  sprintf(name, "mc_geniemv_recoenergy_%i_sig",k);
1345  genie_spects_obsig.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
1346  sprintf(name, "mc_geniemv_recoenergy_%i_bkg",k);
1347  genie_spects_obbkg.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
1348  sprintf(name, "mc_geniemv_recoenergy_%i_tru",k);
1349  genie_spects_tru.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
1350  }
1351  }
1352 
1353  if(param=="pionbdtscore") {
1354  for(int k = 0; k < 100; k++){
1355  char name[50];
1356  sprintf(name, "mc_geniemv_pionscore_%i_sel",k);
1357  genie_spects_sel.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
1358  sprintf(name, "mc_geniemv_pionscore_%i_sig",k);
1359  genie_spects_obsig.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
1360  sprintf(name, "mc_geniemv_pionscore_%i_bkg",k);
1361  genie_spects_obbkg.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
1362  sprintf(name, "mc_geniemv_pionscore_%i_tru",k);
1363  genie_spects_tru.push_back(Spectrum::LoadFrom(filerw->GetDirectory(name)));
1364  }
1365  }
1366 
1367  GenieMultiverseSpectra sMulti_sel = GenieMultiverseSpectra(100,genie_spects_sel,true);
1368  GenieMultiverseSpectra sMulti_sig = GenieMultiverseSpectra(100,genie_spects_obsig,true);
1369  GenieMultiverseSpectra sMulti_bkg = GenieMultiverseSpectra(100,genie_spects_obbkg,true);
1370  GenieMultiverseSpectra sMulti_truth = GenieMultiverseSpectra(100,genie_spects_tru,true);
1371 
1372  /*
1373  for(int ii = 0; ii < kcSels; ii++){
1374  char multiname[50];
1375  char calibupname[50];
1376  char calibdownname[50];
1377  char lightupname[50];
1378  char lightdownname[50];
1379  char calibshapename[50];
1380  char cherenkovname[50];
1381  sprintf(multiname, "genie_multi_%s_%s", all_vars[32].name.c_str(),
1382  selection[ii].name.c_str());
1383  sprintf(calibupname, "cal_up_%s_%s", all_vars[32].name.c_str(),
1384  selection[ii].name.c_str());
1385  sprintf(calibdownname, "cal_down_%s_%s", all_vars[32].name.c_str(),
1386  selection[ii].name.c_str());
1387  sprintf(lightupname, "light_up_%s_%s", all_vars[32].name.c_str(),
1388  selection[ii].name.c_str());
1389  sprintf(lightdownname, "light_down_%s_%s", all_vars[32].name.c_str(),
1390  selection[ii].name.c_str());
1391  sprintf(calibshapename, "calibshape_%s_%s", all_vars[32].name.c_str(),
1392  selection[ii].name.c_str());
1393  sprintf(cherenkovname, "ctherenkov_%s_%s", all_vars[32].name.c_str(),
1394  selection[ii].name.c_str());
1395  multiversespec.push_back(*GenieMultiverseSpectra::LoadFrom(filerw->GetDirectory(multiname)));
1396  calibupspect.push_back(*Spectrum::LoadFrom(file4->GetDirectory(calibupname)));
1397  calibdownspect.push_back(*Spectrum::LoadFrom(file5->GetDirectory(calibdownname)));
1398  lightupspect.push_back(*Spectrum::LoadFrom(file1->GetDirectory(lightupname)));
1399  lightdownspect.push_back(*Spectrum::LoadFrom(file2->GetDirectory(lightdownname)));
1400  calibshapespect.push_back(*Spectrum::LoadFrom(file6->GetDirectory(calibshapename)));
1401  cherenspect.push_back(*Spectrum::LoadFrom(file3->GetDirectory(cherenkovname)));
1402  }
1403  */
1404 
1405  multiversespec.push_back(sMulti_sel);
1406  multiversespec.push_back(sMulti_sig);
1407  multiversespec.push_back(sMulti_bkg);
1408 
1409  calibupspect.push_back(*Spectrum::LoadFrom(file4->GetDirectory(param_sel.c_str())));
1410  calibdownspect.push_back(*Spectrum::LoadFrom(file5->GetDirectory(param_sel.c_str())));
1411  lightupspect.push_back(*Spectrum::LoadFrom(file1->GetDirectory(param_sel.c_str())));
1412  lightdownspect.push_back(*Spectrum::LoadFrom(file2->GetDirectory(param_sel.c_str())));
1413  calibshapespect.push_back(*Spectrum::LoadFrom(file6->GetDirectory(param_sel.c_str())));
1414  cherenspect.push_back(*Spectrum::LoadFrom(file3->GetDirectory(param_sel.c_str())));
1415 
1416  calibupspect.push_back(*Spectrum::LoadFrom(file4->GetDirectory(param_obs_sig.c_str())));
1417  calibdownspect.push_back(*Spectrum::LoadFrom(file5->GetDirectory(param_obs_sig.c_str())));
1418  lightupspect.push_back(*Spectrum::LoadFrom(file1->GetDirectory(param_obs_sig.c_str())));
1419  lightdownspect.push_back(*Spectrum::LoadFrom(file2->GetDirectory(param_obs_sig.c_str())));
1420  calibshapespect.push_back(*Spectrum::LoadFrom(file6->GetDirectory(param_obs_sig.c_str())));
1421  cherenspect.push_back(*Spectrum::LoadFrom(file3->GetDirectory(param_obs_sig.c_str())));
1422 
1423  calibupspect.push_back(*Spectrum::LoadFrom(file4->GetDirectory(param_obs_bkg.c_str())));
1424  calibdownspect.push_back(*Spectrum::LoadFrom(file5->GetDirectory(param_obs_bkg.c_str())));
1425  lightupspect.push_back(*Spectrum::LoadFrom(file1->GetDirectory(param_obs_bkg.c_str())));
1426  lightdownspect.push_back(*Spectrum::LoadFrom(file2->GetDirectory(param_obs_bkg.c_str())));
1427  calibshapespect.push_back(*Spectrum::LoadFrom(file6->GetDirectory(param_obs_bkg.c_str())));
1428  cherenspect.push_back(*Spectrum::LoadFrom(file3->GetDirectory(param_obs_bkg.c_str())));
1429 
1430  auto truth_nominal = *sMulti_truth.Nominal();
1431  auto truth_genieup = *sMulti_truth.UpperSigma();
1432  auto truth_geniedwn = *sMulti_truth.LowerSigma();
1433  auto truth_calup = *Spectrum::LoadFrom(file4->GetDirectory(param_truth.c_str()));
1434  auto truth_caldwn = *Spectrum::LoadFrom(file5->GetDirectory(param_truth.c_str()));
1435  auto truth_lightup = *Spectrum::LoadFrom(file1->GetDirectory(param_truth.c_str()));
1436  auto truth_lightdwn = *Spectrum::LoadFrom(file2->GetDirectory(param_truth.c_str()));
1437  auto truth_calibshape = *Spectrum::LoadFrom(file6->GetDirectory(param_truth.c_str()));
1438  auto truth_cherenkov = *Spectrum::LoadFrom(file3->GetDirectory(param_truth.c_str()));
1439 
1440  TCanvas * c100 = new TCanvas("c100","c100");
1441  auto signal_multiverse = multiversespec[1].Nominal()->ToTH1(pot);
1442  auto background_multiverse = multiversespec[2].Nominal()->ToTH1(pot);
1443  background_multiverse->Draw();
1444  signal_multiverse->Draw("hist same");
1445  c100->SaveAs("OptimizationPlots/debug/nominal_distributions.pdf");
1446  c100->SaveAs("OptimizationPlots/root/nominal_distributions.root");
1447  c100->Close();
1448  delete c100;
1449 
1450  std::cout << truth_nominal.POT() << std::endl;
1451 
1452  std::map<std::string,Spectrum> systs_up_bkg;
1453  std::map<std::string,Spectrum> systs_dwn_bkg;
1454  std::map<std::string,Spectrum> systs_up_signal;
1455  std::map<std::string,Spectrum> systs_dwn_signal;
1456  std::map<std::string,Spectrum> systs_up_truth;
1457  std::map<std::string,Spectrum> systs_dwn_truth;
1458 
1459  //Let's plot the error bands first
1460  int makePlotNumber = 0;
1461  std::string makePlotString = "presel_";
1462  //bool shouldPlot = false;
1463  bool shouldPlot = true;
1464  std::string outname;
1465  std::string outname_root;
1466 
1467  //First Up GENIE
1468  systs_up_truth.insert(std::make_pair("genie",truth_genieup));
1469  systs_dwn_truth.insert(std::make_pair("genie",truth_geniedwn));
1470  systs_up_signal.insert(std::make_pair("genie",*multiversespec[1].UpperSigma()));
1471  systs_dwn_signal.insert(std::make_pair("genie",*multiversespec[1].LowerSigma()));
1472  systs_up_bkg.insert(std::make_pair("genie",*multiversespec[2].UpperSigma()));
1473  systs_dwn_bkg.insert(std::make_pair("genie",*multiversespec[2].LowerSigma()));
1474 
1475  if(shouldPlot){
1476  for(uint i=0; i<3; i++){
1477  if (i==0) makePlotString = "presel_";
1478  if (i==1) makePlotString = "preselsig_";
1479  if (i==2) makePlotString = "preselbkg_";
1480  std::map<std::string,Spectrum> genie_up;
1481  genie_up.insert(std::make_pair("genie",*multiversespec[i].UpperSigma()));
1482  std::map<std::string,Spectrum> genie_dwn;
1483  genie_dwn.insert(std::make_pair("genie",*multiversespec[i].LowerSigma()));
1484 
1485  new TCanvas();
1486  // PlotWithSystErrorBand(*(multiversespec[i].Nominal()),
1487  // genie_up, genie_dwn,pot);
1488  outname = "OptimizationPlots/debug/" + makePlotString + "genie_errorband.pdf";
1489  gPad->SaveAs(outname.c_str());
1490  outname_root = "OptimizationPlots/root/" + makePlotString + "genie_errorband.root";
1491  gPad->SaveAs(outname_root.c_str());
1492  }
1493  }
1494 
1495  //Next Calibration Normalization
1496  systs_up_truth.insert(std::make_pair("calib",truth_calup));
1497  systs_dwn_truth.insert(std::make_pair("calib",truth_caldwn));
1498  systs_up_signal.insert(std::make_pair("calib",calibupspect[1]));
1499  systs_dwn_signal.insert(std::make_pair("calib",calibdownspect[1]));
1500  systs_dwn_bkg.insert(std::make_pair("calib",calibdownspect[2]));
1501  systs_up_bkg.insert(std::make_pair("calib",calibupspect[2]));
1502 
1503  if(shouldPlot){
1504  for(uint i=0; i<3; i++){
1505  if (i==0) makePlotString = "presel_";
1506  if (i==1) makePlotString = "preselsig_";
1507  if (i==2) makePlotString = "preselbkg_";
1508 
1509  std::map<std::string,Spectrum> cal_up;
1510  cal_up.insert(std::make_pair("calib", calibupspect[i]));
1511  std::map<std::string,Spectrum> cal_dwn;
1512  cal_dwn.insert(std::make_pair("calib", calibdownspect[i]));
1513 
1514  new TCanvas();
1515  // PlotWithSystErrorBand(*(multiversespec[i].Nominal()),
1516  // cal_up, cal_dwn,
1517  // pot);
1518  outname = "OptimizationPlots/debug/" + makePlotString + "calibration_errorband.pdf";
1519  gPad->SaveAs(outname.c_str());
1520  outname_root = "OptimizationPlots/root/" + makePlotString + "calibration_errorband.root";
1521  gPad->SaveAs(outname_root.c_str());
1522  new TCanvas();
1523  }
1524  }
1525 
1526  //Next Light Levels
1527 
1528  systs_up_truth.insert(std::make_pair("light",truth_lightup));
1529  systs_dwn_truth.insert(std::make_pair("light",truth_lightdwn));
1530  systs_up_signal.insert(std::make_pair("light",lightupspect[1]));
1531  systs_dwn_signal.insert(std::make_pair("light",lightdownspect[1]));
1532  systs_up_bkg.insert(std::make_pair("light",lightupspect[2]));
1533  systs_dwn_bkg.insert(std::make_pair("light",lightdownspect[2]));
1534 
1535  if(shouldPlot){
1536  for(uint i=0; i<3; i++){
1537  if (i==0) makePlotString = "presel_";
1538  if (i==1) makePlotString = "preselsig_";
1539  if (i==2) makePlotString = "preselbkg_";
1540 
1541  std::map<std::string,Spectrum> light_up;
1542  light_up.insert(std::make_pair("light", lightupspect[i]));
1543  std::map<std::string,Spectrum> light_dwn;
1544  light_dwn.insert(std::make_pair("light", lightdownspect[i]));
1545 
1546  new TCanvas();
1547  // PlotWithSystErrorBand(*(multiversespec[i].Nominal()),
1548  // light_up, light_dwn, pot);
1549  outname = "OptimizationPlots/debug/" + makePlotString + "light_level_errorband.pdf";
1550  gPad->SaveAs(outname.c_str());
1551  outname_root = "OptimizationPlots/root/" + makePlotString + "light_level_errorband.root";
1552  gPad->SaveAs(outname_root.c_str());
1553  }
1554  }
1555 
1556  //Next CalibrationShape
1557  systs_up_truth.insert(std::make_pair("shape",truth_calibshape));
1558  systs_up_signal.insert(std::make_pair("shape",calibshapespect[1]));
1559  systs_up_bkg.insert(std::make_pair("shape",calibshapespect[2]));
1560 
1561  if(shouldPlot){
1562  for(uint i=0; i<3; i++){
1563  if (i==0) makePlotString = "presel_";
1564  if (i==1) makePlotString = "preselsig_";
1565  if (i==2) makePlotString = "preselbkg_";
1566 
1567  std::map<std::string,Spectrum> shape;
1568  shape.insert(std::make_pair("shape", calibshapespect[i]));
1569 
1570  new TCanvas();
1571  // PlotWithSystErrorBand(*(multiversespec[i].Nominal()),
1572  // shape, shape, pot);
1573  outname = "OptimizationPlots/debug/" + makePlotString + "calibration_shape_errorband.pdf";
1574  gPad->SaveAs(outname.c_str());
1575  outname_root = "OptimizationPlots/root/" + makePlotString + "calibration_shape_errorband.root";
1576  gPad->SaveAs(outname_root.c_str());
1577  }
1578  }
1579 
1580  //Next Cherenekov
1581  systs_up_truth.insert(std::make_pair("cherenkov",truth_cherenkov));
1582  systs_up_signal.insert(std::make_pair("cherenkov",cherenspect[1]));
1583  systs_up_bkg.insert(std::make_pair("cherenkov",cherenspect[2]));
1584 
1585  if(shouldPlot){
1586  for(uint i=0; i<3; i++){
1587  if (i==0) makePlotString = "presel_";
1588  if (i==1) makePlotString = "preselsig_";
1589  if (i==2) makePlotString = "preselbkg_";
1590  std::map<std::string,Spectrum> cheren;
1591  cheren.insert(std::make_pair("cherenkov", calibshapespect[i]));
1592 
1593  new TCanvas();
1594  // PlotWithSystErrorBand(*(smultiversespec[i].Nominal()),
1595  // cheren, cheren, pot);
1596  outname = "OptimizationPlots/debug/" + makePlotString + "cherenekov_errorband.pdf";
1597  gPad->SaveAs(outname.c_str());
1598  outname_root = "OptimizationPlots/root/" + makePlotString + "cherenekov_errorband.root";
1599  gPad->SaveAs(outname_root.c_str());
1600  }
1601  }
1602 
1603  TCanvas *c_sel = new TCanvas("c_sel","c_sel");
1604  TH1F* sel_nominal = (TH1F*)multiversespec[0].Nominal()->ToTH1(pot);
1605  TH1F* sel_genie_up = (TH1F*)multiversespec[0].UpperSigma()->ToTH1(pot);
1606  TH1F* sel_genie_down = (TH1F*)multiversespec[0].LowerSigma()->ToTH1(pot);
1607  TH1F* sel_cal_up = (TH1F*)calibupspect[0].ToTH1(pot);
1608  TH1F* sel_cal_down = (TH1F*)calibdownspect[0].ToTH1(pot);
1609  TH1F* sel_light_up = (TH1F*)lightupspect[0].ToTH1(pot);
1610  TH1F* sel_light_down = (TH1F*)lightdownspect[0].ToTH1(pot);
1611  TH1F* sel_shape = (TH1F*)calibshapespect[0].ToTH1(pot);
1612  TH1F* sel_cherenkov = (TH1F*)cherenspect[0].ToTH1(pot);
1613 
1614  sel_nominal->SetLineColor(kBlack);
1615  sel_genie_up->SetLineColor(kRed);
1616  sel_genie_down->SetLineColor(kRed+1);
1617  sel_cal_up->SetLineColor(kGreen);
1618  sel_cal_down->SetLineColor(kGreen+1);
1619  sel_light_up->SetLineColor(kOrange);
1620  sel_light_down->SetLineColor(kOrange+1);
1621  sel_shape->SetLineColor(kMagenta);
1622  sel_cherenkov->SetLineColor(kBlue);
1623 
1624  TLegend *leg = new TLegend(0.5,0.15,0.8,0.35);
1625  leg->AddEntry(sel_nominal,"Nominal");
1626  leg->AddEntry(sel_genie_up,"Genie up");
1627  leg->AddEntry(sel_genie_down,"Genie down");
1628  leg->AddEntry(sel_cal_up,"Cal up");
1629  leg->AddEntry(sel_cal_down,"Cal down");
1630  leg->AddEntry(sel_light_up,"Light up");
1631  leg->AddEntry(sel_light_down,"Light down");
1632  leg->AddEntry(sel_shape,"Shape");
1633  leg->AddEntry(sel_cherenkov,"Cherenkov");
1634 
1635  sel_genie_up->Draw("hist");
1636  sel_genie_down->Draw("histsame");
1637  sel_cal_up->Draw("histsame");
1638  sel_cal_down->Draw("histsame");
1639  sel_light_up->Draw("histsame");
1640  sel_light_down->Draw("histsame");
1641  sel_shape->Draw("histsame");
1642  sel_cherenkov->Draw("histsame");
1643  sel_nominal->Draw("histsame");
1644 
1645  leg->Draw("same");
1646  c_sel->SaveAs("OptimizationPlots/selected_comparison.pdf");
1647  c_sel->SaveAs("OptimizationPlots/root/selected_comparison.root");
1648  c_sel->Close();
1649  delete c_sel;
1650 
1651  TCanvas *c_sig = new TCanvas("c_sig","c_sig");
1652  TH1F* sig_nominal = (TH1F*)multiversespec[1].Nominal()->ToTH1(pot);
1653  TH1F* sig_genie_up = (TH1F*)multiversespec[1].UpperSigma()->ToTH1(pot);
1654  TH1F* sig_genie_down = (TH1F*)multiversespec[1].LowerSigma()->ToTH1(pot);
1655  TH1F* sig_cal_up = (TH1F*)calibupspect[1].ToTH1(pot);
1656  TH1F* sig_cal_down = (TH1F*)calibdownspect[1].ToTH1(pot);
1657  TH1F* sig_light_up = (TH1F*)lightupspect[1].ToTH1(pot);
1658  TH1F* sig_light_down = (TH1F*)lightdownspect[1].ToTH1(pot);
1659  TH1F* sig_shape = (TH1F*)calibshapespect[1].ToTH1(pot);
1660  TH1F* sig_cherenkov = (TH1F*)cherenspect[1].ToTH1(pot);
1661 
1662  sig_nominal->SetLineColor(kBlack);
1663  sig_genie_up->SetLineColor(kRed);
1664  sig_genie_down->SetLineColor(kRed+1);
1665  sig_cal_up->SetLineColor(kGreen);
1666  sig_cal_down->SetLineColor(kGreen+1);
1667  sig_light_up->SetLineColor(kOrange);
1668  sig_light_down->SetLineColor(kOrange+1);
1669  sig_shape->SetLineColor(kMagenta);
1670  sig_cherenkov->SetLineColor(kBlue);
1671 
1672  sig_genie_up->Draw("hist");
1673  sig_genie_down->Draw("histsame");
1674  sig_cal_up->Draw("histsame");
1675  sig_cal_down->Draw("histsame");
1676  sig_light_up->Draw("histsame");
1677  sig_light_down->Draw("histsame");
1678  sig_shape->Draw("histsame");
1679  sig_cherenkov->Draw("histsame");
1680  sig_nominal->Draw("histsame");
1681 
1682  leg->Draw("same");
1683  c_sig->SaveAs("OptimizationPlots/signal_comparison.pdf");
1684  c_sig->SaveAs("OptimizationPlots/root/signal_comparison.root");
1685  c_sig->Close();
1686  delete c_sig;
1687 
1688  TCanvas *c_bkg = new TCanvas("c_bkg","c_bkg");
1689  TH1F* bkg_nominal = (TH1F*)multiversespec[2].Nominal()->ToTH1(pot);
1690  TH1F* bkg_genie_up = (TH1F*)multiversespec[2].UpperSigma()->ToTH1(pot);
1691  TH1F* bkg_genie_down = (TH1F*)multiversespec[2].LowerSigma()->ToTH1(pot);
1692  TH1F* bkg_cal_up = (TH1F*)calibupspect[2].ToTH1(pot);
1693  TH1F* bkg_cal_down = (TH1F*)calibdownspect[2].ToTH1(pot);
1694  TH1F* bkg_light_up = (TH1F*)lightupspect[2].ToTH1(pot);
1695  TH1F* bkg_light_down = (TH1F*)lightdownspect[2].ToTH1(pot);
1696  TH1F* bkg_shape = (TH1F*)calibshapespect[2].ToTH1(pot);
1697  TH1F* bkg_cherenkov = (TH1F*)cherenspect[2].ToTH1(pot);
1698 
1699  bkg_nominal->SetLineColor(kBlack);
1700  bkg_genie_up->SetLineColor(kRed);
1701  bkg_genie_down->SetLineColor(kRed+1);
1702  bkg_cal_up->SetLineColor(kGreen);
1703  bkg_cal_down->SetLineColor(kGreen+1);
1704  bkg_light_up->SetLineColor(kOrange);
1705  bkg_light_down->SetLineColor(kOrange+1);
1706  bkg_shape->SetLineColor(kMagenta);
1707  bkg_cherenkov->SetLineColor(kBlue);
1708 
1709  bkg_genie_up->Draw("hist");
1710  bkg_genie_down->Draw("histsame");
1711  bkg_cal_up->Draw("histsame");
1712  bkg_cal_down->Draw("histsame");
1713  bkg_light_up->Draw("histsame");
1714  bkg_light_down->Draw("histsame");
1715  bkg_shape->Draw("histsame");
1716  bkg_cherenkov->Draw("histsame");
1717  bkg_nominal->Draw("histsame");
1718 
1719  leg->Draw("same");
1720  c_bkg->SaveAs("OptimizationPlots/background_comparison.pdf");
1721  c_bkg->SaveAs("OptimizationPlots/root/background_comparison.root");
1722  c_bkg->Close();
1723  delete c_bkg;
1724 
1725  TH1F* selected_statistical = GetSelectedStatisticalUncertainty(*(multiversespec[0].Nominal()), pot);
1726  TH1F* background_statistical = GetBackgroundStatisticalUncertainty(*(multiversespec[2].Nominal()), mc_pot);
1727  TH1F* denominator = GetDenominator(*(multiversespec[0].Nominal()),
1728  *(multiversespec[2].Nominal()),
1729  pot);
1730  TH1F* background_systematic = GetBackgroundSystematicUncertainty(*(multiversespec[2].Nominal()),denominator, systs_up_bkg, systs_dwn_bkg, pot);
1731  TH1F* efficiency_systematic = GetEfficiencySystematicUncertainty(*(multiversespec[1].Nominal()), truth_nominal, systs_up_signal, systs_dwn_signal, systs_up_truth, systs_dwn_truth, pot);
1732  TH1F* efficiency_denominator = GetEfficiencyDenominator(*(multiversespec[1].Nominal()), systs_up_truth.at("genie"), pot);
1733  // TH1F* efficiency_denominator = GetEfficiencyDenominator(*(multiversespec[1].Nominal()), truth_nominal, pot);
1734 
1735  TCanvas *c1 = new TCanvas("c1","c1");
1736  c1->SetLeftMargin(0.15);
1737  selected_statistical->GetYaxis()->SetTitle("#frac{#sqrt{N_{sel}}}{(N_{sel} - N_{bkg})}");
1738  selected_statistical->GetXaxis()->SetTitle(TString("Cut on ") + multiversespec[1].Nominal()->ToTH1(pot)->GetXaxis()->GetTitle());
1739  selected_statistical->Divide(denominator);
1740  selected_statistical->Draw();
1741  outname = "OptimizationPlots/stat_uncert_selected_contribution.pdf";
1742  Simulation();
1743  c1->SaveAs(outname.c_str());
1744  outname_root = "OptimizationPlots/root/stat_uncert_selected_contribution.root";
1745  c1->SaveAs(outname_root.c_str());
1746  c1->Close();
1747  delete c1;
1748 
1749  TCanvas *c2 = new TCanvas("c2","c2");
1750  c2->SetLeftMargin(0.15);
1751  background_statistical->GetYaxis()->SetTitle("#frac{#sqrt{N_{bkg}}}{(N_{sel} - N_{bkg})}");
1752  background_statistical->GetXaxis()->SetTitle(TString("Cut on ") + multiversespec[1].Nominal()->ToTH1(pot)->GetXaxis()->GetTitle());
1753  background_statistical->Divide(denominator);
1754  background_statistical->Draw();
1755  outname = "OptimizationPlots/stat_uncert_background_contribution.pdf";
1756  Simulation();
1757  c2->SaveAs(outname.c_str());
1758  outname_root = "OptimizationPlots/root/stat_uncert_background_contribution.root";
1759  c2->SaveAs(outname_root.c_str());
1760  c2->Close();
1761  delete c2;
1762 
1763  TCanvas *c3 = new TCanvas("c3","c3");
1764  c3->SetLeftMargin(0.15);
1765  background_systematic->GetYaxis()->SetTitle("#frac{#delta N_{bkg}}{(N_{sel} - N_{bkg})}");
1766  background_systematic->GetXaxis()->SetTitle(TString("Cut on ") + multiversespec[1].Nominal()->ToTH1(pot)->GetXaxis()->GetTitle());
1767  background_systematic->Divide(denominator);
1768  background_systematic->SetMaximum(1);
1769  background_systematic->Draw();
1770  Simulation();
1771  outname = "OptimizationPlots/syst_uncert_background_contribution.pdf";
1772  c3->SaveAs(outname.c_str());
1773  outname_root = "OptimizationPlots/root/syst_uncert_background_contribution.root";
1774  c3->SaveAs(outname_root.c_str());
1775  c3->Close();
1776  delete c3;
1777 
1778  TCanvas *c4 = new TCanvas("c4","c4");
1779  c4->SetLeftMargin(0.15);
1780  efficiency_systematic->GetYaxis()->SetTitle("#frac{#delta #epsilon}{#epsilon}");
1781  efficiency_systematic->GetXaxis()->SetTitle(TString("Cut on ") + multiversespec[1].Nominal()->ToTH1(pot)->GetXaxis()->GetTitle());
1782  efficiency_systematic->Divide(efficiency_denominator);
1783  efficiency_systematic->SetMaximum(1);
1784  efficiency_systematic->Draw();
1785  Simulation();
1786  outname = "OptimizationPlots/syst_uncert_efficiency_contribution.pdf";
1787  c4->SaveAs(outname.c_str());
1788  outname_root = "OptimizationPlots/root/syst_uncert_efficiency_contribution.root";
1789  c4->SaveAs(outname_root.c_str());
1790  c4->Close();
1791  delete c4;
1792 
1793  TH1F* dsigma = new TH1F("dsigma",
1794  TString(";Cut on ") + multiversespec[1].Nominal()->ToTH1(pot)->GetXaxis()->GetTitle() + ";#frac{#delta #sigma}{#sigma}",
1795  efficiency_systematic->GetXaxis()->GetNbins(),
1796  efficiency_systematic->GetXaxis()->GetXmin(),
1797  efficiency_systematic->GetXaxis()->GetXmax());
1798 
1799  for( int ibin = 1; ibin <= efficiency_systematic->GetXaxis()->GetNbins(); ibin++){
1800  double first = selected_statistical->GetBinContent(ibin);
1801  double second = background_statistical->GetBinContent(ibin);
1802  double third = background_systematic->GetBinContent(ibin);
1803  double fourth = efficiency_systematic->GetBinContent(ibin);
1804 
1805  double final = first*first + second*second + third*third + fourth*fourth;
1806 
1807  final = sqrt(final);
1808 
1809  dsigma->SetBinContent(ibin,final);
1810  }
1811 
1812  TCanvas *c5 = new TCanvas("c5","c5");
1813  c5->SetLeftMargin(0.15);
1814  dsigma->SetMaximum(1);
1815  dsigma->Draw();
1816  Simulation();
1817  outname = "OptimizationPlots/xsec_total_uncertainty.pdf";
1818  c5->SaveAs(outname.c_str());
1819  outname_root = "OptimizationPlots/root/xsec_total_uncertainty.root";
1820  c5->SaveAs(outname_root.c_str());
1821  c5->Close();
1822  delete c5;
1823 
1824  file1->Close();
1825  file2->Close();
1826  file3->Close();
1827  file4->Close();
1828  file5->Close();
1829  file6->Close();
1830  filerw->Close();
1831 
1832 }
void Simulation()
Definition: tools.h:16
const XML_Char * name
Definition: expat.h:151
enum BeamMode kOrange
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
const double pot
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1F * GetBackgroundSystematicUncertainty(Spectrum spectrum_nominal, TH1F *denominator, std::map< std::string, Spectrum > ShiftUp, std::map< std::string, Spectrum > ShiftDown, double pot_function)
const double mc_pot
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 * GetEfficiencyDenominator(Spectrum spectrum_nominal, Spectrum spectrum_truth, 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()
const bool up_down_combined
TH1F * GetSelectedStatisticalUncertainty(Spectrum spectrum_nominal, double pot_function)
TH1F * GetBackgroundStatisticalUncertainty(Spectrum spectrum_nominal, double pot_function)
OStream cout
Definition: OStream.cxx:6
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal) const
TGraphAsymmErrors * MyPlotWithSystErrorBand(TH1F *&nom, std::map< std::string, TH1F * > &ups, std::map< std::string, TH1F * > &dns, int col, int errCol, float headroom, bool newaxis)
c1
Definition: demo5.py:24
const Spectrum * LowerSigma(BandOptions opt=kBandFromNominal) const
const Color_t kTotalMCColor
Definition: Style.h:16
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)
TH1F * GetDenominator(Spectrum spectrum_selected, Spectrum spectrum_background, double pot_function)
enum BeamMode kGreen
enum BeamMode kBlue
void xsec_tot_uncert_optimization(std::string param="pionbdtscore")
const Spectrum * Nominal() const
return the nominal universe
enum BeamMode string
unsigned int uint