doTemplateFit.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
3 #include "CAFAna/Core/ISyst.h"
8 #include "3FlavorAna/Cuts/NueCutsSecondAna.h"
15 #include "CAFAna/Vars/XsecTunes.h"
16 
17 #include "TMath.h"
18 #include "Math/ProbFunc.h"
19 
26 
27 using namespace ana;
28 
29 const double pot = 8.09e20;
30 
31 const std::string sDir = "/nova/ana/users/mjudah/CrossSection_Feb2019/";
32 const std::string sTemplates = "fAnalysisSpectra.root";
33 
34 const bool plotFitTemplateResults = false;
35 const bool plotFitAnalysisResults = false;
36 const bool plotSignalPrediction = true;
37 
38 
39 struct axis_info
40 {
41  int color;
45 };
46 
48  {
49  {1, 20,"Data", "data"},
50  {632+2, 1,"Total", "mc"},
51  {600+2, 1,"#nu_{e} CC", "nuecc"},
52  {616, 1,"#bar{#nu}_{e} CC", "nuebarcc"},
53  {416-3, 1,"#nu_{#mu} CC", "numucc"},
54  {416+3, 1,"#bar{#nu}_{#mu} CC", "numubarcc"},
55  {800, 1,"NC", "nc"},
56  {400-3, 1,"Other", "other"},
57  {416-3, 1,"Non-Fiducial Signal","nonfidnue"},
58  };
59 
60 //sort Vector
61 bool sort_in_ascending_order (int i,int j) { return (i<j); }
62 
63 
64 float BinSigma(std::vector<double> events, float nsigmas, float pivot)
65  {
66  int pivotbin = 0;
67  double pivotbincenter = 0;
68 
69  std::sort(events.begin(), events.end());
70 
71  for(unsigned int i = 0; i < events.size()-1; i++)
72  if(pivot >= events[i] && pivot < events[i+1]){
73  pivotbin = i;
74  break;
75  }
76  pivotbincenter = pivotbin+0.5;
77 
78  double count_fraction =
80 
81  int nsideevents = 0;
82  int lastbinindex = (int)events.size() - 1;
83 
84  if (nsigmas >= 0) nsideevents = lastbinindex - pivotbin;
85  else nsideevents = pivotbin;
86 
87  int boundIdx = pivotbincenter + count_fraction*(double)nsideevents;
88  int index = 0;
89 
90  if(nsigmas >= 0) index = min(boundIdx, (int)events.size()-1);
91  else index = max(boundIdx, 0);
92 
93  return events.at(index);
94  }
95 
96 
97 
98 //Return 1 sigma bands from Nominal From Multiverse
99 //1 sigma corresponds to 68% of all universes
100 std::vector<TH3F*> CalculateOneSigmaBandFromMultiverse(TH3F* nominal,
101  std::vector<TH3F*> vHist)
102 {
103  TH3F* hResultUpSigma = (TH3F*)nominal->Clone(ana::UniqueName().c_str());
104  TH3F* hResultDownSigma = (TH3F*)nominal->Clone(ana::UniqueName().c_str());
105 
106  for(int xbin= 1; xbin <= hResultUpSigma->GetXaxis()->GetNbins(); xbin++){
107  for(int ybin= 1; ybin <= hResultUpSigma->GetYaxis()->GetNbins(); ybin++){
108  for(int zbin = 1; zbin <= hResultUpSigma->GetZaxis()->GetNbins(); zbin++){
109 
110  double cv_value = nominal->GetBinContent(xbin,ybin,zbin);
111  std::vector<double> multiverse_vals;
112 
113  for(uint iUniv = 0; iUniv < vHist.size(); iUniv++){
114  multiverse_vals.push_back(vHist[iUniv]->GetBinContent
115  (xbin,ybin,zbin));
116  }
117 
118  std::sort(multiverse_vals.begin(),multiverse_vals.end(),
120 
121 
122  float one_sigma_up = BinSigma(multiverse_vals,1,cv_value);
123  float one_sigma_down = BinSigma(multiverse_vals,-1,cv_value);
124 
125  hResultUpSigma->SetBinContent(xbin,ybin,zbin,one_sigma_up);
126  hResultDownSigma->SetBinContent(xbin,ybin,zbin,one_sigma_down);
127  }
128  }
129  }
130 
131  return {(TH3F*)hResultUpSigma->Clone(ana::UniqueName().c_str()),
132  (TH3F*)hResultUpSigma->Clone(ana::UniqueName().c_str())};
133 
134 }
135 
136 std::vector<TH2F*> CalculateOneSigmaBandFromMultiverse(TH2F* nominal,
137  std::vector<TH2F*> vHist)
138 {
139  TH2F* hResultUpSigma = (TH2F*)nominal->Clone(ana::UniqueName().c_str());
140  TH2F* hResultDownSigma = (TH2F*)nominal->Clone(ana::UniqueName().c_str());
141 
142  for(int xbin= 1; xbin <= hResultUpSigma->GetXaxis()->GetNbins(); xbin++){
143  for(int ybin= 1; ybin <= hResultUpSigma->GetYaxis()->GetNbins(); ybin++){
144  double cv_value = nominal->GetBinContent(xbin,ybin);
145  std::vector<double> multiverse_vals;
146 
147  for(uint iUniv = 0; iUniv < vHist.size(); iUniv++){
148  multiverse_vals.push_back(vHist[iUniv]->GetBinContent
149  (xbin,ybin));
150  }
151 
152  std::sort(multiverse_vals.begin(),multiverse_vals.end(),
154 
155 
156  float one_sigma_up = BinSigma(multiverse_vals,1,cv_value);
157  float one_sigma_down = BinSigma(multiverse_vals,-1,cv_value);
158 
159  hResultUpSigma->SetBinContent(xbin,ybin,one_sigma_up);
160  hResultDownSigma->SetBinContent(xbin,ybin,one_sigma_down);
161  }
162  }
163 
164  return {(TH2F*)hResultUpSigma->Clone(ana::UniqueName().c_str()),
165  (TH2F*)hResultUpSigma->Clone(ana::UniqueName().c_str())};
166 }
167 
168 
169 
170 void MakeTemplatePlotsSyst(TH1F* hData,
171  std::vector<TH1F*> hTemplates,
172  std::vector<TH1F*> vUp,std::vector<TH1F*> vDown,
173  std::string plotType, std::string xTitle,
174  std::string dataName, bool isFitted,
175  int xbin, int ybin)
176 {
177  hData->SetMarkerStyle(kChnsInfo[0].marker_style);
178  hData->SetLineColor(kChnsInfo[0].color);
179  hData->SetMarkerColor(kChnsInfo[0].color);
180 
181  for(int i = 0; i < (int)hTemplates.size(); i++){
182  hTemplates[i]->SetLineColor(kChnsInfo[i+1].color);
183  }
184 
185 
186  hTemplates.front()->GetXaxis()->SetTitle(xTitle.c_str());
187  hTemplates.front()->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
188  hTemplates.front()->GetXaxis()->CenterTitle();
189  hTemplates.front()->GetYaxis()->CenterTitle();
190 
191  hData->SetTitle("");
192  hData->GetXaxis()->SetTitle(xTitle.c_str());
193  hData->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
194  hData->GetXaxis()->CenterTitle();
195  hData->GetYaxis()->CenterTitle();
196 
197  TH1F* nominal_clone = (TH1F*)hTemplates[0]->Clone(ana::UniqueName().c_str());
198 
199  TGraphAsymmErrors* g = PlotSystErrorBand(nominal_clone,
200  vUp, vDown,-1,-1,1.3,false);
201 
202  // std::cout << g->GetErrorYhigh(2) << "\t" << g->GetErrorYlow(2) << std::endl;
203  //std::cout << nominal_clone->GetBinContent(2) << std::endl;
204 
205 
206 
207  TH1F* hRatio = (TH1F*)hData->Clone(ana::UniqueName().c_str());
208  hRatio->Divide(hTemplates[0]);
209  hRatio->SetMarkerStyle(20);
210  hRatio->SetLineColor(kBlack);
211  hRatio->GetYaxis()->SetTitle("Data/MC");
212 
213  //Systematic Ratio
214  std::vector<TH1F*> vUpRatio;
215  std::vector<TH1F*> vDownRatio;
216  TH1F* nominal_clone_2 = (TH1F*)hTemplates[0]->Clone(ana::UniqueName().c_str());
217  nominal_clone_2->Divide(nominal_clone);
218  for(uint i = 0; i < vUp.size(); i++){
219  vUpRatio.push_back( (TH1F*)vUp[i]->Clone(ana::UniqueName().c_str()));
220  vDownRatio.push_back( (TH1F*)vDown[i]->Clone(ana::UniqueName().c_str()));
221 
222  vUpRatio[i]->Divide(nominal_clone);
223  vDownRatio[i]->Divide(nominal_clone);
224  }
225 
226  TGraphAsymmErrors* g_ratio = PlotSystErrorBand(nominal_clone_2,
227  vUpRatio, vDownRatio,
228  -1,-1,1.3,false);
229  //Signal Like Templates
230  hTemplates[1]->Add(hTemplates[2]);
231  hTemplates[1]->Add(hTemplates[7]);
232 
233  //NumuCC Like Templates
234  hTemplates[3]->Add(hTemplates[4]);
235 
236 
237  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
238  c1->SetLeftMargin(0.10);
239  c1->SetRightMargin(0.10);
240  TPad *pad1 = new TPad(ana::UniqueName().c_str(), "pad1",0,0,1,1);
241  TPad *pad2 = new TPad(ana::UniqueName().c_str(), "pad2",0,0,1,1);
242  pad1->SetFillStyle(0);
243  pad2->SetFillStyle(0);
244  pad1->SetBottomMargin(0.30);
245  pad2->SetTopMargin(1-0.30);
246  pad2->SetGridx();
247  pad2->SetGridy();
248  c1->cd();
249  pad1->Draw();
250  pad1->cd();
251  hTemplates[0]->GetXaxis()->SetLabelColor(kWhite);
252  hTemplates[0]->GetYaxis()->SetLabelSize(0.03);
253  hTemplates[0]->GetYaxis()->SetTitleSize(0.035);
254  hTemplates[0]->Draw("hist");
255  g->Draw("e2 same");
256  hTemplates[0]->Draw("hist same");
257  hTemplates[1]->Draw("hist same");
258  hTemplates[3]->Draw("hist same");
259  hTemplates[5]->Draw("hist same");
260  hTemplates[6]->Draw("hist same");
261  hData->Draw("same");
262  TLegend* leg = ana::AutoPlaceLegend(0.2,0.2);
263  if(isFitted)leg->SetHeader("Template Weights Applied");
264  leg->AddEntry(hData, "Fake Data", "p");
265  leg->AddEntry(hTemplates[0], "Nominal MC", "l");
266  leg->AddEntry(hTemplates[1], "Signal-Like Template", "l");
267  leg->AddEntry(hTemplates[3], "#nu_{#mu} CC-Like Template", "l");
268  leg->AddEntry(hTemplates[5], "NC Template", "l");
269  leg->AddEntry(hTemplates[6], "Other Template", "l");
270  leg->Draw();
271  SimulationSide();
272  gPad->RedrawAxis();
273  c1->cd();
274  pad2->Draw();
275  pad2->cd();
276  hRatio->GetYaxis()->SetLabelSize(0.03);
277  hRatio->GetYaxis()->SetTitleSize(0.035);
278  hRatio->GetYaxis()->SetRangeUser(0.4,1.5);
279  hRatio->Draw("same");
280  g_ratio->Draw("e2 same");
281  hRatio->Draw("same");
282  pad2->Draw("same");
283  gPad->RedrawAxis();
284  c1->cd();
285  c1->Update();
286  std::string fitted_string = "nominal_prediction";
287  if(isFitted)fitted_string = "template_fit_prediction";
288  char outname[50];
289  sprintf(outname, "%s%s_%s_%s_%i_%i.png","Plots/", dataName.c_str(),
290  plotType.c_str(), fitted_string.c_str(), xbin,ybin);
291  c1->SaveAs(outname);
292  c1->Close();
293 
294 }
295 
296 void PlotFrom2D(TH2F* hData2D,TH2F* hNominal2D,
297  std::vector<TH2F*> hFitResults2D,
298  std::vector<TH2F*> hSystsUp2D,
299  std::vector<TH2F*> hSystsDown2D,
300  std::string plotType, std::string xTitle,
302 {
303  for(int xbin = 1; xbin <= hData2D->GetXaxis()->GetNbins(); xbin++){
304  TH1F* hData =
305  (TH1F*)hData2D->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
306 
307  if(hData->Integral() == 0) continue;
308 
309  TH1F* hNominal =
310  (TH1F*)hNominal2D->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
311  std::vector<TH1F*> hFitResults1D;
312  for(uint i = 0; i < hFitResults2D.size(); i++){
313  hFitResults1D.push_back((TH1F*)hFitResults2D[i]->ProjectionY(ana::UniqueName().c_str(),xbin,xbin));
314  }
315  std::vector<TH1F*> hSystsUp1D;
316  std::vector<TH1F*> hSystsDown1D;
317  for(uint i = 0; i < hSystsUp2D.size(); i++){
318  hSystsUp1D.push_back((TH1F*)hSystsUp2D[i]->ProjectionY(ana::UniqueName().c_str(),xbin,xbin));
319  hSystsDown1D.push_back((TH1F*)hSystsDown2D[i]->ProjectionY(ana::UniqueName().c_str(),xbin,xbin));
320  }
321 
322  hData->SetMarkerStyle(kChnsInfo[0].marker_style);
323  hData->SetLineColor(kChnsInfo[0].color);
324  hData->SetMarkerColor(kChnsInfo[0].color);
325 
326  hFitResults1D.front()->SetLineColor(kBlue);
327  hFitResults1D.front()->SetMarkerColor(kBlue);
328 
329  hNominal->SetLineColor(kRed);
330  hNominal->SetMarkerColor(kRed);
331 
332  hFitResults1D.front()->GetXaxis()->SetTitle(xTitle.c_str());
333  hFitResults1D.front()->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
334  hFitResults1D.front()->GetXaxis()->CenterTitle();
335  hFitResults1D.front()->GetYaxis()->CenterTitle();
336 
337  hData->SetTitle("");
338  hData->GetXaxis()->SetTitle(xTitle.c_str());
339  hData->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
340  hData->GetXaxis()->CenterTitle();
341  hData->GetYaxis()->CenterTitle();
342 
343  TH1F* nominal_clone = (TH1F*)hNominal->Clone(ana::UniqueName().c_str());
344 
345  TGraphAsymmErrors* g = PlotSystErrorBand(nominal_clone,
346  hSystsUp1D, hSystsDown1D,
347  -1,-1,1.3,false);
348 
349  TH1F* hRatio = (TH1F*)hData->Clone(ana::UniqueName().c_str());
350  hRatio->Divide(hNominal);
351  hRatio->SetMarkerStyle(20);
352  hRatio->SetLineColor(kRed);
353  hRatio->GetYaxis()->SetTitle("Data/MC");
354 
355  TH1F* hRatio2 = (TH1F*)hData->Clone(ana::UniqueName().c_str());
356  hRatio->Divide(hFitResults1D[0]);
357  hRatio->SetMarkerStyle(20);
358  hRatio->SetLineColor(kBlue);
359  hRatio->GetYaxis()->SetTitle("Data/MC");
360 
361 
362  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
363  c1->SetLeftMargin(0.10);
364  c1->SetRightMargin(0.10);
365  TPad *pad1 = new TPad(ana::UniqueName().c_str(), "pad1",0,0,1,1);
366  TPad *pad2 = new TPad(ana::UniqueName().c_str(), "pad2",0,0,1,1);
367  pad1->SetFillStyle(0);
368  pad2->SetFillStyle(0);
369  pad1->SetBottomMargin(0.30);
370  pad2->SetTopMargin(1-0.30);
371  pad2->SetGridx();
372  pad2->SetGridy();
373  c1->cd();
374  pad1->Draw();
375  pad1->cd();
376  hData->GetXaxis()->SetLabelColor(kWhite);
377  hData->GetYaxis()->SetLabelSize(0.03);
378  hData->GetYaxis()->SetTitleSize(0.035);
379  hData->GetXaxis()->SetRangeUser(1.0,10.0);
380  hData->Draw();
381  hNominal->Draw("hist same");
382  g->Draw("e2 same");
383  hNominal->Draw("hist same");
384  hFitResults1D.front()->Draw("same");
385  hData->Draw("same");
386  hFitResults1D.front()->Draw("same");
387  TLegend* leg = ana::AutoPlaceLegend(0.2,0.2);
388  leg->AddEntry(hData, "Fake Data", "p");
389  leg->AddEntry(hNominal, "Nominal MC", "l");
390  leg->AddEntry(hFitResults1D[0], "Weighted MC", "pl");
391  leg->Draw();
392  SimulationSide();
393  gPad->RedrawAxis();
394  c1->cd();
395  pad2->Draw();
396  pad2->cd();
397  hRatio->GetXaxis()->SetRangeUser(1.0,10.0);
398  hRatio->GetYaxis()->SetLabelSize(0.03);
399  hRatio->GetYaxis()->SetTitleSize(0.035);
400  hRatio->GetYaxis()->SetRangeUser(0.4,1.5);
401  hRatio->Draw("same");
402  hRatio2->Draw("same");
403  pad2->Draw("same");
404  gPad->RedrawAxis();
405  c1->cd();
406  c1->Update();
407  char outname[50];
408  sprintf(outname, "%s%s_%s_%i.png","Plots/", dataName.c_str(),
409  plotType.c_str(), xbin);
410  c1->SaveAs(outname);
411  c1->Close();
412  }//loop over xbins
413 }
414 
415 void PlotWithRatio2D(TH2F* hResultA, TH2F* hTruthA, TH2F* hNominalA,
417  float xmin, float xmax,
419  std::string plotType, std::string plotName,
420  std::string varName,
421  std::string legend_results = "Fake Data",
422  std::string legend_truth = "Fake Data Truth",
423  std::string legend_nominal = "Nominal MC")
424 {
425 
426 
427  for(int xbin = 1; xbin <= hResultA->GetXaxis()->GetNbins(); xbin++){
428 
429  TH1F* hResult =
430  (TH1F*)hResultA->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
431  TH1F* hTruth =
432  (TH1F*)hTruthA->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
433  TH1F* hNominal =
434  (TH1F*)hNominalA->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
435 
436  if(hTruth->Integral() == 0) continue;
437 
438  std::stringstream low_X;
439  low_X << std::fixed << std::setprecision(2) <<
440  hResultA->GetXaxis()->GetBinLowEdge(xbin);
441  std::stringstream hi_X;
442  hi_X << std::fixed << std::setprecision(2) <<
443  hResultA->GetXaxis()->GetBinUpEdge(xbin);
444 
445  std::string caption = low_X.str() + " #leq cos #theta_{e} < "+
446  hi_X.str();
447 
448 
449  hResult->GetYaxis()->SetTitle(ytitle.c_str());
450  hResult->GetXaxis()->SetTitle(xtitle.c_str());
451  hResult->SetTitle(caption.c_str());
452  hResult->GetYaxis()->CenterTitle();
453  hResult->GetXaxis()->CenterTitle();
454  hResult->GetYaxis()->SetTitleOffset(1.15);
455  hResult->SetMarkerStyle(20);
456  hResult->SetMarkerColor(kBlack);
457  hTruth->SetLineColor(kRed);
458  hTruth->SetMarkerColor(kRed);
459 
460  hNominal->SetLineColor(kBlue);
461  hNominal->SetMarkerColor(kBlue);
462 
463  TH1F* hResultClone = (TH1F*)hResult->Clone(ana::UniqueName().c_str());
464  TH1F* hNominalClone = (TH1F*)hNominal->Clone(ana::UniqueName().c_str());
465  hResultClone->SetTitle("");
466  hResultClone->Divide(hTruth);
467  hNominalClone->Divide(hTruth);
468 
469  hResult->GetXaxis()->SetRangeUser(xmin,xmax);
470  hResultClone->GetXaxis()->SetRangeUser(xmin,xmax);
471 
472 
473  //Find Y limits
474  double ymax = hResult->GetMaximum();
475 
476  if(ymax < hTruth->GetMaximum())
477  hResult->GetYaxis()->SetRangeUser(0,hTruth->GetMaximum()*1.3);
478  else
479  hResult->GetYaxis()->SetRangeUser(0,ymax*1.3);
480 
481 
482  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str(),
483  ana::UniqueName().c_str());
484  c1->SetLeftMargin(0.10);
485  c1->SetRightMargin(0.10);
486  c1->SetBottomMargin(0.10);
487  TPad *pad1 = new TPad(ana::UniqueName().c_str(), "pad1",0,0,1,1);
488  TPad *pad2 = new TPad(ana::UniqueName().c_str(), "pad2",0,0,1,1);
489  pad1->SetFillStyle(0);
490  pad2->SetFillStyle(0);
491  pad1->SetBottomMargin(0.30);
492  pad2->SetTopMargin(1-0.30);
493  pad2->SetGridx();
494  pad2->SetGridy();
495  c1->cd();
496  pad1->Draw();
497  pad1->cd();
498  hResult->GetXaxis()->SetLabelColor(kWhite);
499  hResult->GetYaxis()->SetLabelSize(0.03);
500  hResult->GetYaxis()->SetTitleSize(0.035);
501  hResult->Draw();
502  // PrintChiSq(caption_chi.c_str());
503  hResult->Draw("same");
504  hTruth->Draw("hist same");
505  hNominal->Draw("same");
506  TLegend* leg = ana::AutoPlaceLegend(0.23,0.23);
507  leg->AddEntry(hResult, legend_results.c_str(), "p");
508  leg->AddEntry(hTruth, legend_truth.c_str(), "l");
509  leg->AddEntry(hNominal, legend_nominal.c_str(), "p");
510  leg->Draw();
511  SimulationSide();
512  gPad->RedrawAxis();
513  c1->cd();
514  pad2->Draw();
515  pad2->cd();
516  hResultClone->GetYaxis()->SetTitle("Ratio to Truth");
517  hResultClone->GetYaxis()->SetLabelSize(0.03);
518  hResultClone->GetYaxis()->SetTitleSize(0.035);
519  hResultClone->GetYaxis()->SetRangeUser(0.4,1.5);
520  hResultClone->Draw("same");
521  hNominalClone->Draw("same");
522  pad2->Draw("same");
523  gPad->RedrawAxis();
524  c1->cd();
525  c1->Update();
526  char outname[50];
527  sprintf(outname, "%s%s_%s_%s_%s_%i.png","Plots/", dataName.c_str(),
528  plotType.c_str(), plotName.c_str(), varName.c_str(),xbin);
529  c1->SaveAs(outname);
530  c1->Close();
531 
532  delete c1;
533  delete leg;
534  delete hResult;
535  delete hTruth;
536  delete hResultClone;
537 
538  }//loop over xbins
539 
540 }
541 
542 
543 void MakeCovarianceCorrelationPlots(TH2F* hCov, TH2F* hCorr,
544  int xbin, int ybin)
545 {
546  hCov->GetXaxis()->CenterTitle();
547  hCov->GetYaxis()->CenterTitle();
548  hCov->GetZaxis()->CenterTitle();
549 
550  hCorr->GetXaxis()->CenterTitle();
551  hCorr->GetYaxis()->CenterTitle();
552  hCorr->GetZaxis()->CenterTitle();
553 
554  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
555  c1->SetLeftMargin(0.15);
556  c1->SetRightMargin(0.15);
557  c1->SetBottomMargin(0.15);
558  gStyle->SetTitleAlign(33);
559  hCov->Draw("COLZ");
560  Simulation();
561  char outname[50];
562  sprintf(outname, "%s%s_%s_%i_%i.png","Plots/",
563  "covariance", "template_bins", xbin,ybin);
564  c1->SaveAs(outname);
565  c1->Close();
566  delete c1;
567 
568  TCanvas *c2 = new TCanvas(ana::UniqueName().c_str());
569  c2->SetLeftMargin(0.15);
570  c2->SetRightMargin(0.15);
571  c2->SetBottomMargin(0.15);
572  gStyle->SetTitleAlign(33);
573  hCorr->Draw("COLZ");
574  Simulation();
575  sprintf(outname, "%s%s_%s_%i_%i.png","Plots/",
576  "correlation", "template_bins", xbin,ybin);
577  c2->SaveAs(outname);
578  c2->Close();
579  delete c2;
580 
581 }
582 
583 
584 
586  bool shouldSave = true)
587 {
588  TFile* fTemplates = new TFile((sDir+sTemplates).c_str(),"read");
589 
590  //////////////////////////////////////////////////////////////////////////
591  //////// Grab Systematic Templates For Covariance Matrix /////////////////
592  //////////////////////////////////////////////////////////////////////////
593 
594  NueCCIncCrossSectionAnalysis nominal_ana =
595  *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("Nominal/Analysis"));
596  NueCCIncCrossSectionAnalysis light_dwn_ana =
597  *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("LightDwn/Analysis"));
598  NueCCIncCrossSectionAnalysis light_up_ana =
599  *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("LightUp/Analysis"));
600  NueCCIncCrossSectionAnalysis calib_dwn_ana =
601  *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("CalibDwn/Analysis"));
602  NueCCIncCrossSectionAnalysis calib_up_ana =
603  *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("CalibUp/Analysis"));
604  NueCCIncCrossSectionAnalysis calib_shape_ana =
605  *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("CalibShape/Analysis"));
606  NueCCIncCrossSectionAnalysis ckv_ana =
607  *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("Ckv/Analysis"));
608 
609  std::vector<std::unique_ptr<NueCCIncCrossSectionAnalysis>> genie_ana;
610  std::vector<std::unique_ptr<NueCCIncCrossSectionAnalysis>> ppfx_ana;
611 
612  for(int i = 0; i < 100; i++){
613  genie_ana.push_back(NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory(Form("GENIE/%s_%i","Universe",i))));
614  ppfx_ana.push_back(NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory(Form("PPFX/%s_%i","Universe",i))));
615  }
616 
617  TH3F* nominal_template = nominal_ana.getTemplate3D(costhetabins,
618  eelecbins,
619  pidbins);
620 
621  TH3F* light_dwn_template = light_dwn_ana.getTemplate3D(costhetabins,
622  eelecbins,
623  pidbins);
624 
625  TH3F* light_up_template = light_up_ana.getTemplate3D(costhetabins,
626  eelecbins,
627  pidbins);
628 
629  TH3F* calib_dwn_template = calib_dwn_ana.getTemplate3D(costhetabins,
630  eelecbins,
631  pidbins);
632 
633  TH3F* calib_up_template = calib_up_ana.getTemplate3D(costhetabins,
634  eelecbins,
635  pidbins);
636 
637  TH3F* calib_shape_template = calib_shape_ana.getTemplate3D(costhetabins,
638  eelecbins,
639  pidbins);
640 
641  TH3F* ckv_template = ckv_ana.getTemplate3D(costhetabins,
642  eelecbins,
643  pidbins);
644 
645  std::vector<TH3F*> genie_templates;
646  std::vector<TH3F*> ppfx_templates;
647 
648  for(int i = 0; i < 100; i++){
649  genie_templates.push_back(genie_ana[i]->getTemplate3D(costhetabins,
650  eelecbins,
651  pidbins));
652  ppfx_templates.push_back(ppfx_ana[i]->getTemplate3D(costhetabins,
653  eelecbins,
654  pidbins));
655  }
656 
657  std::vector<TH3F*> vMultiverse;
658  for(int i = 0; i < 100; i++){
659  vMultiverse.push_back((TH3F*)genie_templates[i]->Clone
660  (ana::UniqueName().c_str()));
661  }
662  for(int i = 0; i < 100; i++){
663  vMultiverse.push_back((TH3F*)ppfx_templates[i]->Clone
664  (ana::UniqueName().c_str()));
665  }
666 
667 
668  //////////////////////////////////////////////////////////////////////////
669  /////////////////// Grab Templates From Nominal///////////////////////////
670  //////////////////////////////////////////////////////////////////////////
671  NueCCIncCrossSectionTemplates fit_templates =
672  *NueCCIncCrossSectionTemplates::LoadFrom(fTemplates->GetDirectory("Nominal/Templates"));
673 
674  std::vector<TH3F*>nominal_fit_templates =
675  fit_templates.getAllTemplates3D(costhetabins,
676  eelecbins,
677  pidbins);
678 
679  std::vector<TH3F*>nominal_fit_analysis =
680  fit_templates.getAllAnalysis3D(costhetabins,
681  eelecbins,
682  enubins);
683 
684  std::vector<TH3F*>nominal_analysis_comparison;
685  for(uint i = 0; i < nominal_fit_analysis.size(); i++){
686  nominal_analysis_comparison.push_back((TH3F*)
687  nominal_fit_analysis[i]->Clone
688  (ana::UniqueName().c_str()));
689  }
690 
691 
692  NueCCIncTemplateFitter
693  fitter((TH3F*)nominal_template->Clone(ana::UniqueName().c_str()),
694  nominal_fit_templates,
695  nominal_fit_analysis,
696  (TH3F*)nominal_template->Clone(ana::UniqueName().c_str()),
697  {(TH3F*)light_up_template->Clone(ana::UniqueName().c_str()),
698  (TH3F*)calib_up_template->Clone(ana::UniqueName().c_str())},
699  {(TH3F*)light_dwn_template->Clone(ana::UniqueName().c_str()),
700  (TH3F*)calib_dwn_template->Clone(ana::UniqueName().c_str())},
701  {(TH3F*)ckv_template->Clone(ana::UniqueName().c_str()),
702  (TH3F*)calib_shape_template->Clone(ana::UniqueName().c_str())},
703  vMultiverse,
704  21);
705 
706  std::vector<TH3F*> hFitResults = fitter.doFullFit();
707 
708 
709  //////////////////////////////////////////////////////////////////////////
710  /////////////////// Plot Template Results ////////////////////////////////
711  //////////////////////////////////////////////////////////////////////////
712 
714  std::vector<TH3F*> hGenieSysts3D =
715  CalculateOneSigmaBandFromMultiverse(nominal_template,
716  genie_templates);
717 
718  std::vector<TH3F*> hPPFXSysts3D =
719  CalculateOneSigmaBandFromMultiverse(nominal_template,
720  ppfx_templates);
721 
722  //Plot Template Fit Results
723  for(int xbin = 1; xbin <= hFitResults[0]->GetXaxis()->GetNbins(); xbin++){
724  for(int ybin = 1; ybin <= hFitResults[0]->GetYaxis()->GetNbins(); ybin++){
725 
726  std::vector<TH1F*> hPID_Nominal = fitter.getPIDTemplates(xbin,ybin);
727  std::vector<TH1F*> hPID_Fitted = fitter.getPIDFitTemplates(xbin,ybin);
728  TH2F* hCovariance = fitter.getCovarianceMatrixTemplateBins(xbin,ybin);
729  TH2F* hCorrelation = fitter.getCorrelationMatrixTemplateBins(xbin,ybin);
730 
731  if(hPID_Nominal[0]->Integral() == 0) continue;
732 
733  std::vector<TH1F*> hSystUpHolder1D = {
734  (TH1F*)light_up_template->ProjectionZ(ana::UniqueName().c_str(),
735  xbin,xbin,ybin,ybin),
736  (TH1F*)calib_up_template->ProjectionZ(ana::UniqueName().c_str(),
737  xbin,xbin,ybin,ybin),
738  (TH1F*)ckv_template->ProjectionZ(ana::UniqueName().c_str(),
739  xbin,xbin,ybin,ybin),
740  (TH1F*)calib_shape_template->ProjectionZ(ana::UniqueName().c_str(),
741  xbin,xbin,ybin,ybin),
742  (TH1F*)hGenieSysts3D[0]->ProjectionZ(ana::UniqueName().c_str(),
743  xbin,xbin,ybin,ybin),
744  (TH1F*)hPPFXSysts3D[0]->ProjectionZ(ana::UniqueName().c_str(),
745  xbin,xbin,ybin,ybin)
746  };
747 
748  std::vector<TH1F*> hSystDownHolder1D = {
749  (TH1F*)light_dwn_template->ProjectionZ(ana::UniqueName().c_str(),
750  xbin,xbin,ybin,ybin),
751  (TH1F*)calib_dwn_template->ProjectionZ(ana::UniqueName().c_str(),
752  xbin,xbin,ybin,ybin),
753  (TH1F*)nominal_template->ProjectionZ(ana::UniqueName().c_str(),
754  xbin,xbin,ybin,ybin),
755  (TH1F*)nominal_template->ProjectionZ(ana::UniqueName().c_str(),
756  xbin,xbin,ybin,ybin),
757  (TH1F*)hGenieSysts3D[1]->ProjectionZ(ana::UniqueName().c_str(),
758  xbin,xbin,ybin,ybin),
759  (TH1F*)hPPFXSysts3D[1]->ProjectionZ(ana::UniqueName().c_str(),
760  xbin,xbin,ybin,ybin)
761  };
762 
763 
764  MakeTemplatePlotsSyst((TH1F*)nominal_template->ProjectionZ
765  (ana::UniqueName().c_str(),xbin,xbin,ybin,ybin),
766  hPID_Nominal,
767  hSystUpHolder1D,
768  hSystDownHolder1D,
769  "templates", "CVNe 2017", "Nominal",
770  false,xbin,ybin);
771 
772  MakeTemplatePlotsSyst((TH1F*)nominal_template->ProjectionZ
773  (ana::UniqueName().c_str(),xbin,xbin,ybin,ybin),
774  hPID_Fitted,
775  hSystUpHolder1D,
776  hSystDownHolder1D,
777  "templates", "CVNe 2017", "Nominal",
778  true,xbin,ybin);
779 
780  MakeCovarianceCorrelationPlots(hCovariance, hCorrelation, xbin,ybin);
781 
782  /*
783  delete hCovariance;
784  delete hCorrelation;
785 
786  for (int i =0; i< hPID_Nominal.size();i++)
787  {
788  delete (hPID_Nominal[i]);
789  delete (hPID_Fitted[i]);
790  }
791  hPID_Nominal.clear();
792  hPID_Fitted.clear();
793  */
794  }
795  }
796  } //if plot fit results
797 
798  //////////////////////////////////////////////////////////////////////////
799  /////////////////// Plot Fit Results in Electron Kinematic Space//////////
800  //////////////////////////////////////////////////////////////////////////
801  //Project all analysis plots down to "yx" axis
802  std::vector<TH2F*> hFitResults2D;
803  for(uint i = 0; i < hFitResults.size();i++){
804  TH2F* hHolder = (TH2F*)hFitResults[i]->Project3D("yx");
805  hHolder->SetName(ana::UniqueName().c_str());
806  hFitResults2D.push_back((TH2F*)hHolder->Clone(ana::UniqueName().c_str()));
807  }
808 
809  TH2F* data_analysis_2D = nominal_ana.getAnalysis2D(costhetabins,
810  eelecbins,
811  enubins);
812 
813 
814  TH2F* nominal_analysis_2D = nominal_ana.getAnalysis2D(costhetabins,
815  eelecbins,
816  enubins);
818  TH2F* light_dwn_analysis_2D = light_dwn_ana.getAnalysis2D(costhetabins,
819  eelecbins,
820  enubins);
821 
822  TH2F* light_up_analysis_2D = light_up_ana.getAnalysis2D(costhetabins,
823  eelecbins,
824  enubins);
825 
826  TH2F* calib_dwn_analysis_2D = calib_dwn_ana.getAnalysis2D(costhetabins,
827  eelecbins,
828  enubins);
829 
830  TH2F* calib_up_analysis_2D = calib_up_ana.getAnalysis2D(costhetabins,
831  eelecbins,
832  enubins);
833 
834  TH2F* calib_shape_analysis_2D = calib_shape_ana.getAnalysis2D(costhetabins,
835  eelecbins,
836  enubins);
837 
838  TH2F* ckv_analysis_2D = ckv_ana.getAnalysis2D(costhetabins,
839  eelecbins,
840  enubins);
841 
842  std::vector<TH2F*> genie_analysis_2D;
843  std::vector<TH2F*> ppfx_analysis_2D;
844 
845  for(int i = 0; i < 100; i++){
846  genie_analysis_2D.push_back(genie_ana[i]->getAnalysis2D(costhetabins,
847  eelecbins,
848  enubins));
849  ppfx_analysis_2D.push_back(ppfx_ana[i]->getAnalysis2D(costhetabins,
850  eelecbins,
851  enubins));
852  }
853  //if(plotFitAnalysisResults){
854 
855  std::vector<TH2F*> hGenieSysts2D =
856  CalculateOneSigmaBandFromMultiverse(nominal_analysis_2D,
857  genie_analysis_2D);
858 
859  std::vector<TH2F*> hPPFXSysts2D =
860  CalculateOneSigmaBandFromMultiverse(nominal_analysis_2D,
861  ppfx_analysis_2D);
862 
863 
864  std::vector<TH2F*> hSystUpHolder2D = {
865  (TH2F*)light_up_analysis_2D->Clone(ana::UniqueName().c_str()),
866  (TH2F*)calib_up_analysis_2D->Clone(ana::UniqueName().c_str()),
867  (TH2F*)ckv_analysis_2D->Clone(ana::UniqueName().c_str()),
868  (TH2F*)calib_shape_analysis_2D->Clone(ana::UniqueName().c_str()),
869  (TH2F*)hGenieSysts2D[0]->Clone(ana::UniqueName().c_str()),
870  (TH2F*)hPPFXSysts2D[0]->Clone(ana::UniqueName().c_str()),
871  };
872 
873  std::vector<TH2F*> hSystDownHolder2D = {
874  (TH2F*)light_dwn_analysis_2D->Clone(ana::UniqueName().c_str()),
875  (TH2F*)calib_dwn_analysis_2D->Clone(ana::UniqueName().c_str()),
876  (TH2F*)nominal_analysis_2D->Clone(ana::UniqueName().c_str()),
877  (TH2F*)nominal_analysis_2D->Clone(ana::UniqueName().c_str()),
878  (TH2F*)hGenieSysts2D[1]->Clone(ana::UniqueName().c_str()),
879  (TH2F*)hPPFXSysts2D[1]->Clone(ana::UniqueName().c_str()),
880  };
881 
882 
883  PlotFrom2D(data_analysis_2D,
884  nominal_analysis_2D,
885  hFitResults2D,
886  hSystUpHolder2D,
887  hSystDownHolder2D,
888  "analysis",
889  "Reconstructed Electron Energy, E_{e} (GeV)",
890  "Nominal");
891 
892  }//if plot analysis fit results
893 
894  //////////////////////////////////////////////////////////////////////////
895  /////////////////// Unfold Signal Estimate to True Space /////////////////
896  //////////////////////////////////////////////////////////////////////////
897  TH2F* hSignalEst = (TH2F*)hFitResults2D[1]->Clone("hSignalEst");
898 
899  TH2F* hUnfoldedSignalEst = nominal_ana.doUnfolding2D(hSignalEst, 6);
900 
901 
902  TH2F* hTrueRecoSignal =
903  (TH2F*)nominal_ana.getSelectedSignalPrediction2D(costhetabins,
904  eelecbins);
905 
906  TH2F* hTrueRecoSignal_nominal =
907  (TH2F*)nominal_ana.getSelectedSignalPrediction2D(costhetabins,
908  eelecbins);
909 
910 
911 
913  PlotWithRatio2D(hUnfoldedSignalEst, hTrueRecoSignal,hTrueRecoSignal_nominal,
914  "Signal Events/8.09 #times 10^{20} POT",
915  "True Electron Energy, E_{e} (GeV)",
916  1.0, 10.0,
917  dataName,
918  "reco_signal", "comparison",
919  "double", "Weighted MC - Signal Prediction",
920  "True Signal from Fake Data",
921  "Nominal MC - Signal Prediction");
922  }
923 
924  //////////////////////////////////////////////////////////////////////////
925  /////////////////// Efficiency //////////////////////////////////////////
926  //////////////////////////////////////////////////////////////////////////
927  TH2F* nominal_efficiency_2D =
928  nominal_ana.getEfficiency2D(costhetabins,eelecbins);
929  TH2F* light_dwn_efficiency_2D =
930  light_dwn_ana.getEfficiency2D(costhetabins,eelecbins);
931 
932  TH2F* light_up_efficiency_2D =
933  light_up_ana.getEfficiency2D(costhetabins,eelecbins);
934 
935  TH2F* calib_dwn_efficiency_2D =
936  calib_dwn_ana.getEfficiency2D(costhetabins,eelecbins);
937 
938  TH2F* calib_up_efficiency_2D =
939  calib_up_ana.getEfficiency2D(costhetabins,eelecbins);
940 
941  TH2F* calib_shape_efficiency_2D =
942  calib_shape_ana.getEfficiency2D(costhetabins,eelecbins);
943 
944  TH2F* ckv_efficiency_2D = ckv_ana.getEfficiency2D(costhetabins,eelecbins);
945 
946  std::vector<TH2F*> genie_efficiency_2D;
947  std::vector<TH2F*> ppfx_efficiency_2D;
948 
949  for(int i = 0; i < 100; i++){
950  genie_efficiency_2D.push_back(genie_ana[i]->getEfficiency2D(costhetabins,
951  eelecbins));
952  ppfx_efficiency_2D.push_back(ppfx_ana[i]->getEfficiency2D(costhetabins,
953  eelecbins));
954  }
955 
956  //////////////////////////////////////////////////////////////////////////
957  /////////////////// Efficiency - MRE Correction /////////////////////////
958  //////////////////////////////////////////////////////////////////////////
959 
960 
961 
962 
963  std::cout << "DONE" << std::endl;
964  fTemplates->Close();
965 }
void doTemplateFit(std::string dataName="nominal", bool shouldSave=true)
void Simulation()
Definition: tools.h:16
const bool plotFitAnalysisResults
Definition: doTemplateFit.C:35
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
void PlotFrom2D(TH2F *hData2D, TH2F *hNominal2D, std::vector< TH2F * > hFitResults2D, std::vector< TH2F * > hSystsUp2D, std::vector< TH2F * > hSystsDown2D, std::string plotType, std::string xTitle, std::string dataName)
const ana::Binning eelecbins
double Integral(const Spectrum &s, const double pot, int cvnregion)
bool sort_in_ascending_order(int i, int j)
Definition: doTemplateFit.C:61
const ana::Binning costhetabins
void PlotWithRatio2D(TH2F *hResultA, TH2F *hTruthA, TH2F *hNominalA, std::string ytitle, std::string xtitle, float xmin, float xmax, std::string dataName, std::string plotType, std::string plotName, std::string varName, std::string legend_results="Fake Data", std::string legend_truth="Fake Data Truth", std::string legend_nominal="Nominal MC")
std::vector< TH3F * > CalculateOneSigmaBandFromMultiverse(TH3F *nominal, std::vector< TH3F * > vHist)
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
const ana::Binning pidbins
Double_t ymax
Definition: plot.C:25
c2
Definition: demo5.py:33
const Binning enubins
const std::string sTemplates
Definition: doTemplateFit.C:32
std::string title
Definition: doTemplateFit.C:43
const double j
Definition: BetheBloch.cxx:29
void MakeTemplatePlotsSyst(TH1F *hData, std::vector< TH1F * > hTemplates, std::vector< TH1F * > vUp, std::vector< TH1F * > vDown, std::string plotType, std::string xTitle, std::string dataName, bool isFitted, int xbin, int ybin)
void MakeCovarianceCorrelationPlots(TH2F *hCov, TH2F *hCorr, int xbin, int ybin)
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
const bool plotFitTemplateResults
Definition: doTemplateFit.C:34
double BinSigma(std::vector< double > events, double nsigma, double pivot)
void SimulationSide()
Definition: rootlogon.C:137
const double pot
Definition: doTemplateFit.C:29
Eigen::ArrayXd ProjectionY(const Eigen::MatrixXd &mat)
Helper for WeightingVariable.
const int kcNumChns
Definition: NueCCIncCuts.h:293
std::string cutName
Definition: doTemplateFit.C:44
const std::string sDir
Definition: doTemplateFit.C:31
c1
Definition: demo5.py:24
const bool plotSignalPrediction
Definition: doTemplateFit.C:36
std::string dataName
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
TPad * pad2
Definition: analysis.C:13
const axis_info kChnsInfo[kcNumChns+1]
Definition: doTemplateFit.C:47
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
enum BeamMode kBlue
return_type< T_y, T_loc, T_scale >::type normal_cdf(const T_y &y, const T_loc &mu, const T_scale &sigma)
Definition: normal_cdf.hpp:39
void events(int which)
Definition: Cana.C:52
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
TGraphAsymmErrors * PlotSystErrorBand(TH1F *&nom, std::vector< TH1F * > &ups, std::vector< TH1F * > &dns, int col, int errCol, float headroom, bool newaxis)
int marker_style
Definition: doTemplateFit.C:42
TPad * pad1
Definition: analysis.C:13
enum BeamMode string
unsigned int uint