make_muonid_opt.C
Go to the documentation of this file.
1 
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include "Utilities/rootlogon.C"
8 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Analysis/Plots.h"
11 #include "CAFAna/Analysis/Style.h"
13 #include "CAFAna/Core/Ratio.h"
15 
16 
17 #include "TH1.h"
18 #include "TH2.h"
19 #include "TCanvas.h"
20 #include "TLegend.h"
21 
22 #include "math.h"
23 
24 using namespace ana;
25 
26 
27 const double pot = 8.09e20; // POT
28 const double potmc = 3.54e21; // POT
29 // const double potmc = 8.09e20; // POT
30 const unsigned int N_UNIVERSES = 1000;
31 
32 
34 
35 void FindMinimum(TH1*& h, const int minx, const int maxx)
36 {
37  int nx = h->GetXaxis()->GetNbins();
38 
39  double min = 100;
40  int minxbin = 0;
41 
42  for(int ix = minx; ix <= maxx; ix++){
43  if( h->GetBinContent(ix) < min){
44  min = h->GetBinContent(ix);
45  minxbin = ix;
46  }
47  }// end loop over x bins
48 
49  std::cout<<"Best Muon ID cut: "<<h->GetXaxis()->GetBinLowEdge(minxbin)
50  <<" Uncertainty: "<< h->GetBinContent(minxbin)<<"\n";
51 }
52 
53 //------------------------------------------------------------
54 //------------------------------------------------------------
55 //------------------------------------------------------------
56 
57 void PrintPlots(TH1*& h,
58  TString name)
59 {
60  TCanvas *c = new TCanvas();
61  c->SetLeftMargin(0.18);
62  h->GetXaxis()->CenterTitle();
63  h->GetYaxis()->CenterTitle();
64  h->GetYaxis()->SetTitleOffset(1.2);
65  // h->GetYaxis()->SetTitleSize(0.02);
66  h->SetLineColor(kRed);
67  h->SetLineWidth(2);
68  h->Draw("hist ][");
69  //c->SetLogz(1);
70  c->SaveAs(string(outputPlotDir + "/" + name).c_str());
71 
72  // Store the plot in the output root file, making sure to cut off the ".pdf" extension
73  h->Write(string(string(name).substr(0, string(name).find_last_of("."))).c_str());
74 }
75 //------------------------------------------------------------
76 //------------------------------------------------------------
77 //------------------------------------------------------------
78 
80 {
81  TLatex* prelim = new TLatex(.9, .95, "NOvA ND Simulation");
82  prelim->SetTextColor(kGray+1);
83  prelim->SetNDC();
84  prelim->SetTextSize(2/30.);
85  prelim->SetTextAlign(32);
86  prelim->Draw();
87 }
88 
89 //------------------------------------------------------------
90 //------------------------------------------------------------
91 //------------------------------------------------------------
92 
93 void ComputeErrors(std::vector<TH1*>& bg,
94  std::vector<TH1*>& sig,
95  std::vector<TH1*>& eff,
96  std::vector<TH1*>& bg_nom,
97  TString xylabel,
98  TString figname)
99 {
100 
101  TH1* nom = sig[0];
102  TH1* bgnom = bg[0];
103  TH1* effnom = eff[0];
104  TH1* bgstat = bg_nom[0];
105 
106  int nx = nom->GetNbinsX();
107 
108  //int minx = nom->GetXaxis()->FindBin(-1);
109  //int maxx = nom->GetXaxis()->FindBin(1);
110 
111  // std::cout <<" minx "<< minx << " maxx " << maxx << " nbinx " << nx <<std::endl;
112 
113  TH1* hsig_stat_uncert = (TH1*) nom->Clone("hsig_stat_uncert");
114  hsig_stat_uncert->SetTitle("Relative Stat. Uncertainty on Sel" + xylabel + "#frac{#sqrt{N_{sel}}}{N_{sel}-N_{bkg}}");
115  // hsig_stat_uncert->GetYaxis()->SetTitleSize(0.04);
116 
117  TH1* hbg_stat_uncert = (TH1*) nom->Clone("hbg_stat_uncert");
118  hbg_stat_uncert->SetTitle("Relative Stat. Uncertainty on Bkg" + xylabel + "#frac{#sqrt{N_{bkg}}}{N_{sel}-N_{bkg}}");
119  // hbg_stat_uncert->GetYaxis()->SetTitleSize(0.04);
120 
121  TH1* hbg_syst_uncert = (TH1*) nom->Clone("hbg_syst_uncert");
122  hbg_syst_uncert->SetTitle("Relative Syst. Uncertainty on Bkg" + xylabel + "#frac{#deltaN^{syst}_{bkg}}{N_{sel}-N_{bkg}}");
123  // hbg_syst_uncert->GetYaxis()->SetTitleSize(0.04);
124 
125  TH1* heff_syst_uncert = (TH1*) nom->Clone("heff_syst_uncert");
126  heff_syst_uncert->SetTitle("Relative Syst. Uncertainty on Eff" + xylabel + "#frac{#delta#varepsilon}{#varepsilon}");
127 
128  TH1* hxsec_uncert = (TH1*) nom->Clone("hxsec_uncert");
129  hxsec_uncert->SetTitle("Relative Uncertainty on Cross-section" + xylabel + "#frac{#delta#sigma}{#sigma}");
130 
131  TH1* heff_uncert = (TH1*) nom->Clone("heff_uncert");
132  heff_uncert->SetTitle("Uncertainty on Efficiency" + xylabel + "#delta#varepsilon");
133 
134  //TH1* heff = (TH1*) effnom->Clone("heff");
135  //heff->SetTitle("Efficiency" + xylabel + "Efficiency");
136 
137  //TH1* heff = (TH1*) nom->Clone("heff");
138  // heff->SetTitle("Efficiency" + xylabel + "#varepsilon");
139  // hxsec_uncert->SetMaximum(0.3);
140  //hsig_stat_uncert->SetMaximum(0.05);
141  //hbg_stat_uncert->SetMaximum(0.05);
142  //hbg_syst_uncert->SetMaximum(0.1);
143  //heff_syst_uncert->SetMaximum(0.2);
144  // }
145  int nsysts = (sig.size() - 1) / 2;
146 
147  // Set everything to 0
148  for(int ix = 1; ix <= nx; ++ix){
149  hbg_syst_uncert->SetBinContent( ix, 0);
150  hbg_stat_uncert->SetBinContent( ix, 0);
151  heff_syst_uncert->SetBinContent( ix, 0);
152  hxsec_uncert->SetBinContent( ix, 0);
153  heff_uncert->SetBinContent( ix, 0);
154  hsig_stat_uncert->SetBinContent( ix, 0);
155  }
156 
157  for(int ix = 1; ix <= nx; ++ix){
158 
159  const double nsel = nom->GetBinContent(ix);
160  const double nbg = bgnom->GetBinContent(ix);
161  const double neff = effnom->GetBinContent(ix);
162  double nbgstat = bgstat->GetBinContent(ix);
163  double bg_error_syst = 0;
164  double eff_error_syst = 0;
165 
166 
167  // compute systematic errors in each bin for
168  // background, signal, and efficiency
169 
170  for(int isyst = 1; isyst <= nsysts; ++isyst){
171  // std :: cout << isyst << "......." << nsysts + isyst << std::endl;
172  double ieff_down = eff[isyst]->GetBinContent(ix);
173  double ieff_up = eff[nsysts + isyst]->GetBinContent(ix);
174  double ibg_down = bg[isyst]->GetBinContent(ix);
175  double ibg_up = bg[nsysts + isyst]->GetBinContent(ix);
176 
177  double bg_error_isyst = std::max( std::abs(ibg_up - nbg), std::abs(ibg_down - nbg) );
178  double eff_error_isyst = std::max( std::abs(ieff_up - neff), std::abs(ieff_down - neff) );
179 
180  bg_error_syst += std::pow(bg_error_isyst, 2);
181  eff_error_syst += std::pow(eff_error_isyst, 2);
182  } // end for isyst
183 
184  double bg_error_stat = nbgstat;
185  double sel_error_stat = nsel;
186  double totbg_error = sel_error_stat + bg_error_stat + bg_error_syst;
187 
188 
189  // Don't divide by zero.
190  double bg_error_term = 0;
191  if(std::pow( nsel - nbg, 2) != 0)
192  bg_error_term = totbg_error / std::pow( nsel - nbg, 2);
193 
194 
195  double eff_error_term = 0;
196  if(std::pow( neff, 2) != 0)
197  eff_error_term = eff_error_syst / std::pow( neff, 2);
198 
199 
200  double xsec_uncert = eff_error_term + bg_error_term;
201 
202  if (nsel-nbg != 0){
203  hbg_syst_uncert->SetBinContent( ix, std::sqrt(bg_error_syst) / (nsel - nbg));
204  hbg_stat_uncert->SetBinContent( ix, std::sqrt(bg_error_stat) / (nsel - nbg));
205  hsig_stat_uncert->SetBinContent( ix, std::sqrt(sel_error_stat) / (nsel - nbg));
206  }
207 
208  heff_syst_uncert->SetBinContent( ix, std::sqrt(eff_error_term));
209 
210  hxsec_uncert->SetBinContent( ix, std::sqrt(xsec_uncert));
211 
212  heff_uncert->SetBinContent( ix, std::sqrt(eff_error_syst));
213 
214  //heff->SetBinContent( ix, neff);
215  /*
216  std::cout<<" vertex "<<nom->GetXaxis()->GetBinCenter(ix)
217  <<" nsel "<<nsel<<" nbkg "<< nbg << " eff " << neff
218  <<" syst bkg " << std::sqrt(bg_error_syst)
219  << " nsel - nbkg " << nsel - nbg
220  << "bg_syst_uncert " << std::sqrt(bg_error_syst) / (nsel - nbg)
221  <<" xsec error "<<std::sqrt(xsec_uncert)<<std::endl;
222  */
223  } // end loop for x bins
224 
225  // NDSimulation();
226  PrintPlots(hbg_syst_uncert, "bg_syst_uncert_"+figname+".pdf");
227  PrintPlots(hbg_stat_uncert, "bg_stat_uncert_"+figname+".pdf");
228  PrintPlots(heff_syst_uncert, "eff_syst_uncert_"+figname+".pdf");
229  PrintPlots(hxsec_uncert, "xsec_uncert_"+figname+".pdf");
230  PrintPlots(heff_uncert, "eff_uncert_"+figname+".pdf");
231  PrintPlots(hsig_stat_uncert, "sig_stat_uncert_"+figname+".pdf");
232  // PrintPlots(heff, "efficiency"+figname+".pdf");
233 
234  TCanvas * allUncertCanv = new TCanvas("allUncertCanv", "allUncertCanv", 1440, 1080);
235  allUncertCanv->Divide(2,2);
236  allUncertCanv->cd(1);
237  gPad->SetLeftMargin(0.18);
238  hbg_syst_uncert->Draw("hist ][");
239  allUncertCanv->cd(2);
240  gPad->SetLeftMargin(0.18);
241  hbg_stat_uncert->Draw("hist ][");
242  allUncertCanv->cd(3);
243  gPad->SetLeftMargin(0.18);
244  heff_syst_uncert->Draw("hist ][");
245  allUncertCanv->cd(4);
246  gPad->SetLeftMargin(0.18);
247  hxsec_uncert->Draw("hist ][");
248  allUncertCanv->SaveAs(string(outputPlotDir + "/AllUncertainties_" + figname + ".pdf").c_str());
249  allUncertCanv->SaveAs(string(outputPlotDir + "/AllUncertainties_" + figname + ".png").c_str());
250  delete allUncertCanv;
251 
252  // FindMinimum(hxsec_uncert, minx, maxx);
253  FindMinimum(hxsec_uncert, 1, 100);
254 }//end of plotwith error band
255 
256 //------------------------------------------------------------
257 //------------------------------------------------------------
258 //------------------------------------------------------------
259 // A function to compute the efficiency from a cumulative histogram and a denomenator histogram
260 TH1* ComputeEfficiency(TH1* num_cumu_hist, TH1* denom_hist, std::string name){
261  double denom_val = denom_hist->Integral(0,-1);
262 
263  // Create new TH1 from the orignal, and divide by the integral of the denomenator
264  TH1* eff = (TH1*)num_cumu_hist->Clone(("eff_"+name).c_str());
265  eff->Scale(1/denom_val);
266  eff->GetYaxis()->SetTitle("Efficiency");
267 
268  return eff;
269 }
270 
271 std::vector<TH1*> ComputeEfficiency(std::vector<TH1*> num_cumu_hists,
272  std::vector<TH1*> denom_hists,
273  std::vector<std::string> names)
274 {
275  // Loop through a vector of num/den, and calculate the efficiency of each pair.
276  std::vector<TH1*> effs;
277  for(unsigned int ihist = 0; ihist < num_cumu_hists.size(); ++ihist)
278  effs.push_back(ComputeEfficiency(num_cumu_hists[ihist], denom_hists[ihist], names[ihist]));
279 
280  return effs;
281 }
282 
284  // Copy the numerator (don't actually change the original). Not bothering changing labels for now.
285  ana::Spectrum* scaled_spec = new ana::Spectrum(*num);
286 
287  // Scale the numerator (a cumulative spectrum of selected signal) by the integrated denomenator (true signal)
288  double denom_integral = denom->Integral(denom->POT());
289  scaled_spec->Scale(1/denom_integral);
290 
291  return scaled_spec;
292 }
293 
294 // Alternative function to loop through vectors of cumulative numerators and (non-cumulated) denomenator Spectrums to ComputeEfficiency
295 std::vector<ana::Spectrum*> ComputeEfficiency(const std::vector<ana::Spectrum*> nums, const std::vector<ana::Spectrum*> denoms){
296  // Loop through a set of Spectrums to calculate num/denom
297  std::vector<ana::Spectrum*> eff_specs;
298  for (unsigned int ispec = 0; ispec < nums.size(); ++ispec)
299  eff_specs.push_back(ComputeEfficiency(nums[ispec], denoms[ispec]));
300 
301  return eff_specs;
302 }
303 
304 
305 std::vector<TH1*> ComputeRatio(std::vector<TH1*> num,
306  std::vector<TH1*> denom,
307  // int rebin = 1,
308  std::vector<std::string> names)
309 {
310 
311  std::vector<TH1*> hists_cum ;
312  int nx = num[0]->GetNbinsX()+1;
313  int ny = num[0]->GetNbinsY()+1;
314  int minx = num[0]->GetXaxis()->FindBin(-1);
315  int maxx = num[0]->GetXaxis()->FindBin(1);
316  for(unsigned int ihist = 0; ihist < num.size(); ++ihist){
317 
318  TH1* num_cum = (TH1*)num[ihist]->Clone(("num_" + names[ihist]).c_str());
319  TH1* denum_cum = (TH1*)denom[ihist]->Clone(("denum_" + names[ihist]).c_str());
320  // for(int ix = nx + 1; ix > -1 ; --ix){
321  for(int ix = minx; ix < maxx ; ++ix){
322  const float nnum = num[ihist]->Integral(ix, -1);
323  const float ndenum = denom[ihist]->Integral(ix, -1);
324  // const float n = hists[ihist]->Integral(ix, maxx);
325  //const float n = hists[ihist]->GetBinContent(ix);
326 
327  num_cum->SetBinContent( ix, nnum);
328  denum_cum->SetBinContent( ix, ndenum);
329 
330  //cout << ix << "\t" << nnum << "\t" << ndenum << "\t"<< nnum/ndenum << endl;
331  }
332  num_cum->Divide(denum_cum);
333  num_cum->SetTitle(";MuonID Cut; Fractional Background Contamination");
334  // }// end loop over x bins
335 
336  hists_cum.push_back(num_cum);
337 
338  }// end loop over histograms
339  return hists_cum;
340 }
341 
343  std::vector<TH1*>& ups,
344  std::vector<TH1*>& dns,
345  int col, int errCol,
346  float headroom, bool newaxis,
347  std::string histName = "")
348 {
349  if(col == -1){
350  col = kTotalMCColor;
351  errCol = kTotalMCErrorBandColor;
352  }
353  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
354 
355  nom->SetLineColor(col);
356  nom->GetXaxis()->CenterTitle();
357  nom->GetYaxis()->CenterTitle();
358  if(newaxis) nom->Draw("hist ]["); // Set the axes up
359 
360  double yMax = nom->GetBinContent(nom->GetMaximumBin());
361  //cout << nom->GetMaximumBin() << "\t" << yMax << endl;
362  TGraphAsymmErrors* g = new TGraphAsymmErrors;
363 
364  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
365  const double y = nom->GetBinContent(binIdx);
366  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
367 
368  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
369 
370  double errUp = 0;
371  double errDn = 0;
372  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
373 
374  double hi = (ups[systIdx]->GetBinContent(binIdx)-y);
375  double lo = (y-dns[systIdx]->GetBinContent(binIdx));
376 
377  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
378 
379  errUp += hi*hi;
380  errDn += lo*lo;
381 
382  } // end for systIdx
383 
384  // if ( sqrt(errDn) < 0 || sqrt(errUp) < 0) continue;
385  //std::cout << sqrt(errDn) << "\t"<<sqrt(errUp) << endl;
386  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
387  } // end for i
388 
389  //std::cout << "reached ...................\t"<< endl;
390  g->SetFillColor(errCol);
391  g->Draw("e2 same");
392  g->GetYaxis()->SetRangeUser(0.00001, headroom*yMax);
393  nom->GetYaxis()->SetRangeUser(0.0001, headroom*yMax);
394  gPad->RedrawAxis();
395  nom->Draw("hist ][ same");
396 
397  if (!histName.empty()){
398  g->Write(string(histName + "_errors").c_str());
399  nom->Write(histName.c_str());
400  }
401 
402  }//end of plotwith error band
403 
404 
405 //------------------------------------------------------------
406 //------------------------------------------------------------
407 //------------------------------------------------------------
408 TH1* Cumulative(TH1* hist, std::string name, std::string qualifier){
409  // Create new histogram from the original
410  TH1* hist_cumulative = (TH1*) hist->Clone(("cumu_" + qualifier + name).c_str());
411 
412  // For each bin, calculate the integral from that bin up to the end (ignoring overflow)
413  for (int ix = 0; ix <= hist->GetNbinsX(); ++ix){
414  const float integralAbove = hist->Integral(ix, hist->GetNbinsX());
415  hist_cumulative->SetBinContent(ix, integralAbove);
416  } // End loop over x bins
417 
418  return hist_cumulative;
419 }
420 
421 // Given a hist compute the cumulative hists starting from the
422 // top right corner, ie max x and y bins
423 std::vector<TH1*> Cumulative(std::vector<TH1*> hists,
424  std::vector<std::string> names,
425  std::string qualifier)
426 {
427  std::vector<TH1*> hists_cumu;
428  for(unsigned int ihist = 0; ihist < hists.size(); ++ihist)
429  hists_cumu.push_back(Cumulative(hists[ihist], names[ihist], qualifier));
430 
431  return hists_cumu;
432 }
433 
434 // Function to cumulate a Spectrum. Needs to temporarily convert to TH1 for the logic,
435 // and make a new Spectrum based on that.
437  // Read the input Spectrum
438  double specPOT = spec->POT();
439  double specLivetime = spec->Livetime();
440  TH1 * specHist = spec->ToTH1(specPOT);
441 
442  // Read in the original hist info
443  TH1* hist_cumulative = (TH1*) specHist->Clone();
444  for(int ix = 0; ix <= specHist->GetNbinsX(); ++ix){
445  const float integralAbove = specHist->Integral(ix, specHist->GetNbinsX());
446  hist_cumulative->SetBinContent( ix, integralAbove);
447  }// end loop over x bins
448 
449  ana::Spectrum* cumu_spec = new ana::Spectrum(hist_cumulative, specPOT, specLivetime);
450  return cumu_spec;
451 }
452 
453 // Logic to Cumulative over a vector of Spectrums
454 std::vector<ana::Spectrum*> Cumulative(const std::vector<ana::Spectrum*> specs){
455  std::vector<ana::Spectrum*> cumulative_specs;
456  for (unsigned int i = 0; i < specs.size(); ++i)
457  cumulative_specs.push_back(Cumulative(specs[i]));
458 
459  return cumulative_specs;
460 }
461 
462 //------------------------------------------------------------
463 //------------------------------------------------------------
464 //------------------------------------------------------------
465 
466 void make_muonid_opt( std::string fname = "fOptMuonID.root",
467  std::string outputRoot = "fMuonIDOutput.root",
468  std::string odir = "./")
469 {
470 
471  // Load the input and output files
472  TFile fin(fname.c_str());
473  if (fin.IsZombie() || !fin.IsOpen()) throw std::runtime_error("Could not open input file" + fname);
474  TFile fout(outputRoot.c_str(), "RECREATE");
475  if (fout.IsZombie() || !fout.IsOpen()) throw std::runtime_error("Could not open output file: " + outputRoot);
476  fin.cd();
477 
478  // Set the global outputDir string variable to avoid having to pass it around.
480 
481  std::vector<TString> xtitle = {"; Muon ID Cut ; "};
482 
483  std::vector<std::string> cut_labels = {"all_signal", "sel_signal", "sel_bkgd", "selection"};
484  std::vector<std::string> var_labels = {"muonid"};
485 
486  // gStyle->SetOptStat(111100);
487  // int ivar = 2;
488  for(int ivar = 0; ivar < 1; ivar++){
489  // Switch to input file to load Spectrums
490  fin.cd();
491 
492  // vectors for different cut levels
493  std::vector<ana::Spectrum*> calibup, calibdown;
494  std::vector<ana::Spectrum*> lightup, lightdown;
495  std::vector<ana::Spectrum*> fluxup, fluxdown;
496  std::vector<ana::Spectrum*> ckv, calib_shape;
497  std::vector<ana::Spectrum*> nominal;
498  std::vector<const ana::Spectrum*> xsecup, xsecdown;
499  std::vector<std::vector<ana::Spectrum*> > xsec_all_univs;
500 
501  for(int icut = 0; icut < 4; icut++){
502 
503  nominal.push_back(LoadFromFile<Spectrum>(fname, "nominal/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
504 
505  calibup.push_back(LoadFromFile<Spectrum>(fname, "calibpos/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
506  calibdown.push_back(LoadFromFile<Spectrum>(fname, "calibneg/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
507 
508  lightup.push_back(LoadFromFile<Spectrum>(fname, "lightup/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
509  lightdown.push_back(LoadFromFile<Spectrum>(fname, "lightdown/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
510 
511  ckv.push_back(LoadFromFile<Spectrum>(fname, "ckv/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
512  calib_shape.push_back(LoadFromFile<Spectrum>(fname, "calibshape/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
513 
514  // get the univ spectra
515  std::vector<std::unique_ptr<ana::Spectrum> > xsec_univs;
516 
517  std::cout<<"loading xsec universes for "<<cut_labels[icut]<<"\n";
518 
519  for(unsigned int iuniv = 0; iuniv < N_UNIVERSES; iuniv++)
520  xsec_univs.push_back(Spectrum::LoadFrom((TDirectory*)fin.Get(("xsec/"+var_labels[ivar]+"_"+cut_labels[icut]+"_"+std::to_string(iuniv)).c_str() )));
521 
522  // for xsec errors
523  ana::MultiverseSpectra verses(xsec_univs);
524  xsecup.push_back( verses.UpperSigma(MultiverseSpectra::kBandFromNominal));
525  xsecdown.push_back( verses.LowerSigma(MultiverseSpectra::kBandFromNominal));
526 
527  xsec_all_univs.emplace_back();
528  for(unsigned int iuniv = 0; iuniv < N_UNIVERSES; iuniv++)
529  xsec_all_univs.back().push_back(Spectrum::LoadFrom((TDirectory*)fin.Get(("xsec/"+var_labels[ivar]+"_"+cut_labels[icut]+"_"+std::to_string(iuniv)).c_str() )).release());
530 
531  std::cout<<"Computed xsec error for "<<cut_labels[icut]<<"\n";
532  }// end loop over cuts
533 
534 
535  // Done loading from input file, switch to output for writing.
536  fout.cd();
537 
538 
539  std::vector<std::string> names = {"nominal",
540  "xsecdown",
541  "lightdown", "calibdown", "ckvdown",
542  "calibshapedown",
543  "xsecup",
544  "lightup", "calibup", "ckvup", "calibshapeup"};
545 
546  std::vector<std::string> xsecnames = {"nominal", "xsecdown", "xsecup"};
547  std::vector<std::string> lightnames = {"nominal", "lightdown", "lightup"};
548  std::vector<std::string> calibnames = {"nominal","calibdown", "calibup"};
549  //
550  std::cout<<"Built background cumulative spectra\n";
551 
552  std::vector<TH1*> bg_xsec_cumulative = Cumulative({nominal[2]->ToTH1(pot), xsecdown[2]->ToTH1(pot), xsecup[2]->ToTH1(pot)}, xsecnames, "xsec_bg_");
553 
554  std::vector<TH1*> bg_light_cumulative = Cumulative({nominal[2]->ToTH1(pot), lightdown[2]->ToTH1(pot),lightup[2]->ToTH1(pot)}, lightnames, "light_bg_");
555 
556  std::vector<TH1*> bg_calib_cumulative = Cumulative({nominal[2]->ToTH1(pot), calibdown[2]->ToTH1(pot), calibup[2]->ToTH1(pot)}, calibnames, "calib_bg_");
557 
558  std::vector<TH1*> sel_xsec_cumulative = Cumulative({nominal[3]->ToTH1(pot), xsecdown[3]->ToTH1(pot), xsecup[3]->ToTH1(pot)}, xsecnames, "xsec_sel_");
559 
560  std::vector<TH1*> sel_light_cumulative = Cumulative({nominal[3]->ToTH1(pot), lightdown[3]->ToTH1(pot),lightup[3]->ToTH1(pot)}, lightnames, "light_sel_");
561 
562  std::vector<TH1*> sel_calib_cumulative = Cumulative({nominal[3]->ToTH1(pot), calibdown[3]->ToTH1(pot), calibup[3]->ToTH1(pot)}, calibnames, "calib_sel_");
563 
564  std::vector<TH1*> sig_xsec_cumulative = Cumulative({nominal[1]->ToTH1(pot), xsecdown[1]->ToTH1(pot), xsecup[1]->ToTH1(pot)}, xsecnames, "xsec_sig_");
565 
566  std::vector<TH1*> sig_light_cumulative = Cumulative({nominal[1]->ToTH1(pot), lightdown[1]->ToTH1(pot),lightup[1]->ToTH1(pot)}, lightnames, "light_sig_");
567 
568  std::vector<TH1*> sig_calib_cumulative = Cumulative({nominal[1]->ToTH1(pot), calibdown[1]->ToTH1(pot), calibup[1]->ToTH1(pot)}, calibnames, "calib_sig_");
569 
570  // Do Genie universe logic universe-by-universe
571  std::vector<ana::Spectrum*> sig_genie_univ_cumulative = Cumulative(xsec_all_univs[1]); // Calculate the total number of signal events we see for each possible MuonID cut
572  std::vector<ana::Spectrum*> eff_genie_univ_cumulative = ComputeEfficiency(sig_genie_univ_cumulative, xsec_all_univs[0]); // For each MuonID cut, take the total signal seen / total true signal (without a cut)
573  std::vector<std::unique_ptr<ana::Spectrum> > eff_genie_univ_cumu_unique;
574  for (unsigned int iuniv = 0; iuniv < N_UNIVERSES; ++iuniv)
575  eff_genie_univ_cumu_unique.emplace_back(new ana::Spectrum(*eff_genie_univ_cumulative[iuniv]));
576  ana::MultiverseSpectra eff_univs(eff_genie_univ_cumu_unique);
577 
578  std::vector<TH1*> eff_xsec_cumulative{
579  eff_univs.Nominal()->ToTH1(eff_univs.Nominal()->POT()),
580  eff_univs.LowerSigma(MultiverseSpectra::kBandFromNominal)->ToTH1(eff_univs.Nominal()->POT()),
582  };
583 
584  TCanvas * tempCanv = new TCanvas("tempcanv", "tempcanv", 1600, 900);
585  std::vector<const ana::Spectrum*> all_genie_effs = eff_univs.AllUniverses();
586  double eff_pot = eff_univs.Nominal()->POT();
587  TH1 * nom_hist = all_genie_effs[0]->ToTH1(eff_pot);
588  TH1 * low_hist = eff_xsec_cumulative[1];
589  TH1 * upp_hist = eff_xsec_cumulative[2];
590  nom_hist->SetLineColor(kRed);
591  nom_hist->SetLineWidth(2);
592  low_hist->SetLineColor(kGreen+1);
593  low_hist->SetLineWidth(2);
594  upp_hist->SetLineColor(kMagenta);
595  upp_hist->SetLineWidth(2);
596  TH2D * eff2D = new TH2D("genie_effs_dist", ";MuonID Cut;Efficiency;Num of Universes", 100, -1, 1, 100, 0, 0.25);
597  for (unsigned int iuniv = 1; iuniv < N_UNIVERSES; ++iuniv){
598  TH1 * univ_hist = all_genie_effs[iuniv]->ToTH1(eff_pot);
599  for (int ibin = 1; ibin < univ_hist->GetNbinsX(); ++ibin){
600  eff2D->Fill(univ_hist->GetBinCenter(ibin), univ_hist->GetBinContent(ibin));
601  }
602  univ_hist->SetLineColorAlpha(kBlue, 0.3);
603  univ_hist->SetLineWidth(1);
604  if(iuniv == 1){
605  univ_hist->GetYaxis()->SetRangeUser(0, 0.22);
606  univ_hist->SetTitle(";MuonID Cut;Efficiency");
607  univ_hist->Draw("hist ][");
608  }
609  else
610  univ_hist->Draw("hist ][ SAME");
611  }
612 
613  nom_hist->Draw("hist ][ SAME");
614  low_hist->Draw("hist ][ SAME");
615  upp_hist->Draw("hist ][ SAME");
616  TLegend * tmpLeg = new TLegend();
617  tmpLeg->AddEntry(nom_hist, "Nominal");
618  tmpLeg->AddEntry(low_hist, "Lower #sigma");
619  tmpLeg->AddEntry(upp_hist, "Upper #sigma");
620  tmpLeg->Draw();
621  tempCanv->SaveAs((outputPlotDir + "/genie_effs_all.png").c_str());
622  eff2D->Draw("COLZ");
623  nom_hist->Draw("hist ][ SAME");
624  tempCanv->SaveAs((outputPlotDir + "/genie_effs_dist.png").c_str());
625  delete tempCanv;
626 
627  // Can't do it this way!! Need to compute efficiencies universe by universe
628  // std::vector<TH1*> eff_xsec_cumulative = ComputeEfficiency(sig_xsec_cumulative,
629  // {nominal[0]->ToTH1(pot), xsecdown[0]->ToTH1(pot), xsecup[0]->ToTH1(pot)}, xsecnames);
630 
631  std::vector<TH1*> eff_calib_cumulative = ComputeEfficiency(sig_calib_cumulative,
632  {nominal[0]->ToTH1(pot), calibdown[0]->ToTH1(pot), calibup[0]->ToTH1(pot)}, calibnames);
633 
634  std::vector<TH1*> eff_light_cumulative = ComputeEfficiency(sig_light_cumulative,
635  {nominal[0]->ToTH1(pot), lightdown[0]->ToTH1(pot), lightup[0]->ToTH1(pot)}, lightnames);
636 
637  std::vector<TH1*> bg_nom = {nominal[2]->ToTH1(potmc)};
638 
639  std::vector<TH1*> bg_cumulative =
640  Cumulative({nominal[2]->ToTH1(pot),
641  xsecdown[2]->ToTH1(pot),
642  lightdown[2]->ToTH1(pot),
643  calibdown[2]->ToTH1(pot),
644  ckv[2]->ToTH1(pot),
645  calib_shape[2]->ToTH1(pot),
646  xsecup[2]->ToTH1(pot),
647  lightup[2]->ToTH1(pot),
648  calibup[2]->ToTH1(pot),
649  ckv[2]->ToTH1(pot),
650  calib_shape[2]->ToTH1(pot)},
651  names, "bg_");
652 
653  std::cout<<"Built selected cumulative spectra\n";
654 
655  std::vector<TH1*> sel_cumulative =
656  Cumulative({nominal[3]->ToTH1(pot),
657  xsecdown[3]->ToTH1(pot),
658  lightdown[3]->ToTH1(pot),
659  calibdown[3]->ToTH1(pot),
660  ckv[3]->ToTH1(pot),
661  calib_shape[3]->ToTH1(pot),
662  xsecup[3]->ToTH1(pot),
663  lightup[3]->ToTH1(pot),
664  calibup[3]->ToTH1(pot),
665  ckv[3]->ToTH1(pot),
666  calib_shape[3]->ToTH1(pot)},
667  names, "sel_");
668 
669  std::cout<<"Built signal cumulative spectra\n";
670 
671  std::vector<TH1*> sig_cumulative =
672  Cumulative({nominal[1]->ToTH1(pot),
673  xsecdown[1]->ToTH1(pot),
674  lightdown[1]->ToTH1(pot),
675  calibdown[1]->ToTH1(pot),
676  ckv[1]->ToTH1(pot),
677  calib_shape[1]->ToTH1(pot),
678  xsecup[1]->ToTH1(pot),
679  lightup[1]->ToTH1(pot),
680  calibup[1]->ToTH1(pot),
681  ckv[1]->ToTH1(pot),
682  calib_shape[1]->ToTH1(pot)},
683  names, "sig_");
684 
685  std::cout<<"Computing efficiencies\n";
686 
687  std::vector<TH1*> eff_cumulative =
688  ComputeEfficiency(sig_cumulative,
689  {nominal[0]->ToTH1(pot),
690  xsecdown[0]->ToTH1(pot),
691  lightdown[0]->ToTH1(pot),
692  calibdown[0]->ToTH1(pot),
693  ckv[0]->ToTH1(pot),
694  calib_shape[0]->ToTH1(pot),
695  xsecup[0]->ToTH1(pot),
696  lightup[0]->ToTH1(pot),
697  calibup[0]->ToTH1(pot),
698  ckv[0]->ToTH1(pot),
699  calib_shape[0]->ToTH1(pot)},
700  names);
701  // Override the genie systematic efficiencies to be the ones calculated universe-by-universe.
702  eff_cumulative[1] = eff_xsec_cumulative[1]; // LowerSigma
703  eff_cumulative[6] = eff_xsec_cumulative[2]; // UpperSigma
704  tempCanv = new TCanvas("tempcanv", "tempcanv", 1600, 900);
705  eff_xsec_cumulative[0]->SetLineColor(kRed);
706  eff_xsec_cumulative[1]->SetLineColor(kBlue);
707  eff_xsec_cumulative[2]->SetLineColor(kGreen);
708  eff_xsec_cumulative[2]->SetTitle(";MuonID Cut;Efficiency");
709  eff_xsec_cumulative[2]->GetYaxis()->SetRangeUser(0, 0.22);
710  eff_xsec_cumulative[2]->Draw("hist ]["); // Draw the upper sigma first.
711  eff_xsec_cumulative[0]->Draw("hist ][ SAME");
712  eff_xsec_cumulative[1]->Draw("hist ][ SAME");
713  tmpLeg = new TLegend();
714  tmpLeg->AddEntry(eff_xsec_cumulative[0], "Nominal");
715  tmpLeg->AddEntry(eff_xsec_cumulative[1], "Lower #sigma");
716  tmpLeg->AddEntry(eff_xsec_cumulative[2], "Upper #sigma");
717  tmpLeg->Draw();
718  tempCanv->SaveAs((outputPlotDir + "/genie_effs_sigmas.png").c_str());
719  delete tempCanv;
720 
721  // K, we're gonna print this bin/by/bin for checking
722  cout << "\n\n\nBin Checking:" << endl;
723  for (int ibin = 1; ibin <= eff_xsec_cumulative[0]->GetNbinsX(); ibin++){
724  double nom = eff_xsec_cumulative[0]->GetBinContent(ibin);
725  double low = eff_xsec_cumulative[1]->GetBinContent(ibin);
726  double upp = eff_xsec_cumulative[2]->GetBinContent(ibin);
727 
728  double fraclow = std::abs(low - nom) / nom;
729  double fracupp = std::abs(upp - nom) / nom;
730  cout << "\tbin: " << to_string(ibin) <<
731  "\tnom: " << nom <<
732  "\tlow err: " << std::abs(low - nom) <<
733  "\tfraclow: " << fraclow <<
734  "\tupp err: " << std::abs(upp - nom) <<
735  "\tfracupp: " << fracupp <<
736  "\tmax err: " << std::max(fraclow, fracupp) << endl;
737  }
738  cout << "\n\n\n\n" << endl;
739 
740  std::cout<<"Computing uncertainties\n";
741 
742 
743  std::vector<TH1*> nominal_sel_bg = {nominal[3]->ToTH1(pot), nominal[2]->ToTH1(pot), nominal[1]->ToTH1(pot)};
744 
745 
746  std::vector<TH1*> bg_up =
747  {xsecup[2]->ToTH1(pot),
748  lightup[2]->ToTH1(pot),
749  calibup[2]->ToTH1(pot),
750  ckv[2]->ToTH1(pot),
751  calib_shape[2]->ToTH1(pot)};
752 
753  std::vector<TH1*> bg_down =
754  {xsecdown[2]->ToTH1(pot),
755  lightdown[2]->ToTH1(pot),
756  calibdown[2]->ToTH1(pot),
757  ckv[2]->ToTH1(pot),
758  calib_shape[2]->ToTH1(pot)};
759 
760  std::cout<<"Built selected cumulative spectra\n";
761 
762  std::vector<TH1*> sel_up =
763  {xsecup[3]->ToTH1(pot),
764  lightup[3]->ToTH1(pot),
765  calibup[3]->ToTH1(pot),
766  ckv[3]->ToTH1(pot),
767  calib_shape[3]->ToTH1(pot)};
768 
769  std::vector<TH1*> sel_down =
770  {xsecdown[3]->ToTH1(pot),
771  lightdown[3]->ToTH1(pot),
772  calibdown[3]->ToTH1(pot),
773  ckv[3]->ToTH1(pot),
774  calib_shape[3]->ToTH1(pot)};
775 
776  std::cout<<"Built signal cumulative spectra\n";
777 
778  std::vector<TH1*> truesig_up =
779  {xsecup[1]->ToTH1(pot),
780  lightup[1]->ToTH1(pot),
781  calibup[1]->ToTH1(pot),
782  ckv[1]->ToTH1(pot),
783  calib_shape[1]->ToTH1(pot)};
784 
785  std::vector<TH1*> truesig_down =
786  {xsecdown[1]->ToTH1(pot),
787  lightdown[1]->ToTH1(pot),
788  calibdown[1]->ToTH1(pot),
789  ckv[1]->ToTH1(pot),
790  calib_shape[1]->ToTH1(pot)};
791 
792  std::cout<<"Built true signal cumulative spectra\n";
793 
794  std::vector<TH1*> rup_eff =
795  {eff_cumulative[6],
796  eff_cumulative[7],
797  eff_cumulative[8],
798  eff_cumulative[9],
799  eff_cumulative[10]};
800 
801  std::vector<TH1*> rdn_eff =
802  {eff_cumulative[1],
803  eff_cumulative[2],
804  eff_cumulative[3],
805  eff_cumulative[4],
806  eff_cumulative[5]
807  };
808 
809  std::vector<TH1*> rup_sel =
810  {sel_cumulative[6],
811  sel_cumulative[7],
812  sel_cumulative[8],
813  sel_cumulative[9],
814  sel_cumulative[10]};
815 
816  std::vector<TH1*> rdn_sel =
817  {sel_cumulative[1],
818  sel_cumulative[2],
819  sel_cumulative[3],
820  sel_cumulative[4],
821  sel_cumulative[5]
822  };
823 
824  std::vector<TH1*> rup_bg =
825  {bg_cumulative[6],
826  bg_cumulative[7],
827  bg_cumulative[8],
828  bg_cumulative[9],
829  bg_cumulative[10]};
830 
831  std::vector<TH1*> rdn_bg =
832  {bg_cumulative[1],
833  bg_cumulative[2],
834  bg_cumulative[3],
835  bg_cumulative[4],
836  bg_cumulative[5]
837  };
838 
839 
840  ComputeErrors(bg_cumulative,
841  sel_cumulative,
842  eff_cumulative,
843  bg_nom,
844  "; Muon ID Cut; ",
845  "muonid");
846 
847 
848  ComputeErrors(bg_xsec_cumulative,
849  sel_xsec_cumulative,
850  eff_xsec_cumulative,
851  bg_nom,
852  "; Muon ID Cut; ",
853  "xsec_muonid");
854 
855  ComputeErrors(bg_light_cumulative,
856  sel_light_cumulative,
857  eff_light_cumulative,
858  bg_nom,
859  "; Muon ID Cut; ",
860  "light_muonid");
861 
862  ComputeErrors(bg_calib_cumulative,
863  sel_calib_cumulative,
864  eff_calib_cumulative,
865  bg_nom,
866  "; Muon ID Cut; ",
867  "calib_muonid");
868 
869  std::vector<TH1*> ratio =
870  ComputeRatio({nominal[2]->ToTH1(pot),
871  calibdown[2]->ToTH1(pot),
872  calibup[2]->ToTH1(pot),
873  lightdown[2]->ToTH1(pot),
874  lightup[2]->ToTH1(pot),ckv[2]->ToTH1(pot), calib_shape[2]->ToTH1(pot), xsecup[2]->ToTH1(pot), xsecdown[2]->ToTH1(pot)},
875  {nominal[3]->ToTH1(pot),
876  calibdown[3]->ToTH1(pot),
877  calibup[3]->ToTH1(pot),
878  lightdown[3]->ToTH1(pot),
879  lightup[3]->ToTH1(pot),ckv[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), xsecup[3]->ToTH1(pot), xsecdown[3]->ToTH1(pot)}, names);
880 
881  std::cout<<"Compute ratio\n";
882 
883  std::vector<std::string> names_ratio = {"calib",
884  "light", "ckv",
885  "calibshape",
886  "xsec"};
887 
888 
889  std::vector<TH1*> rup = {ratio[2], ratio[4], ratio[5], ratio[6], ratio[7]};
890 
891  std::vector<TH1*> rdn = {ratio[1], ratio[3], ratio[5], ratio[6], ratio[8]};
892 
893  std::cout<<"Compute ratio.............\n";
894  // int minx = ratio[0]->GetXaxis()->FindBin(-1);
895  //int maxx = ratio[0]->GetXaxis()->FindBin(1);
896 
897  // TH1* hrbkg = (TH1*)ratio[0]->Clone("cum_");
898  //hrbkg->Draw("hist");
899  /* // for(int ix = nx + 1; ix > -1 ; --ix){
900  for(int ix = minx; ix < maxx ; ++ix){
901  const float n = ratio[0]->Integral(ix, -1);
902  hrbkg->SetBinContent( ix, n);
903  }
904  */
905  ratio[0]->GetYaxis()->SetTitleSize(0.048);
906  ratio[0]->GetYaxis()->SetTitleSize(0.048);
907  std::cout<<"Compute ratioxxxxxxxxxxxxxxxxxxxxxxxx\n";
908  new TCanvas();
909  gPad->SetLeftMargin(0.12);
910  // TCanvas* ceff_cum = new TCanvas("ceff_cum", "ceff_cum");
911  PlotWithSystErrorBandWidth(ratio[0], rup, rdn, -1, -1, 1.3f, true, "bkgfractosig");
912  Simulation();
913  gPad->SaveAs(string(outputPlotDir + "/bkgfractosig.pdf").c_str());
914 
915 
916  sel_cumulative[0]->SetTitle(" " + xtitle[ivar] + "Cumulative Selected Events / 8.09 #times 10^{20} POT");
917  bg_cumulative[0]->SetTitle(" " + xtitle[ivar] + "Cumulative Background Events / 8.09 #times 10^{20} POT");
918  eff_cumulative[0]->SetTitle(" " + xtitle[ivar] + "Efficiency");
919  bg_cumulative[0]->GetYaxis()->SetTitleSize(0.048);
920  sel_cumulative[0]->GetYaxis()->SetTitleSize(0.048);
921  nominal_sel_bg[0]-> SetTitle(" ; Muon ID ; Selected Events / 8.09 #times 10^{20} POT");
922  nominal_sel_bg[1]-> SetTitle(" ; Muon ID ; Background Events / 8.09 #times 10^{20} POT");
923  nominal_sel_bg[2]-> SetTitle(" ; Muon ID ; True Sig Events / 8.09 #times 10^{20} POT");
924  nominal_sel_bg[0]-> GetYaxis()->SetTitleSize(0.04);
925  nominal_sel_bg[1]-> GetYaxis()->SetTitleSize(0.04);
926  nominal_sel_bg[2]-> GetYaxis()->SetTitleSize(0.04);
927 
928  TCanvas* ceff = new TCanvas("ceff", "ceff");
929  eff_cumulative[0]->Draw("hist ][");
930  eff_cumulative[0]->SetLineColor(2);
931  ceff->SaveAs(string(outputPlotDir + "/efficiency.pdf").c_str());
932 
933  TCanvas* csel_cum = new TCanvas("csel_cum", "csel_cum");
934  PlotWithSystErrorBandWidth(sel_cumulative[0], rup_sel, rdn_sel, -1, -1, 1.3f, true, "muonid_sel_cumu");
935  csel_cum->SaveAs(string(outputPlotDir + "/muonid_sel_cum.pdf").c_str());
936 
937  TCanvas* cbkg_cum = new TCanvas("cbkg_cum", "cbkg_cum");
938  PlotWithSystErrorBandWidth(bg_cumulative[0], rup_bg, rdn_bg, -1, -1, 1.3f, true, "muonid_bkg_cumu");
939  cbkg_cum->SaveAs(string(outputPlotDir + "/muonid_bkg_cum.pdf").c_str());
940 
941  TCanvas* ceff_cum = new TCanvas("ceff_cum", "ceff_cum");
942  PlotWithSystErrorBandWidth(eff_cumulative[0], rup_eff, rdn_eff, -1, -1, 1.3f, true, "muonid_eff_cumu");
943  ceff_cum->SaveAs(string(outputPlotDir + "/muonid_eff_cum.pdf").c_str());
944 
945  TCanvas* csel = new TCanvas("csel", "csel");
946  PlotWithSystErrorBandWidth(nominal_sel_bg[0], sel_up, sel_down, -1, -1, 1.3f, true, "muonid_sel");
947  csel->SaveAs(string(outputPlotDir + "/muonid_sel.pdf").c_str());
948 
949  TCanvas* cbkg = new TCanvas("cbkg", "cbkg");
950  PlotWithSystErrorBandWidth(nominal_sel_bg[1], bg_up, bg_down, -1, -1, 1.3f, true, "muonid_bkg");
951  cbkg->SaveAs(string(outputPlotDir + "/muonid_bkg.pdf").c_str());
952 
953  TCanvas * ctruesig = new TCanvas("ctruesig", "ctruesig");
954  PlotWithSystErrorBandWidth(nominal_sel_bg[2], truesig_up, truesig_down, -1, -1, 1.3f, true, "muonid_truesig");
955  ctruesig->SaveAs(string(outputPlotDir + "/muonid_truesig.pdf").c_str());
956 
957  // TCanvas* ceff_cum = new TCanvas("ceff_cum", "ceff_cum");
958  //PlotWithSystErrorBandWidth(eff_cumulative[0], rup_eff, rdn_eff, -1, -1, 1.3f, true);
959  //ceff_cum->SaveAs("opt_plots/muonid_eff_cum.pdf");
960 
961  // Save all of the systematic histograms
962  for (unsigned int i = 0; i < cut_labels.size(); i++){
963  nominal[i]->ToTH1(pot)->Write(string("cut_" + cut_labels[i] + "_nominal").c_str());
964  xsecdown[i]->ToTH1(pot)->Write(string("cut_" + cut_labels[i] + "_xsecdown").c_str());
965  lightdown[i]->ToTH1(pot)->Write(string("cut_" + cut_labels[i] + "_lightdown").c_str());
966  calibdown[i]->ToTH1(pot)->Write(string("cut_" + cut_labels[i] + "_calibdown").c_str());
967  ckv[i]->ToTH1(pot)->Write(string("cut_" + cut_labels[i] + "_ckvdown").c_str());
968  calib_shape[i]->ToTH1(pot)->Write(string("cut_" + cut_labels[i] + "_calibshapedown").c_str());
969  xsecup[i]->ToTH1(pot)->Write(string("cut_" + cut_labels[i] + "_xsecup").c_str());
970  lightup[i]->ToTH1(pot)->Write(string("cut_" + cut_labels[i] + "_lightup").c_str());
971  calibup[i]->ToTH1(pot)->Write(string("cut_" + cut_labels[i] + "_calibup").c_str());
972  ckv[i]->ToTH1(pot)->Write(string("cut_" + cut_labels[i] + "_ckvup").c_str());
973  calib_shape[i]->ToTH1(pot)->Write(string("cut_" + cut_labels[i] + "_calibshapeup").c_str());
974  }
975 
976  for(unsigned int i = 0; i < eff_cumulative.size(); ++i){
977  bg_cumulative[i]->Write(string("cumu_bkg_" + names[i]).c_str());
978  eff_cumulative[i]->Write(string("cumu_eff_" + names[i]).c_str());
979  sig_cumulative[i]->Write(string("cumu_sig_" + names[i]).c_str());
980  sel_cumulative[i]->Write(string("cumu_sel_" + names[i]).c_str());
981  }// end loop over effcum
982  } // end of var loop
983 } //end of void
984 
985 
986 
987 
988 
void Simulation()
Definition: tools.h:16
T max(const caf::Proxy< T > &a, T b)
TString fin
Definition: Style.C:24
const XML_Char * name
Definition: expat.h:151
TSpline3 lo("lo", xlo, ylo, 12,"0")
void NDSimulation()
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::string outputPlotDir
correl_xv GetYaxis() -> SetDecimals()
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
constexpr T pow(T x)
Definition: pow.h:75
std::vector< TH1 * > ComputeRatio(std::vector< TH1 * > num, std::vector< TH1 * > denom, std::vector< std::string > names)
TH1 * ComputeEfficiency(TH1 *num_cumu_hist, TH1 *denom_hist, std::string name)
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
TString hists[nhists]
Definition: bdt_com.C:3
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
gargamelle SetTitle("Gargamelle #nu_{e} CC data")
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:237
float abs(float number)
Definition: d0nt_math.hpp:39
std::vector< std::string > cut_labels
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
void PrintPlots(TH1 *&h, TString name)
TSpline3 hi("hi", xhi, yhi, 18,"0")
const double pot
void PlotWithSystErrorBandWidth(TH1 *&nom, std::vector< TH1 * > &ups, std::vector< TH1 * > &dns, int col, int errCol, float headroom, bool newaxis, std::string histName="")
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
Int_t col[ntarg]
Definition: Style.C:29
std::vector< float > Spectrum
Definition: Constants.h:739
TLatex * prelim
Definition: Xsec_final.C:133
double POT() const
Definition: Spectrum.h:227
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
void make_muonid_opt(string inputFilename, string outputDir=".")
void FindMinimum(TH1 *&h, const int minx, const int maxx)
std::vector< TH1D * > nom_hist
int num
Definition: f2_nu.C:119
TH1 * Cumulative(TH1 *hist, std::string name, std::string qualifier)
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal) const
const Spectrum * LowerSigma(BandOptions opt=kBandFromNominal) const
std::vector< std::vector< TH1D * > > univ_hist
const Color_t kTotalMCColor
Definition: Style.h:16
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
std::vector< const Spectrum * > AllUniverses() const
return all universes
void ComputeErrors(std::vector< TH1 * > &bg, std::vector< TH1 * > &sig, std::vector< TH1 * > &eff, std::vector< TH1 * > &bg_nom, TString xylabel, TString figname)
enum BeamMode kGreen
enum BeamMode kBlue
Float_t w
Definition: plot.C:20
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
const double potmc
const unsigned int N_UNIVERSES
const Spectrum * Nominal() const
return the nominal universe
enum BeamMode string