make_vertex_optimiz.C
Go to the documentation of this file.
1 
2 
3 #include <stdio.h>
4 #include <stdlib.h>
5 #include "Utilities/rootlogon.C"
9 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Cuts/TimingCuts.h"
11 #include "CAFAna/Analysis/Plots.h"
12 #include "CAFAna/Analysis/Style.h"
14 #include "CAFAna/Core/Ratio.h"
16 
17 
18 #include "TH1.h"
19 #include "TCanvas.h"
20 
21 #include "math.h"
22 
23 using namespace ana;
24 
25 
26 const double pot = 8.09e20; // POT
27 const double potmc = 3.54e21; // POT
28 
29 const unsigned int N_UNIVERSES = 100;
30 
32 
33 void FindMinimum(TH1*& h, const int minx, const int maxx)
34 {
35  int nx = h->GetXaxis()->GetNbins();
36 
37  double min = 100;
38  int minxbin = 0;
39 
40  for(int ix = minx; ix <= maxx; ix++){
41  if( h->GetBinContent(ix) < min){
42  min = h->GetBinContent(ix);
43  minxbin = ix;
44  }
45  }// end loop over x bins
46 
47  std::cout<<"Best vtx cut "<<h->GetXaxis()->GetBinLowEdge(minxbin)<<"\n";
48 }
49 
50 //------------------------------------------------------------
51 //------------------------------------------------------------
52 //------------------------------------------------------------
53 
54 void PrintPlots(TH1*& h,
55  TString name)
56 {
57  TCanvas *c = new TCanvas();
58  c->SetLeftMargin(0.18);
59  h->GetXaxis()->CenterTitle();
60  h->GetYaxis()->CenterTitle();
61  h->GetYaxis()->SetTitleOffset(1.2);
62  h->SetLineColor(kRed);
63  h->SetLineWidth(2);
64  h->Draw("hist ][");
65  Simulation();
66  //c->SetLogz(1);
67  c->SaveAs(outputPlotDir + "/" + name);
68 }
69 //------------------------------------------------------------
70 //------------------------------------------------------------
71 //------------------------------------------------------------
72 
74 {
75  TLatex* prelim = new TLatex(.9, .95, "NOvA ND Simulation");
76  prelim->SetTextColor(kGray+1);
77  prelim->SetNDC();
78  prelim->SetTextSize(2/30.);
79  prelim->SetTextAlign(32);
80  prelim->Draw();
81 }
82 
83 //------------------------------------------------------------
84 //------------------------------------------------------------
85 //------------------------------------------------------------
86 
87 void ComputeErrors(std::vector<TH1*>& bg,
88  std::vector<TH1*>& sig,
89  std::vector<TH1*>& eff,
90  std::vector<TH1*>& bg_nom,
91  TString xylabel,
92  TString vtx)
93 
94 {
95 
96  TH1* nom = sig[0];
97  TH1* bgnom = bg[0];
98  TH1* effnom = eff[0];
99  TH1* bgstat = bg_nom[0];
100  int nx = nom->GetNbinsX();
101  //new TCanvas;
102  //nom->SetLineWidth(2);
103  //nom->SetLineColor(kRed);
104  //nom->GetXaxis()->CenterTitle();
105  //nom->GetYaxis()->CenterTitle();
106  //nom->GetYaxis()->SetTitle("Events / 8.09 x 10^{20} POT");
107  //nom->Draw("hist");
108  //PrintPlots(nom, vtx+".png");
109 
110  // only explore the space between -5 and +5 LLR
111 
112  //int minx = 0;
113  //int maxx = nom->GetXaxis()->GetNbins();
114 
115  // TString xylabel = "; Vertex Z ; ";
116  //TString vtx = "vtxz";
117  // for (int ivar = 0; ivar < 3, ivar++){
118  //int minx = nom->GetXaxis()->FindBin(-200);
119  //int maxx = nom->GetXaxis()->FindBin(200);
120  // xtitle[ivar].c_str()
121  // std::cout <<" minx "<< minx << " maxx " << maxx << " nbinx " << nx <<std::endl;
122  std::cout << " nbinx " << nx <<std::endl;
123 
124  //Relative Stat. Uncertainty on Sel
125  TH1* hsig_stat_uncert = (TH1*) nom->Clone("hsig_stat_uncert");
126  hsig_stat_uncert->SetTitle(" " + xylabel + "#frac{#sqrt{N_{sel}}}{N_{sel}-N_{bkg}}");
127 
128  //Relative Stat. Uncertainty on Bkg
129  TH1* hbg_stat_uncert = (TH1*) nom->Clone("hbg_stat_uncert");
130  hbg_stat_uncert->SetTitle(" " + xylabel + "#frac{#sqrt{N_{bkg}}}{N_{sel}-N_{bkg}}");
131 
132  //Relative Syst. Uncertainty on Bkg
133  TH1* hbg_syst_uncert = (TH1*) nom->Clone("hbg_syst_uncert");
134  hbg_syst_uncert->SetTitle(" " + xylabel + "#frac{#deltaN^{syst}_{bkg}}{N_{sel}-N_{bkg}}");
135 
136  //Relative Syst. Uncertainty on Eff
137  TH1* heff_syst_uncert = (TH1*) nom->Clone("heff_syst_uncert");
138  heff_syst_uncert->SetTitle(" " + xylabel + "#frac{#delta#varepsilon}{#varepsilon}");
139 
140  TH1* hxsec_uncert = (TH1*) nom->Clone("hxsec_uncert");
141  hxsec_uncert->SetTitle(" " + xylabel + "Relative Cross Section Uncertainty");
142 
143  //Uncertainty on Efficiency
144  TH1* heff_uncert = (TH1*) nom->Clone("heff_uncert");
145  heff_uncert->SetTitle(" " + xylabel + "#delta#varepsilon");
146 
147  //TH1* heff = (TH1*) nom->Clone("heff");
148  // heff->SetTitle("Efficiency" + xylabel + "#varepsilon");
149  hxsec_uncert->SetMaximum(0.3);
150  hsig_stat_uncert->SetMaximum(0.05);
151  hbg_stat_uncert->SetMaximum(0.05);
152  hbg_syst_uncert->SetMaximum(0.1);
153  heff_syst_uncert->SetMaximum(0.15);
154  // }
155  int nsysts = (sig.size() - 1) / 2;
156 
157  // Set everything to 0
158  for(int ix = 1; ix <= nx; ++ix){
159  hbg_syst_uncert->SetBinContent( ix, 0);
160  hbg_stat_uncert->SetBinContent( ix, 0);
161  heff_syst_uncert->SetBinContent( ix, 0);
162  hxsec_uncert->SetBinContent( ix, 0);
163  heff_uncert->SetBinContent( ix, 0);
164  hsig_stat_uncert->SetBinContent( ix, 0);
165  }
166 
167 
168  // for(int ix = minx; ix <= maxx; ++ix){
169 
170  for(int ix = 1; ix <= nx; ++ix){
171 
172  const double nsel = nom->GetBinContent(ix);
173  const double nbg = bgnom->GetBinContent(ix);
174  const double neff = effnom->GetBinContent(ix);
175  double nbgstat = bgstat->GetBinContent(ix);
176  double bg_error_syst = 0;
177  double eff_error_syst = 0;
178 
179 
180  // compute systematic errors in each bin for
181  // background, signal, and efficiency
182 
183  for(int isyst = 1; isyst <= nsysts; ++isyst){
184  // std :: cout << isyst << "......." << nsysts + isyst << std::endl;
185  double ieff_down = eff[isyst]->GetBinContent(ix);
186  double ieff_up = eff[nsysts + isyst]->GetBinContent(ix);
187  double ibg_down = bg[isyst]->GetBinContent(ix);
188  double ibg_up = bg[nsysts + isyst]->GetBinContent(ix);
189 
190  // double bg_error_isyst = std::max( std::fabs(ibg_up - nbg), std::fabs(ibg_down - nbg) );
191  // double eff_error_isyst = std::max( std::fabs(ieff_up - neff), std::fabs(ieff_down - neff) );
192 
193  double bg_error_isyst = std::max(std::abs(ibg_up - nbg), std::abs(ibg_down - nbg));
194  double eff_error_isyst = std::max(std::abs(ieff_up - neff), std::abs(ieff_down - neff));
195 
196  bg_error_syst += std::pow(bg_error_isyst, 2);
197  eff_error_syst += std::pow(eff_error_isyst, 2);
198  } // end for isyst
199 
200  double bg_error_stat = nbgstat;
201  double sel_error_stat = nsel;
202  double totbg_error = sel_error_stat + bg_error_stat + bg_error_syst;
203 
204 
205  // Don't divide by zero.
206  double bg_error_term = 0;
207  if(std::pow( nsel - nbg, 2) != 0)
208  bg_error_term = totbg_error / std::pow( nsel - nbg, 2);
209 
210 
211  double eff_error_term = 0;
212  if(std::pow( neff, 2) != 0)
213  eff_error_term = eff_error_syst / std::pow( neff, 2);
214 
215 
216  double xsec_uncert = eff_error_term + bg_error_term;
217 
218  if (nsel-nbg != 0){
219  hbg_syst_uncert->SetBinContent( ix, std::sqrt(bg_error_syst) / (nsel - nbg));
220  hbg_stat_uncert->SetBinContent( ix, std::sqrt(bg_error_stat) / (nsel - nbg));
221  hsig_stat_uncert->SetBinContent( ix, std::sqrt(sel_error_stat) / (nsel - nbg));
222  }
223 
224  heff_syst_uncert->SetBinContent( ix, std::sqrt(eff_error_term));
225 
226  hxsec_uncert->SetBinContent( ix, std::sqrt(xsec_uncert));
227 
228  heff_uncert->SetBinContent( ix, std::sqrt(eff_error_syst));
229 
230  //heff->SetBinContent( ix, neff);
231  /*
232  std::cout<<" vertex "<<nom->GetXaxis()->GetBinCenter(ix)
233  <<" nsel "<<nsel<<" nbkg "<< nbg << " eff " << neff
234  <<" syst bkg " << std::sqrt(bg_error_syst)
235  << " nsel - nbkg " << nsel - nbg
236  << "bg_syst_uncert " << std::sqrt(bg_error_syst) / (nsel - nbg)
237  <<" xsec error "<<std::sqrt(xsec_uncert)<<std::endl;
238  */
239  } // end loop for x bins
240 
241  // NDSimulation();
242  PrintPlots(hbg_syst_uncert, "bg_syst_uncert_"+vtx+".png");
243  PrintPlots(hbg_stat_uncert, "bg_stat_uncert_"+vtx+".png");
244  PrintPlots(heff_syst_uncert, "eff_syst_uncert_"+vtx+".png");
245  /*
246  TCanvas *c = new TCanvas();
247  c->SetLeftMargin(0.18);
248  h->GetXaxis()->CenterTitle();
249  h->GetYaxis()->CenterTitle();
250  h->GetYaxis()->SetTitleOffset(1.2);
251  h->SetLineColor(kRed);
252  h->SetLineWidth(2);
253  hxsec_uncert->Draw("hist ][");
254  Simulation();
255  //c->SetLogz(1);
256  c->SaveAs( "opt_plots/" + name);
257  */
258  PrintPlots(hxsec_uncert, "xsec_uncert_"+vtx+".png");
259  PrintPlots(heff_uncert, "eff_uncert_"+vtx+".png");
260  PrintPlots(hsig_stat_uncert, "sig_stat_uncert_"+vtx+".png");
261 
262  //FindMinimum(hxsec_uncert, minx, maxx);
263 }//end of plotwith error band
264 
265 
266 //------------------------------------------------------------
267 //------------------------------------------------------------
268 //------------------------------------------------------------
269 // A function to compute the efficiency from a cumulative histogram and a denomenator histogram
270 TH1* ComputeEfficiency(TH1* num_cumu_hist, TH1* denom_hist, std::string name){
271  // Create new TH1 from the orignal, and divide by the integral of the denomenator
272  TH1* eff = (TH1*)num_cumu_hist->Clone(("eff_"+name).c_str());
273  eff->Divide(denom_hist);
274  eff->GetYaxis()->SetTitle("Efficiency");
275  eff->GetZaxis()->SetTitle("Efficiency");
276  return eff;
277 }
278 
279 std::vector<TH1*> ComputeEfficiency(std::vector<TH1*> num_cumu_hists,
280  std::vector<TH1*> denom_hists,
281  std::vector<std::string> names)
282 {
283  // Loop through a vector of num/den, and calculate the efficiency of each pair.
284  std::vector<TH1*> effs;
285  for(unsigned int ihist = 0; ihist < num_cumu_hists.size(); ++ihist)
286  effs.push_back(ComputeEfficiency(num_cumu_hists[ihist], denom_hists[ihist], names[ihist]));
287 
288  return effs;
289 }
290 
292  // Copy the numerator (don't actually change the original). Not bothering changing labels for now.
293  ana::Spectrum* scaled_spec = new ana::Spectrum(*num);
294 
295  // Scale the numerator (a cumulative spectrum of selected signal) by the integrated denomenator (true signal)
296  scaled_spec->Hist()->Divide(denom->Hist());
297 
298  return scaled_spec;
299 }
300 
301 // Alternative function to loop through vectors of cumulative numerators and (non-cumulated) denomenator Spectrums to ComputeEfficiency
302 std::vector<ana::Spectrum*> ComputeEfficiency(const std::vector<ana::Spectrum*> nums, const std::vector<ana::Spectrum*> denoms){
303  // Loop through a set of Spectrums to calculate num/denom
304  std::vector<ana::Spectrum*> eff_specs;
305  for (unsigned int ispec = 0; ispec < nums.size(); ++ispec)
306  eff_specs.push_back(ComputeEfficiency(nums[ispec], denoms[ispec]));
307 
308  return eff_specs;
309 }
310 
311 
313  std::vector<TH1*>& ups,
314  std::vector<TH1*>& dns,
315  int col, int errCol,
316  float headroom, bool newaxis)
317 {
318  if(col == -1){
319  col = kTotalMCColor;
320  errCol = kTotalMCErrorBandColor;
321  }
322  else if(errCol == -1) errCol = col-7; // hopefully a lighter version
323 
324  nom->SetLineColor(col);
325  nom->GetXaxis()->CenterTitle();
326  nom->GetYaxis()->CenterTitle();
327  if(newaxis) nom->Draw("hist ]["); // Set the axes up
328 
329  double yMax = nom->GetBinContent(nom->GetMaximumBin());
330  //cout << nom->GetMaximumBin() << "\t" << yMax << endl;
331  TGraphAsymmErrors* g = new TGraphAsymmErrors;
332 
333  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
334  const double y = nom->GetBinContent(binIdx);
335  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
336 
337  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
338 
339  double errUp = 0;
340  double errDn = 0;
341  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
342 
343  double hi = (ups[systIdx]->GetBinContent(binIdx)-y);
344  double lo = (y-dns[systIdx]->GetBinContent(binIdx));
345 
346  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
347 
348  errUp += hi*hi;
349  errDn += lo*lo;
350 
351  } // end for systIdx
352 
353  // if ( sqrt(errDn) < 0 || sqrt(errUp) < 0) continue;
354  //std::cout << sqrt(errDn) << "\t"<<sqrt(errUp) << endl;
355  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
356  } // end for i
357 
358  //std::cout << "reached ...................\t"<< endl;
359  g->SetFillColor(errCol);
360  g->Draw("e2 same");
361  g->GetYaxis()->SetRangeUser(0.00001, headroom*yMax);
362  nom->GetYaxis()->SetRangeUser(0.0001, headroom*yMax);
363  gPad->RedrawAxis();
364  nom->Draw("hist ][ same");
365 
366  // PrintPlots(nom, vtx+".png");
367 
368 }//end of plotwith error band
369 
370 
371 
372 
373 void make_vertex_optimiz(std::string filename = "/nova/ana/users/bbehera/IncXsec/VtxOpt/fOpt_Vertex.root", std::string odir = ".")
374 {
375 
376  // std::string filename = "/nova/ana/users/bbehera/IncXsec/VtxOpt/fvertex_Opt_ndrespin.root";
377  // "fvertex_Optimize.root";
378  TFile fin(filename.c_str());
379 
381 
382  // std::vector<std::string> xtitle = {"; Vertex X (cm); ", "; Vertex Y (cm); ","; Vertex Z (cm); "} ;
383  std::vector<TString> xtitle = {"; Vertex X (cm); ", "; Vertex Y (cm); ","; Vertex Z (cm); "} ;
384 
385  std::vector<std::string> cut_labels = {"signal", "sel_sig", "sel_bkg", "selection"};
386  std::vector<std::string> var_labels = {"vtxx", "vtxy", "vtxz"};
387  // int ivar = 2;
388  for(int ivar = 0; ivar < 3; ivar++){
389  // vectors for different cut levels
390  std::vector<ana::Spectrum*> calibup, calibdown;
391  std::vector<ana::Spectrum*> lightup, lightdown;
392  std::vector<ana::Spectrum*> fluxup, fluxdown;
393  std::vector<ana::Spectrum*> ckv, calib_shape;
394  std::vector<ana::Spectrum*> nominal;
395  std::vector<const ana::Spectrum*> xsecup, xsecdown;
396  std::vector<std::vector<ana::Spectrum*> > xsec_all_univs;
397 
398  for(int icut = 0; icut < 4; icut++){
399 
400  nominal.push_back(LoadFromFile<Spectrum>(filename, "nominal/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
401 
402  calibup.push_back(LoadFromFile<Spectrum>(filename, "calibpos/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
403  calibdown.push_back(LoadFromFile<Spectrum>(filename, "calibneg/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
404 
405  lightup.push_back(LoadFromFile<Spectrum>(filename, "lightup/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
406  lightdown.push_back(LoadFromFile<Spectrum>(filename, "lightdown/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
407 
408  ckv.push_back(LoadFromFile<Spectrum>(filename, "ckv/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
409  calib_shape.push_back(LoadFromFile<Spectrum>(filename, "calibshape/"+var_labels[ivar]+"_"+cut_labels[icut]).release());
410 
411  // get the univ spectra
412  // std::vector<ana::Spectrum*> flux_univs;
413  std::vector<std::unique_ptr<ana::Spectrum> > xsec_univs;
414 
415  std::cout<<"loading xsec universes for "<<cut_labels[icut]<<"\n";
416 
417  for(unsigned int iuniv = 0; iuniv < N_UNIVERSES; iuniv++){
418  xsec_univs.push_back(Spectrum::LoadFrom((TDirectory*)fin.Get(("xsec/"+var_labels[ivar]+"_"+cut_labels[icut]+"_"+ std::to_string(iuniv)).c_str() )));
419  }// end loop over universes
420 
421  // for flux error
422  //std::vector<TH1*> fupdown;
423  //GetFluxError( nominal[icut], flux_univs, fupdown);
424  //fluxup.push_back(fupdown[0]->Rebin(4));
425  //fluxdown.push_back(fupdown[1]->Rebin(4));
426 
427  // std::cout<<"Computed flux error for "<<cut_labels[icut]<<"\n";
428 
429  // for xsec errors
430  ana::MultiverseSpectra verses(xsec_univs);
431  xsecup.push_back( verses.UpperSigma(MultiverseSpectra::kBandFromNominal));
432  xsecdown.push_back( verses.LowerSigma(MultiverseSpectra::kBandFromNominal));
433 
434  xsec_all_univs.emplace_back();
435  for(unsigned int iuniv = 0; iuniv < N_UNIVERSES; iuniv++)
436  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());
437 
438  std::cout<<"Computed xsec error for "<<cut_labels[icut]<<"\n";
439  }// end loop over cuts
440 
441 
442  std::vector<std::string> names = {"nominal",
443  "xsecdown",// "fluxdown",
444  "lightdown", "calibdown", "ckvdown",
445  "calibshapedown",
446  "xsecup",// "fluxup",
447  "lightup", "calibup", "ckvup", "calibshapeup"};
448 
449  std::vector<std::string> xsecnames = {"nominal", "xsecdown", "xsecup"}; // "fluxdown",
450  std::vector<std::string> lightnames = {"nominal", "lightdown", "lightup"};
451  std::vector<std::string> calibnames = {"nominal","calibdown", "calibup"}; //"ckvdown",
452  std::vector<std::string> ckvnames = {"nominal", "ckvdown", "ckvup"};
453  std::vector<std::string> cshapenames = {"nominal", "calibshapedown", "calibshapeup"};
454  // "calibshapedown",
455 
456  // "lightup", "calibup", "ckvup", "calibshapeup"};
457 
458  std::cout<<"Built background cumulative spectra\n";
459 
460  std::vector<TH1*> bg_xsec = {nominal[2]->ToTH1(pot), xsecdown[2]->ToTH1(pot), xsecup[2]->ToTH1(pot)};
461 
462  std::vector<TH1*> bg_light = {nominal[2]->ToTH1(pot), lightdown[2]->ToTH1(pot),lightup[2]->ToTH1(pot)};
463 
464  std::vector<TH1*> bg_calib = {nominal[2]->ToTH1(pot), calibdown[2]->ToTH1(pot), calibup[2]->ToTH1(pot)};
465 
466  std::vector<TH1*> sel_xsec = {nominal[3]->ToTH1(pot), xsecdown[3]->ToTH1(pot), xsecup[3]->ToTH1(pot)};
467 
468  std::vector<TH1*> sel_light = {nominal[3]->ToTH1(pot), lightdown[3]->ToTH1(pot),lightup[3]->ToTH1(pot)};
469 
470  std::vector<TH1*> sel_calib = {nominal[3]->ToTH1(pot), calibdown[3]->ToTH1(pot), calibup[3]->ToTH1(pot)};
471 
472  std::vector<TH1*> sig_xsec = {nominal[1]->ToTH1(pot), xsecdown[1]->ToTH1(pot), xsecup[1]->ToTH1(pot)};
473 
474  std::vector<TH1*> sig_light = {nominal[1]->ToTH1(pot), lightdown[1]->ToTH1(pot),lightup[1]->ToTH1(pot)};
475 
476  std::vector<TH1*> sig_calib = {nominal[1]->ToTH1(pot), calibdown[1]->ToTH1(pot), calibup[1]->ToTH1(pot)};
477 
478  std::vector<TH1*> eff_xsec = ComputeEfficiency(sig_xsec, {nominal[0]->ToTH1(pot), xsecdown[0]->ToTH1(pot), xsecup[0]->ToTH1(pot)}, xsecnames);
479 
480  std::vector<TH1*> eff_calib = ComputeEfficiency(sig_calib, {nominal[0]->ToTH1(pot), calibdown[0]->ToTH1(pot), calibup[0]->ToTH1(pot)}, calibnames);
481 
482  std::vector<TH1*> eff_light = ComputeEfficiency(sig_light, {nominal[0]->ToTH1(pot), lightdown[0]->ToTH1(pot), lightup[0]->ToTH1(pot)}, lightnames);
483 
484  std::vector<TH1*> bg_nom = {nominal[2]->ToTH1(potmc)};
485 
486  std::vector<TH1*> bg_ckv = {nominal[2]->ToTH1(pot), ckv[2]->ToTH1(pot), ckv[2]->ToTH1(pot)};
487  std::vector<TH1*> bg_cshape = {nominal[2]->ToTH1(pot), calib_shape[2]->ToTH1(pot), calib_shape[2]->ToTH1(pot)};
488  std::vector<TH1*> sel_ckv = {nominal[3]->ToTH1(pot), ckv[3]->ToTH1(pot), ckv[3]->ToTH1(pot)};
489  std::vector<TH1*> sel_cshape = {nominal[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot), calib_shape[3]->ToTH1(pot)};
490  std::vector<TH1*> sig_ckv = {nominal[1]->ToTH1(pot), ckv[1]->ToTH1(pot), ckv[1]->ToTH1(pot)};
491  std::vector<TH1*> sig_cshape = {nominal[1]->ToTH1(pot), calib_shape[1]->ToTH1(pot), calib_shape[1]->ToTH1(pot)};
492  std::vector<TH1*> eff_ckv = ComputeEfficiency(sig_ckv, {nominal[0]->ToTH1(pot), ckv[0]->ToTH1(pot), ckv[0]->ToTH1(pot)}, ckvnames);
493  std::vector<TH1*> eff_cshape = ComputeEfficiency(sig_cshape, {nominal[0]->ToTH1(pot), calib_shape[0]->ToTH1(pot), calib_shape[0]->ToTH1(pot)}, cshapenames);
494 
495  // Do Genie universe logic universe-by-universe
496  // 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
497  std::vector<ana::Spectrum*> eff_genie_univ_cumulative = ComputeEfficiency(xsec_all_univs[1], xsec_all_univs[0]); // For each MuonID cut, take the total signal seen / total true signal (without a cut)
498  std::vector<std::unique_ptr<ana::Spectrum> > eff_genie_univ_cumu_unique;
499  for (unsigned int iuniv = 0; iuniv < N_UNIVERSES; ++iuniv)
500  eff_genie_univ_cumu_unique.emplace_back(new ana::Spectrum(*eff_genie_univ_cumulative[iuniv]));
501  ana::MultiverseSpectra eff_univs(eff_genie_univ_cumu_unique);
502 
503  std::vector<TH1*> eff_xsec_cumulative{
504  eff_univs.Nominal()->ToTH1(eff_univs.Nominal()->POT()),
505  eff_univs.LowerSigma(MultiverseSpectra::kBandFromNominal)->ToTH1(eff_univs.Nominal()->POT()),
507  };
508 
509  std::vector<const ana::Spectrum*> all_eff_univs = eff_univs.AllUniverses();
510 
511  TCanvas * tempCanv = new TCanvas("temp", "temp", 1600, 900);
512  eff_xsec_cumulative[2]->SetLineColor(kRed);
513  eff_xsec_cumulative[2]->GetYaxis()->SetRangeUser(0.0, 0.20);
514  eff_xsec_cumulative[2]->Draw("hist ][");
515  for (unsigned int iuniv = 1; iuniv < N_UNIVERSES; iuniv++){
516  TH1* hist = all_eff_univs[iuniv]->ToTH1(eff_univs.Nominal()->POT());
517  hist->SetLineColor(kBlue);
518  hist->Draw("hist ][ SAME");
519  }
520  eff_xsec_cumulative[0]->SetLineColor(kGreen);
521  eff_xsec_cumulative[1]->SetLineColor(kRed-1);
522  eff_xsec_cumulative[2]->SetLineColor(kRed);
523  eff_xsec_cumulative[2]->Draw("hist ][ SAME");
524  eff_xsec_cumulative[1]->Draw("hist ][ SAME");
525  eff_xsec_cumulative[0]->Draw("hist ][ SAME");
526  tempCanv->SaveAs((outputPlotDir + "/genie_eff_" + var_labels[ivar] +".png").c_str());
527  delete tempCanv;
528 
529  std::vector<TH1*> bg =
530  {nominal[2]->ToTH1(pot),
531  xsecdown[2]->ToTH1(pot),
532  //fluxdown[2]->ToTH1(pot),
533  lightdown[2]->ToTH1(pot),
534  calibdown[2]->ToTH1(pot),
535  ckv[2]->ToTH1(pot),
536  calib_shape[2]->ToTH1(pot),
537  xsecup[2]->ToTH1(pot),
538  //fluxup[2]->ToTH1(pot),
539  lightup[2]->ToTH1(pot),
540  calibup[2]->ToTH1(pot),
541  ckv[2]->ToTH1(pot),
542  calib_shape[2]->ToTH1(pot)};
543  // names, "bg_");
544 
545  std::cout<<"Built selected cumulative spectra\n";
546 
547  std::vector<TH1*> sel =
548  {nominal[3]->ToTH1(pot),
549  xsecdown[3]->ToTH1(pot),
550  //fluxdown[3]->ToTH1(pot),
551  lightdown[3]->ToTH1(pot),
552  calibdown[3]->ToTH1(pot),
553  ckv[3]->ToTH1(pot),
554  calib_shape[3]->ToTH1(pot),
555  xsecup[3]->ToTH1(pot),
556  //fluxup[3]->ToTH1(pot),
557  lightup[3]->ToTH1(pot),
558  calibup[3]->ToTH1(pot),
559  ckv[3]->ToTH1(pot),
560  calib_shape[3]->ToTH1(pot)};
561  // names, "sel_");
562 
563  std::cout<<"Built signal cumulative spectra\n";
564 
565  std::vector<TH1*> sig =
566  {nominal[1]->ToTH1(pot),
567  xsecdown[1]->ToTH1(pot),
568  //fluxdown[1]->ToTH1(pot),
569  lightdown[1]->ToTH1(pot),
570  calibdown[1]->ToTH1(pot),
571  ckv[1]->ToTH1(pot),
572  calib_shape[1]->ToTH1(pot),
573  xsecup[1]->ToTH1(pot),
574  //fluxup[1]->ToTH1(pot),
575  lightup[1]->ToTH1(pot),
576  calibup[1]->ToTH1(pot),
577  ckv[1]->ToTH1(pot),
578  calib_shape[1]->ToTH1(pot)};
579  // names, "sig_");
580 
581  std::cout<<"Computing efficiencies\n";
582 
583  std::vector<TH1*> eff =
584  ComputeEfficiency(sig,
585  {nominal[0]->ToTH1(pot),
586  xsecdown[0]->ToTH1(pot),
587  //fluxdown[0]->ToTH1(pot),
588  lightdown[0]->ToTH1(pot),
589  calibdown[0]->ToTH1(pot),
590  ckv[0]->ToTH1(pot),
591  calib_shape[0]->ToTH1(pot),
592  xsecup[0]->ToTH1(pot),
593  //fluxup[0]->ToTH1(pot),
594  lightup[0]->ToTH1(pot),
595  calibup[0]->ToTH1(pot),
596  ckv[0]->ToTH1(pot),
597  calib_shape[0]->ToTH1(pot)},
598  names);
599 
600  eff[1] = eff_xsec_cumulative[1];
601  eff[6] = eff_xsec_cumulative[2];
602 
603  std::cout<<"Computing uncertainties\n";
604 
605  std::vector<TH1*> rup_eff =
606  {eff[6],
607  eff[7],
608  eff[8],
609  eff[9],
610  eff[10]};
611 
612  std::vector<TH1*> rdn_eff =
613  {eff[1],
614  eff[2],
615  eff[3],
616  eff[4],
617  eff[5]
618  };
619 
620  std::vector<TH1*> rup_sel =
621  {sel[6],
622  sel[7],
623  sel[8],
624  sel[9],
625  sel[10]};
626 
627  std::vector<TH1*> rdn_sel =
628  {sel[1],
629  sel[2],
630  sel[3],
631  sel[4],
632  sel[5]
633  };
634 
635  std::vector<TH1*> rup_bg =
636  {bg[6],
637  bg[7],
638  bg[8],
639  bg[9],
640  bg[10]};
641 
642  std::vector<TH1*> rdn_bg =
643  {bg[1],
644  bg[2],
645  bg[3],
646  bg[4],
647  bg[5]
648  };
649 
650  if (ivar == 0){
651 
652  ComputeErrors(bg,
653  sel,
654  eff,
655  bg_nom,
656  "; Vertex X (cm); ",
657  "vtxx");
658 
659  ComputeErrors(bg_xsec,
660  sel_xsec,
661  eff_xsec,
662  bg_nom,
663  "; Vertex X (cm); ",
664  "xsec_vtxx");
665 
666  ComputeErrors(bg_light,
667  sel_light,
668  eff_light,
669  bg_nom,
670  "; Vertex X (cm); ",
671  "light_vtxx");
672 
673  ComputeErrors(bg_calib,
674  sel_calib,
675  eff_calib,
676  bg_nom,
677  "; Vertex X (cm); ",
678  "calib_vtxx");
679 
680  ComputeErrors(bg_ckv,
681  sel_ckv,
682  eff_ckv,
683  bg_nom,
684  "; Vertex X (cm); ",
685  "ckv_vtxx");
686 
687  ComputeErrors(bg_cshape,
688  sel_cshape,
689  eff_cshape,
690  bg_nom,
691  "; Vertex X (cm); ",
692  "calibshape_vtxx");
693 
694  }else if (ivar == 1){
695 
696  ComputeErrors(bg,
697  sel,
698  eff,
699  bg_nom,
700  "; Vertex Y (cm); ",
701  "vtxy");
702 
703  ComputeErrors(bg_xsec,
704  sel_xsec,
705  eff_xsec,
706  bg_nom,
707  "; Vertex Y (cm); ",
708  "xsec_vtxy");
709 
710  ComputeErrors(bg_light,
711  sel_light,
712  eff_light,
713  bg_nom,
714  "; Vertex Y (cm); ",
715  "light_vtxy");
716 
717  ComputeErrors(bg_calib,
718  sel_calib,
719  eff_calib,
720  bg_nom,
721  "; Vertex Y (cm); ",
722  "calib_vtxy");
723 
724  ComputeErrors(bg_ckv,
725  sel_ckv,
726  eff_ckv,
727  bg_nom,
728  "; Vertex Y (cm); ",
729  "ckv_vtxy");
730 
731  ComputeErrors(bg_cshape,
732  sel_cshape,
733  eff_cshape,
734  bg_nom,
735  "; Vertex Y (cm); ",
736  "calibshape_vtxy");
737 
738  }else{
739 
740  ComputeErrors(bg,
741  sel,
742  eff,
743  bg_nom,
744  "; Vertex Z (cm) ; ",
745  "vtxz");
746 
747  ComputeErrors(bg_xsec,
748  sel_xsec,
749  eff_xsec,
750  bg_nom,
751  "; Vertex Z (cm); ",
752  "xsec_vtxz");
753 
754  ComputeErrors(bg_light,
755  sel_light,
756  eff_light,
757  bg_nom,
758  "; Vertex Z (cm); ",
759  "light_vtxz");
760 
761  ComputeErrors(bg_calib,
762  sel_calib,
763  eff_calib,
764  bg_nom,
765  "; Vertex Z (cm); ",
766  "calib_vtxz");
767 
768  ComputeErrors(bg_ckv,
769  sel_ckv,
770  eff_ckv,
771  bg_nom,
772  "; Vertex Z (cm); ",
773  "ckv_vtxz");
774 
775  ComputeErrors(bg_cshape,
776  sel_cshape,
777  eff_cshape,
778  bg_nom,
779  "; Vertex Z (cm); ",
780  "calibshape_vtxz");
781 
782 
783  } // end of if
784 
785  sel[0]->SetTitle(" " + xtitle[ivar] + "Selected Events / 8.09 #times 10^{20} POT");
786  bg[0]->SetTitle(" " + xtitle[ivar] + "Background Events / 8.09 #times 10^{20} POT");
787  eff[0]->SetTitle(" " + xtitle[ivar] + "Efficiency");
788 
789  TCanvas * c1 = new TCanvas("c1", "c1", 1600, 1200);
790 
791  if (ivar == 0){
792 
793  PlotWithSystErrorBandWidth(sel[0], rup_sel, rdn_sel, -1, -1, 1.3f, true);
794  c1->SaveAs((outputPlotDir + "/sel_vtxx.png").c_str());
795 
796  PlotWithSystErrorBandWidth(bg[0], rup_bg, rdn_bg, -1, -1, 1.3f, true);
797  c1->SaveAs((outputPlotDir + "/bkg_vtxx.png").c_str());
798 
799  PlotWithSystErrorBandWidth(eff[0], rup_eff, rdn_eff, -1, -1, 1.3f, true);
800  c1->SaveAs((outputPlotDir + "/eff_vtxx.png").c_str());
801 
802  }else if (ivar == 1){
803 
804  PlotWithSystErrorBandWidth(sel[0], rup_sel, rdn_sel, -1, -1, 1.3f, true);
805  c1->SaveAs((outputPlotDir + "/sel_vtxy.png").c_str());
806 
807  PlotWithSystErrorBandWidth(bg[0], rup_bg, rdn_bg, -1, -1, 1.3f, true);
808  c1->SaveAs((outputPlotDir + "/bkg_vtxy.png").c_str());
809 
810  PlotWithSystErrorBandWidth(eff[0], rup_eff, rdn_eff, -1, -1, 1.3f, true);
811  c1->SaveAs((outputPlotDir + "/eff_vtxy.png").c_str());
812 
813  }else{
814  PlotWithSystErrorBandWidth(sel[0], rup_sel, rdn_sel, -1, -1, 1.3f, true);
815  c1->SaveAs((outputPlotDir + "/sel_vtxz.png").c_str());
816 
817  PlotWithSystErrorBandWidth(bg[0], rup_bg, rdn_bg, -1, -1, 1.3f, true);
818  c1->SaveAs((outputPlotDir + "/bkg_vtxz.png").c_str());
819 
820  PlotWithSystErrorBandWidth(eff[0], rup_eff, rdn_eff, -1, -1, 1.3f, true);
821  c1->SaveAs((outputPlotDir + "/eff_vtxz.png").c_str());
822  }
823 
824 
825  } // end of var loop
826 
827  /*
828  /// Save some histograms to file for debugging, and
829  /// for gaining intuition about the data.
830  TFile fout("fout_pid_optimization_llr.root", "RECREATE");
831 
832  nominal[3]->ToTH1(pot)->Write();
833  xsecdown[3]->ToTH1(pot)->Write();
834  fluxdown[3]->ToTH1(pot)->Write();
835  lightdown[3]->ToTH1(pot)->Write();
836  calibdown[3]->ToTH1(pot)->Write();
837  ckv[3]->ToTH1(pot)->Write();
838  calib_shape[3]->ToTH1(pot)->Write();
839  xsecup[3]->ToTH1(pot)->Write();
840  fluxup[3]->ToTH1(pot)->Write();
841  lightup[3]->ToTH1(pot)->Write();
842  calibup[3]->ToTH1(pot)->Write();
843  ckv[3]->ToTH1(pot)->Write();
844  calib_shape[3]->ToTH1(pot)->Write();
845 
846 
847  nominal[2]->ToTH1(pot)->Write();
848  xsecdown[2]->ToTH1(pot)->Write();
849  fluxdown[2]->ToTH1(pot)->Write();
850  lightdown[2]->ToTH1(pot)->Write();
851  calibdown[2]->ToTH1(pot)->Write();
852  ckv[2]->ToTH1(pot)->Write();
853  calib_shape[2]->ToTH1(pot)->Write();
854  xsecup[2]->ToTH1(pot)->Write();
855  fluxup[2]->ToTH1(pot)->Write();
856  lightup[2]->ToTH1(pot)->Write();
857  calibup[2]->ToTH1(pot)->Write();
858  ckv[2]->ToTH1(pot)->Write();
859  calib_shape[2]->ToTH1(pot)->Write();
860 
861 
862  // nominal[1]->ToTH1(pot)->Write();
863  // xsecdown[1]->ToTH1(pot)->Write();
864  // fluxdown[1]->ToTH1(pot)->Write();
865  // lightdown[1]->ToTH1(pot)->Write();
866  // calibdown[1]->ToTH1(pot)->Write();
867  // ckv[1]->ToTH1(pot)->Write();
868  // calib_shape[1]->ToTH1(pot)->Write();
869  // xsecup[1]->ToTH1(pot)->Write();
870  // fluxup[1]->ToTH1(pot)->Write();
871  // lightup[1]->ToTH1(pot)->Write();
872  // calibup[1]->ToTH1(pot)->Write();
873  // ckv[1]->ToTH1(pot)->Write();
874  // calib_shape[1]->ToTH1(pot)->Write();
875 
876  // nominal[0]->ToTH1(pot)->Write();
877  // xsecdown[0]->ToTH1(pot)->Write();
878  // fluxdown[0]->ToTH1(pot)->Write();
879  // lightdown[0]->ToTH1(pot)->Write();
880  // calibdown[0]->ToTH1(pot)->Write();
881  // ckv[0]->ToTH1(pot)->Write();
882  // calib_shape[0]->ToTH1(pot)->Write();
883  // xsecup[0]->ToTH1(pot)->Write();
884  // fluxup[0]->ToTH1(pot)->Write();
885  // lightup[0]->ToTH1(pot)->Write();
886  // calibup[0]->ToTH1(pot)->Write();
887  // ckv[0]->ToTH1(pot)->Write();
888  // calib_shape[0]->ToTH1(pot)->Write();
889 
890  for(unsigned int i = 0; i < eff.size(); ++i){
891  bg[i]->Write();
892  eff[i]->Write();
893  sig[i]->Write();
894  sel[i]->Write();
895  }// end loop over effcum
896  */
897 }
898 
899 
900 
901 
902 
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
void PlotWithSystErrorBandWidth(TH1 *&nom, std::vector< TH1 * > &ups, std::vector< TH1 * > &dns, int col, int errCol, float headroom, bool newaxis)
TSpline3 lo("lo", xlo, ylo, 12,"0")
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void NDSimulation()
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
const double pot
string filename
Definition: shutoffs.py:106
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
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)
TSpline3 hi("hi", xhi, yhi, 18,"0")
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
Int_t col[ntarg]
Definition: Style.C:29
const unsigned int N_UNIVERSES
const double potmc
void make_vertex_optimiz(std::string filename="/nova/ana/users/bbehera/IncXsec/VtxOpt/fOpt_Vertex.root", std::string odir=".")
std::vector< float > Spectrum
Definition: Constants.h:610
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 PrintPlots(TH1 *&h, TString name)
void ComputeErrors(std::vector< TH1 * > &bg, std::vector< TH1 * > &sig, std::vector< TH1 * > &eff, std::vector< TH1 * > &bg_nom, TString xylabel, TString vtx)
int num
Definition: f2_nu.C:119
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal) const
c1
Definition: demo5.py:24
const Spectrum * LowerSigma(BandOptions opt=kBandFromNominal) const
const Color_t kTotalMCColor
Definition: Style.h:16
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
std::vector< const Spectrum * > AllUniverses() const
return all universes
enum BeamMode kGreen
enum BeamMode kBlue
Float_t w
Definition: plot.C:20
TH1 * ComputeEfficiency(TH1 *num_cumu_hist, TH1 *denom_hist, std::string name)
std::string outputPlotDir
void FindMinimum(TH1 *&h, const int minx, const int maxx)
const Spectrum * Nominal() const
return the nominal universe
enum BeamMode string