makeXSecPlots1D.C
Go to the documentation of this file.
5 #include "NDAna/nuebarcc_inc/NuebarCCIncCrossSectionFunctions.h"
6 
7 #include "RooUnfoldResponse.h"
8 #include "RooUnfoldBayes.h"
9 
10 #include "Utilities/rootlogon.C"
13 #include "TSpline.h"
14 using namespace ana;
15 
16 //const double pot = 8.09e20;
17 const double pot = kAna2018FHCPOT;
18 const double nnucleons = 3.02225e+31;
19 
20 const std::string sDir = "/nova/ana/users/ddoyle/NuebarCrossSection_020319/";
21 const std::string sSignalEst = "CrossSection_TemplateFitResults_fakedata_development_stat.root";
22 const std::string sAnalysis = "SpectraNominal/fCrossSection_Nominal_FakeData_prongcvn_new.root";
23 const std::string sSysts = "SpectraSystematics/fCrossSection_EfficiencyFlux_Systs.root";
24 const std::string sGenie = "/cvmfs/nova.opensciencegrid.org/externals/genie_xsec/v2_12_10/NULL/DefaultPlusMECWithNC/data/xsec_graphs.root";
25 void MakePlots1D(TH1F* hData, TH1F* hFit, TH1F* hMC,
27  std::string dataName, std::string xsecResult,
28  std::string varName)
29 {
30  TH1F* hD = (TH1F*)hData->Clone(ana::UniqueName().c_str());
31  TH1F* hF = (TH1F*)hFit->Clone(ana::UniqueName().c_str());
32  TH1F* hN = (TH1F*)hMC->Clone(ana::UniqueName().c_str());
33 
34  hD->SetMarkerStyle(20);
35  hF->SetMarkerStyle(22);
36  hF->SetMarkerColor(kBlue);
37  hF->SetLineColor(kBlue);
38  hN->SetLineColor(kRed);
39  hF->SetTitle("");
40  hN->SetTitle("");
41  hD->GetYaxis()->SetTitle(ytitle.c_str());
42  hD->GetXaxis()->SetTitle(xtitle.c_str());
43 
44 
45  double chi_fit = 0;
46  int ndf_fit = 0, igood_fit = 0;
47  hD->Chi2TestX(hF,chi_fit, ndf_fit, igood_fit,"WW");
48  double chi_nom = 0;
49  int ndf_nom = 0, igood_nom = 0;
50  hD->Chi2TestX(hN,chi_nom, ndf_nom, igood_nom,"WW");
51 
52  std::string caption = "Fitted #chi^{2}: " + std::to_string(chi_fit) + "\n"
53  + "Nominal #chi^{2}: " + std::to_string(chi_nom);
54 
55 
56  TH1F* hratio_fit = (TH1F*)hF->Clone(ana::UniqueName().c_str());
57  TH1F* hratio_nom = (TH1F*)hN->Clone(ana::UniqueName().c_str());
58 
59  hD->GetXaxis()->SetRangeUser(1.0,10.0);
60  hF->GetXaxis()->SetRangeUser(1.0,10.0);
61  hN->GetXaxis()->SetRangeUser(1.0,10.0);
62  hratio_fit->GetXaxis()->SetRangeUser(1.0,10.0);
63  hratio_nom->GetXaxis()->SetRangeUser(1.0,10.0);
64 
65  hratio_fit->Divide(hD);
66  hratio_nom->Divide(hD);
67 
68  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str(), "c1");
69  c1->SetLeftMargin(0.10);
70  c1->SetRightMargin(0.10);
71  TPad *pad1 = new TPad(ana::UniqueName().c_str(), "pad1",0,0,1,1);
72  TPad *pad2 = new TPad(ana::UniqueName().c_str(), "pad2",0,0,1,1);
73  pad1->SetFillStyle(0);
74  pad2->SetFillStyle(0);
75  pad1->SetBottomMargin(0.30);
76  pad2->SetTopMargin(1-0.30);
77  pad2->SetGridx();
78  pad2->SetGridy();
79  c1->cd();
80  pad1->Draw();
81  pad1->cd();
82  hD->SetTitle("");
83  hD->GetXaxis()->SetLabelColor(kWhite);
84  hD->GetYaxis()->SetLabelSize(0.03);
85  hD->GetYaxis()->SetTitleSize(0.035);
86  hD->SetTitle("");
87  hD->Draw();
88  PrintChiSq(caption.c_str());
89  hD->Draw("hist same");
90  hF->Draw("same");
91  hN->Draw("same");
92  TLegend* leg = ana::AutoPlaceLegend(0.3,0.3);
93  leg->AddEntry(hD, "Fake Data", "p");
94  leg->AddEntry(hN, "Nominal MC", "l");
95  leg->AddEntry(hF, "Fitted MC", "p");
96  leg->Draw();
97  Simulation();
98  gPad->RedrawAxis();
99  c1->cd();
100  pad2->Draw();
101  pad2->cd();
102  hratio_fit->GetYaxis()->SetTitle("MC/Data");
103  hratio_fit->GetXaxis()->SetTitle(xtitle.c_str());
104  hratio_fit->GetXaxis()->CenterTitle();
105  hratio_fit->GetYaxis()->SetLabelSize(0.03);
106  hratio_fit->GetYaxis()->SetTitleSize(0.035);
107  hratio_fit->GetYaxis()->SetRangeUser(0.9,1.1);
108  hratio_fit->Draw("same");
109  hratio_nom->Draw("hist same");
110  pad2->Draw("same");
111  gPad->RedrawAxis();
112  c1->cd();
113  c1->Update();
114  char outname[50];
115  sprintf(outname, "%s%s_%s.png", "AnalysisPlots/",dataName.c_str(),
116  varName.c_str());
117  c1->SaveAs(outname);
118  c1->Close();
119 }
120 
121 void MakePlots(TH1F* hFitted1D, std::vector<TH1F*> hSysts, TH1F* hTruth1D,
123  std::string dataName, std::string xsecResult,
124  std::string varName, bool calcCov = false)
125 {
126 
127  TH1F* hFitted = (TH1F*)hFitted1D->Clone(ana::UniqueName().c_str());
128  TH1F* hTruth = (TH1F*)hTruth1D->Clone(ana::UniqueName().c_str());
129 
130  std::cout << hSysts.size() << std::endl;
131 
132  hFitted->SetLineColor(kBlack);
133  hFitted->SetMarkerColor(kBlack);
134  hFitted->SetMarkerStyle(20);
135  hFitted->GetXaxis()->SetRangeUser(1.0,10.0);
136  hFitted->SetTitle("");
137  hTruth->SetLineColor(kRed);
138  hTruth->SetMarkerColor(kRed);
139 
140 
141  hFitted->GetYaxis()->SetTitle(ytitle.c_str());
142  hFitted->GetXaxis()->SetTitle(xtitle.c_str());
143 
144  std::vector<TH1F*> hRelUncert;
145  for(int i = 0; i <= (int)hSysts.size();i++){
146  TH1F* hHolder = (TH1F*)hFitted->Clone(ana::UniqueName().c_str());
147  for(int ibin = 0; ibin <= hFitted->GetXaxis()->GetNbins(); ibin++){
148 
149  if(i < (int)hSysts.size()){
150  double uncert = hSysts[i]->GetBinError(ibin);
151  double total = hSysts[i]->GetBinContent(ibin);
152 
153  if(total == 0)
154  hHolder->SetBinContent(ibin,0);
155  else
156  hHolder->SetBinContent(ibin,uncert/total);
157  }
158  if( i == (int) hSysts.size()){
159  double uncert = hFitted->GetBinError(ibin);
160  double total = hFitted->GetBinContent(ibin);
161 
162 
163  if(total == 0)
164  hHolder->SetBinContent(ibin,0);
165  else
166  hHolder->SetBinContent(ibin,uncert/total);
167  }
168  }
169  hRelUncert.push_back(hHolder);
170  }
171 
172  hRelUncert[0]->SetLineColor(kGreen+2); //calibshape
173  hRelUncert[1]->SetLineColor(kOrange-3); //ckv
174  hRelUncert[2]->SetLineColor(kRed - 2); //genie
175  hRelUncert[3]->SetLineColor(kCyan-3); //ppfx
176  hRelUncert[4]->SetLineColor(kBlue-4); //light
177  hRelUncert[5]->SetLineColor(kMagenta+1); //calib
178  hRelUncert[6]->SetLineColor(kBlack); //tot
179 
180  hSysts[0]->SetLineColor(kGreen+2); //calibshape
181  hSysts[1]->SetLineColor(kOrange-3); //ckv
182  hSysts[2]->SetLineColor(kRed - 2); //genie
183  hSysts[3]->SetLineColor(kCyan-3); //ppfx
184  hSysts[4]->SetLineColor(kBlue-4); //light
185  hSysts[5]->SetLineColor(kMagenta+1); //calib
186  hSysts[0]->SetMarkerColor(kGreen+2); //calibshape
187  hSysts[1]->SetMarkerColor(kOrange-3); //ckv
188  hSysts[2]->SetMarkerColor(kRed - 2); //genie
189  hSysts[3]->SetMarkerColor(kCyan-3); //ppfx
190  hSysts[4]->SetMarkerColor(kBlue-4); //light
191  hSysts[5]->SetMarkerColor(kMagenta+1); //calib
192 
193  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str(),"c1");
194  hFitted->Draw();
195  hTruth->Draw("hist same");
196  //hSysts[0]->Draw("same");
197  //hSysts[1]->Draw("same");
198  //hSysts[2]->Draw("same");
199  //hSysts[3]->Draw("same");
200  //hSysts[4]->Draw("same");
201  //hSysts[5]->Draw("same");
202  TLegend* leg = ana::AutoPlaceLegend(0.3,0.3);
203  leg->AddEntry(hFitted, "Fitted Signal Estimate", "pl");
204  leg->AddEntry(hTruth,"Truth", "l");
205  leg->Draw();
206  char outname[50];
207  sprintf(outname, "%s%s_%s_%s.png","AnalysisPlots/", dataName.c_str(),
208  varName.c_str(), xsecResult.c_str());
209  c1->SaveAs(outname);
210  c1->Close();
211  //delete c1;
212  //delete leg;
213 
214  TCanvas* c2 = new TCanvas(ana::UniqueName().c_str(),"c2");
215  hRelUncert[6]->GetYaxis()->SetTitle("Relative Uncertainty");
216  hRelUncert[6]->GetYaxis()->SetRangeUser(0,0.3);
217  hRelUncert[6]->Draw("hist");
218  hRelUncert[5]->Draw("hist same");
219  hRelUncert[4]->Draw("hist same");
220  hRelUncert[3]->Draw("hist same");
221  hRelUncert[2]->Draw("hist same");
222  hRelUncert[1]->Draw("hist same");
223  hRelUncert[0]->Draw("hist same");
224  hRelUncert[6]->Draw("hist same");
225  TLegend* leg_systs = ana::AutoPlaceLegend(0.3,0.3);
226  leg_systs->AddEntry(hRelUncert[6], "Total", "l");
227  leg_systs->AddEntry(hRelUncert[5], "Calibration Normalization", "l");
228  leg_systs->AddEntry(hRelUncert[4], "Light Level", "l");
229  leg_systs->AddEntry(hRelUncert[3], "Flux", "l");
230  leg_systs->AddEntry(hRelUncert[2], "#nu-Interaction", "l");
231  leg_systs->AddEntry(hRelUncert[1], "Cherenkov Variation", "l");
232  leg_systs->AddEntry(hRelUncert[0], "Calibration Shape", "l");
233  leg_systs->Draw();
234  sprintf(outname, "%s%s_%s_%s_%s.png","AnalysisPlots/", dataName.c_str(),
235  varName.c_str(), "relative_uncertantiy",xsecResult.c_str());
236  c2->SaveAs(outname);
237  //delete c2;
238  //delete leg_systs;
239 }
240 
242  std::string xsecResult = "energy")
243 {
244 
245  TFile *fSignalEst = new TFile((sDir+sSignalEst).c_str(),"read");
246  TFile* fAnalysis = new TFile((sDir+sAnalysis).c_str(),"read");
247  TFile* fSysts = new TFile((sDir+sSysts).c_str(),"read");
248  TFile* fGenie = new TFile((sGenie).c_str(),"read");
249  ////////////////////////////////////////////////////////////////////
250  ////////////////////// Compare Fit Results /////////////////////////
251  ////////////////////////////////////////////////////////////////////
252  char name[50];
253  sprintf(name, "mc_%s_%s_%s_postfit",dataName.c_str(),
254  "tot", chns[0].name.c_str()); // mc
255  Spectrum sFitTotal = *Spectrum::LoadFrom(fSignalEst, name);
256  TH3F* hFitTotal3D = (TH3F*)ana::ToTH3(sFitTotal,pot,kPOT,
259  TH1F* hFitTotal1D = (TH1F*)hFitTotal3D->ProjectionZ("hFitTotal1D");
260 
261 
262  //Get Total Fake Data
263  sprintf(name, "mc_%s_%s_%s", dataName.c_str(), "analysis", "mc");
264  Spectrum sDataTotal = *Spectrum::LoadFrom(fAnalysis, name);
265  TH3F* hDataTotal3D = (TH3F*)ana::ToTH3(sDataTotal,pot,kPOT,
268  TH1F* hDataTotal1D = (TH1F*)hDataTotal3D->ProjectionZ("hDataTotal1D");
269 
270  //Get Total Nominal MC
271  sprintf(name, "mc_%s_%s_%s", "nominal", "analysis", "mc");
272  Spectrum sMCTotal = *Spectrum::LoadFrom(fAnalysis, name);
273  TH3F* hMCTotal3D = (TH3F*)ana::ToTH3(sMCTotal,pot,kPOT,
276  TH1F* hMCTotal1D = (TH1F*)hMCTotal3D->ProjectionZ("hMCTotal1D");
277 
278  MakePlots1D(hDataTotal1D, hFitTotal1D, hMCTotal1D,
279  "Events/8.09 #times 10^{20} POT",
280  "Reconstructed Neutrino Energy, E_{#nu} (GeV)",
281  dataName, xsecResult,"total_fit_comparison");
282 
283 
284  ////////////////////////////////////////////////////////////////////
285  ////////////////////// Signal Estimate /////////////////////////////
286  ////////////////////////////////////////////////////////////////////
287  sprintf(name, "mc_%s_%s_signal_postfit", dataName.c_str(), "tot");
288  Spectrum sSignalEst = *Spectrum::LoadFrom(fSignalEst, name);
289  TH3F* hSignalEst3D = (TH3F*)ana::ToTH3(sSignalEst, pot, kPOT,
292  TH1F* hSignalEst = (TH1F*)hSignalEst3D->ProjectionZ("hSignalEst");
293 
294  //Get Individual Systs
295  std::vector<std::unique_ptr<Spectrum>> sSystsReco;
296  std::vector<TH3F*> hSystsReco3D;
297  std::vector<TH1F*> hSystsReco1D;
298  std::vector<std::string> sSystsName = {"calibshape", "ckv", "genie", "ppfx",
299  "light", "calib"};
300 
301  for(uint i = 0; i < sSystsName.size();i++){
302  sprintf(name, "mc_%s_%s_signal_postfit", dataName.c_str(),
303  sSystsName[i].c_str());
304  sSystsReco.push_back(Spectrum::LoadFrom(fSignalEst, name));
305  hSystsReco3D.push_back((TH3F*)ana::ToTH3(*sSystsReco[i],pot,kPOT,
308  sprintf(name, "mc_%s_%s_signal_postfit_1d", dataName.c_str(),
309  sSystsName[i].c_str());
310  hSystsReco1D.push_back((TH1F*)hSystsReco3D[i]->ProjectionZ(name));
311  }//loop over systs
312 
313  //Get Fake Data Truth
314  sprintf(name, "mc_%s_analysis_%s", dataName.c_str(),chns[2].name.c_str()); // nuebar
315  Spectrum sTrueRecoSignal =
316  *Spectrum::LoadFrom(fAnalysis, name);
317  TH3F* hTrueRecoSignal3D = (TH3F*)ana::ToTH3(sTrueRecoSignal,
318  pot,kPOT,
321  TH1F* hTrueRecoSignal =
322  (TH1F*)hTrueRecoSignal3D->ProjectionZ("hTrueRecoSignal");
323 
324  //Get Nominal MC
325  sprintf(name, "mc_nominal_analysis_%s", chns[2].name.c_str()); // nuebar
326  Spectrum sTrueRecoSignal_nofit =
327  *Spectrum::LoadFrom(fAnalysis, name);
328  TH3F* hTrueRecoSignal_nofit3D = (TH3F*)ana::ToTH3(sTrueRecoSignal_nofit,
329  pot,kPOT,
332  TH1F* hTrueRecoSignal_nominal =
333  (TH1F*)hTrueRecoSignal_nofit3D->ProjectionZ("hTrueRecoSignal_nominal");
334 
335  /*
336  MakePlots(hSignalEst, hSystsReco1D, hTrueRecoSignal,
337  "Signal Events/8.09 #times 10^{20}",
338  "Reconstructed Neutrino Energy, E_{#nu} (GeV)",
339  dataName, xsecResult,"reco_signal_estimate");
340  */
341 
342  ////////////////////////////////////////////////////////////////////
343  ////////////////////// Efficiency //////////////////////////////////
344  ////////////////////////////////////////////////////////////////////
345  sprintf(name, "%s_%s_%s_%s", "mc", "nominal", "efficiency", "num_1d");
346  Spectrum eff_num_nom = *Spectrum::LoadFrom(fSysts, name);
347  sprintf(name, "%s_%s_%s_%s", "mc", "nominal", "efficiency", "denom_1d");
348  Spectrum eff_denom_nom = *Spectrum::LoadFrom(fSysts, name);
349 
350 
351  std::vector<std::unique_ptr<Spectrum>> vSystsEfficiencyNum;
352  std::vector<std::unique_ptr<Spectrum>> vSystsEfficiencyDenom;
353  std::vector<std::string> dataset_labels = {"lightdown",
354  "lightup","ckv", "calibneg",
355  "calibpos","calibshape",
356  "genie","ppfx"};
357  for(uint idataset = 0; idataset < dataset_labels.size(); idataset++){
358  if(dataset_labels[idataset] == "genie" ||
359  dataset_labels[idataset] == "ppfx"){
360 
361  for(uint i = 0; i < 100; i++){
362  sprintf(name, "%s_%s_%i_%s_%s", "mc",
363  dataset_labels[idataset].c_str(),i,
364  "efficiency","num_1d");
365  vSystsEfficiencyNum.push_back(Spectrum::LoadFrom
366  (fSysts->GetDirectory(name)));
367  sprintf(name, "%s_%s_%i_%s_%s", "mc",
368  dataset_labels[idataset].c_str(),i,
369  "efficiency","denom_1d");
370  vSystsEfficiencyDenom.push_back(Spectrum::LoadFrom
371  (fSysts->GetDirectory(name)));
372  }
373  }
374  else{
375  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[idataset].c_str(),
376  "efficiency","num_1d");
377  vSystsEfficiencyNum.push_back(Spectrum::LoadFrom
378  (fSysts->GetDirectory(name)));
379  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[idataset].c_str(),
380  "efficiency","denom_1d");
381  vSystsEfficiencyDenom.push_back(Spectrum::LoadFrom
382  (fSysts->GetDirectory(name)));
383  }
384 
385  }
386 
387  TH1F* cv_eff =(TH1F*)eff_num_nom.ToTH1(pot);
388  cv_eff->Divide((TH1F*)eff_denom_nom.ToTH1(pot));
389  std::vector<TH1F*> vSystsEff;
390  for(uint i = 0; i < vSystsEfficiencyNum.size();i++){
391  TH1F* hHolder = (TH1F*)vSystsEfficiencyNum[i]->ToTH1(pot);
392  hHolder->SetName(ana::UniqueName().c_str());
393  TH1F* hHolder2 = (TH1F*)vSystsEfficiencyDenom[i]->ToTH1(pot);
394  hHolder2->SetName(ana::UniqueName().c_str());
395  hHolder->Divide(hHolder2);
396  vSystsEff.push_back(hHolder);
397  }
398 
399  TH1F* hLightDwnEfficiency =
400  (TH1F*)vSystsEff[0]->Clone(ana::UniqueName().c_str());
401  TH1F* hLightUpEfficiency =
402  (TH1F*)vSystsEff[1]->Clone(ana::UniqueName().c_str());
403  TH1F* hCkvEfficiency =
404  (TH1F*)vSystsEff[2]->Clone(ana::UniqueName().c_str());
405  TH1F* hCalibDwnEfficiency =
406  (TH1F*)vSystsEff[3]->Clone(ana::UniqueName().c_str());
407  TH1F* hCalibUpEfficiency =
408  (TH1F*)vSystsEff[4]->Clone(ana::UniqueName().c_str());
409  TH1F* hCalibShapeEfficiency =
410  (TH1F*)vSystsEff[5]->Clone(ana::UniqueName().c_str());
411 
412  std::vector<Spectrum*> spects_genie;
413  std::vector<std::unique_ptr<Spectrum>> genie_hists;
414  for(int i = 6; i < 100+6; i++){
415  Spectrum* spect = new Spectrum(vSystsEff[i],pot,0);
416  spects_genie.push_back(spect);
417  genie_hists.push_back((std::unique_ptr<Spectrum>)spects_genie[i-6]);
418  }
419  std::vector<Spectrum*> spects_ppfx;
420  std::vector<std::unique_ptr<Spectrum>> ppfx_hists;
421  for(int i = 100+6; i < 200+6; i++){
422  Spectrum* spect = new Spectrum(vSystsEff[i],pot,0);
423  spects_ppfx.push_back(spect);
424  ppfx_hists.push_back((std::unique_ptr<Spectrum>)spects_ppfx[i-106]);
425  }
426 
427  GenieMultiverseSpectra genie_eff =
428  GenieMultiverseSpectra(100,genie_hists,true);
429  GenieMultiverseSpectra ppfx_eff =
430  GenieMultiverseSpectra(100,genie_hists,true);
431 
432  ////////////////////////////////////////////////////////////////////
433  ////////////////////// Flux ///////////////////////////////////////
434  ////////////////////////////////////////////////////////////////////
435  sprintf(name, "%s_%s_%s", "mc", "nominal", "flux");
436  Spectrum sflux_nom = *Spectrum::LoadFrom(fSysts, name);
437  TH1F* hFlux = (TH1F*)sflux_nom.ToTH1(pot);
438  hFlux->Scale(1e-4); //nu/cm^2 from nu/m^2
439  std::vector<std::unique_ptr<Spectrum>> vsSystsFlux;
440  std::vector<TH1F*> vSystsFlux;
441  for(uint i = 0; i < 100; i++){
442  sprintf(name, "%s_%s_%i_%s", "mc","ppfx",i,"flux");
443  vsSystsFlux.push_back(Spectrum::LoadFrom
444  (fSysts->GetDirectory(name)));
445  vSystsFlux.push_back((TH1F*)vsSystsFlux[i]->ToTH1(pot));
446  vSystsFlux[i]->Scale(1e-4); //nu/cm^2 from nu/m^2
447  }
448 
449  ////////////////////////////////////////////////////////////////////
450  ////////////////////// Unfolding////////////////////////////////////
451  ////////////////////////////////////////////////////////////////////
452  sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), "unfolding", "1d");
453  TH2F* hResponseMatrix =
454  (TH2F*)(Spectrum::LoadFrom(fAnalysis, name)))->ToTH2(pot;
455 
456  // TH1F* hUnfolded_SignalEst = DoUnfolding(hSignalEst, hResponseMatrix,8);
457  TH1F* hUnfolded_SignalEst = DoUnfolding(hSignalEst, hResponseMatrix,1); // try with one iteration to see if we get input back
458  std::vector<TH1F*> hUnfoldedSysts;
459  for(int i = 0; i < (int)hSystsReco1D.size(); i++){
460  // hUnfoldedSysts.push_back((TH1F*)DoUnfolding(hSystsReco1D[i],
461  // hResponseMatrix,8));
462  hUnfoldedSysts.push_back((TH1F*)DoUnfolding(hSystsReco1D[i],
463  hResponseMatrix,1));
464  }
465 
466  ////////////////////////////////////////////////////////////////////
467  ////////////////////// Calculate XSec //////////////////////////////
468  ////////////////////////////////////////////////////////////////////
469  TH1F* hXsec_cv = (TH1F*)hUnfolded_SignalEst->Clone(ana::UniqueName().c_str());
470  hXsec_cv->SetTitle("Cross Section");
471 
472  std::cout << "Here" << std::endl;
473  hXsec_cv->Divide(cv_eff);
474  hXsec_cv->Divide(hFlux);
475  hXsec_cv->Scale(1./nnucleons);
476  hXsec_cv->GetXaxis()->SetTitle("Neutrino Energy, E_{#nu} (GeV)");
477  hXsec_cv->GetYaxis()->SetTitle("#sigma (cm^{2}/nucleon)");
478  hXsec_cv->GetXaxis()->CenterTitle();
479  hXsec_cv->GetYaxis()->CenterTitle();
480 
481  std::vector<TH1F*>hXSec_systs;
482  for(int i = 0; i < (int)vSystsEff.size(); i++){
483  TH1F* hHolder =
484  (TH1F*)hUnfolded_SignalEst->Clone(ana::UniqueName().c_str());
485 
486  hHolder->Scale(1./nnucleons);
487  hHolder->Divide(vSystsEff[i]);
488 
489  if(i > 105){
490  std::cout << i << std::endl;
491  hHolder->Divide(vSystsFlux[i-106]);
492  }
493  else
494  hHolder->Divide(hFlux);
495  hHolder->SetLineColor(kBlue);
496  hXSec_systs.push_back(hHolder);
497  }
498 
499  //Get Systs Ready
500  //Order: {"lightdown","lightup","ckv", "calibneg","calibpos",
501  // "calibshape","genie","ppfx"};
502  hXSec_systs[0]->SetLineColor(kRed+3);
503  hXSec_systs[0]->SetLineStyle(7);
504  hXSec_systs[1]->SetLineColor(kRed-3);
505  hXSec_systs[1]->SetLineStyle(7);
506  hXSec_systs[2]->SetLineColor(kGreen);
507  hXSec_systs[3]->SetLineColor(kBlue-3);
508  hXSec_systs[4]->SetLineColor(kBlue+3);
509  hXSec_systs[5]->SetLineColor(kGreen);
510  hXSec_systs[5]->SetLineStyle(7);
511 
512  //Now For Multiverse Stuff
513  std::vector<Spectrum*> spects_xsec_genie;
514  std::vector<std::unique_ptr<Spectrum>> vspects_xsec_genie;
515  for(int i = 6; i < 100+6; i++){
516  Spectrum* spect = new Spectrum(hXSec_systs[i],pot,0);
517  spects_xsec_genie.push_back(spect);
518  vspects_xsec_genie.push_back((std::unique_ptr<Spectrum>)
519  spects_xsec_genie[i-6]);
520  }
521 
522  std::vector<Spectrum*> spects_xsec_ppfx;
523  std::vector<std::unique_ptr<Spectrum>> vspects_xsec_ppfx;
524  for(int i = 106; i < 200+6; i++){
525  Spectrum* spect = new Spectrum(hXSec_systs[i],pot,0);
526  spects_xsec_ppfx.push_back(spect);
527  vspects_xsec_ppfx.push_back((std::unique_ptr<Spectrum>)
528  spects_xsec_ppfx[i-106]);
529  }
530 
531  GenieMultiverseSpectra genie_xsec_syst =
532  GenieMultiverseSpectra(100,vspects_xsec_genie,true);
533  GenieMultiverseSpectra ppfx_xsec_syst =
534  GenieMultiverseSpectra(100,vspects_xsec_ppfx,true);
535 
536  TH1F* hGenie_up = (TH1F*)(genie_xsec_syst.UpperSigma())->ToTH1(pot);
537  TH1F* hGenie_dwn = (TH1F*)(genie_xsec_syst.LowerSigma())->ToTH1(pot);
538  TH1F* hPPFX_up = (TH1F*)(ppfx_xsec_syst.UpperSigma())->ToTH1(pot);
539  TH1F* hPPFX_dwn = (TH1F*)(ppfx_xsec_syst.LowerSigma())->ToTH1(pot);
540 
541  hGenie_up->SetLineColor(kOrange+3);
542  hGenie_dwn->SetLineColor(kOrange-3);
543  hPPFX_up->SetLineColor(kOrange+3);
544  hPPFX_dwn->SetLineColor(kOrange-3);
545  hPPFX_up->SetLineStyle(7);
546  hPPFX_dwn->SetLineStyle(7);
547 
548  //Order: {"lightdown","lightup","ckv", "calibneg","calibpos",
549  // "calibshape","genie","ppfx"};
550  hXSec_systs[0]->SetLineColor(kRed+3);
551  hXSec_systs[0]->SetLineStyle(7);
552  hXSec_systs[1]->SetLineColor(kRed-3);
553  hXSec_systs[1]->SetLineStyle(7);
554  hXSec_systs[2]->SetLineColor(kGreen);
555  hXSec_systs[3]->SetLineColor(kBlue-3);
556  hXSec_systs[4]->SetLineColor(kBlue+3);
557  hXSec_systs[5]->SetLineColor(kGreen);
558  hXSec_systs[5]->SetLineStyle(7);
559 
560  std::vector<TH1F*> systs_up =
561  {(TH1F*)hGenie_up->Clone(ana::UniqueName().c_str()),
562  (TH1F*)hPPFX_up->Clone(ana::UniqueName().c_str()),
563  (TH1F*)hXSec_systs[1]->Clone(ana::UniqueName().c_str()),
564  (TH1F*)hXSec_systs[4]->Clone(ana::UniqueName().c_str()),
565  (TH1F*)hXSec_systs[2]->Clone(ana::UniqueName().c_str()),
566  (TH1F*)hXSec_systs[5]->Clone(ana::UniqueName().c_str())};
567 
568  std::vector<TH1F*> systs_down =
569  {(TH1F*)hGenie_dwn->Clone(ana::UniqueName().c_str()),
570  (TH1F*)hPPFX_dwn->Clone(ana::UniqueName().c_str()),
571  (TH1F*)hXSec_systs[0]->Clone(ana::UniqueName().c_str()),
572  (TH1F*)hXSec_systs[3]->Clone(ana::UniqueName().c_str()),
573  (TH1F*)hXsec_cv->Clone(ana::UniqueName().c_str()),
574  (TH1F*)hXsec_cv->Clone(ana::UniqueName().c_str())};
575 
576  TH1F* hCVHolder = (TH1F*)hXsec_cv->Clone(ana::UniqueName().c_str());
577 
578  TGraphAsymmErrors* g =
579  PlotSystErrorBand(hCVHolder,
580  systs_up,systs_down,
581  -1,-1,1.3,false);
582 
583  //Get True Cross Section From Fake Data
584  sprintf(name, "%s_%s_%s_%s","mc",dataName.c_str(),"efficiency", "denom_1d");
585  TH1F* hTrueSignal =
586  (TH1F*)(Spectrum::LoadFrom(fAnalysis, name)))->ToTH1(pot;
587  hTrueSignal->Divide(hFlux);
588  hTrueSignal->Scale(1./nnucleons);
589  hTrueSignal->SetLineColor(kRed);
590 
591  //Propagate Uncertainties
592  //CV
593  hXsec_cv->SetMarkerColor(kBlack);
594  hXsec_cv->SetLineColor(kBlack);
595  TH1F* hUncert_tot = (TH1F*)hXsec_cv->Clone(ana::UniqueName().c_str());
596  for(int i = 1; i <= (int)hUncert_tot->GetXaxis()->GetNbins();i++){
597  double uncert = hUnfolded_SignalEst->GetBinError(i);
598  double total = hUncert_tot->GetBinContent(i);
599  double eff_val = cv_eff->GetBinContent(i);
600  double flux_val = hFlux->GetBinContent(i);
601  double selected_tot = hUnfolded_SignalEst->GetBinContent(i);
602  double err_hi = g->GetErrorYhigh(i);
603  double err_lo = g->GetErrorYlow(i);
604 
605  // double err_fluxeff = (err_hi)/total;
606  double err_fluxeff = (err_hi > err_lo)? err_hi/total : err_lo/total;
607 
608 
609  double tot_err = total*sqrt(pow(uncert/selected_tot,2) +
610  pow(err_fluxeff,2));
611 
612  hUncert_tot->SetBinError(i,tot_err);
613  std::cout << i << "\t" << total << "\t" <<
614  uncert << "\t" << selected_tot << "\t" <<
615  err_fluxeff << "\t" <<
616  sqrt(pow(uncert/selected_tot,2) +
617  pow(err_fluxeff,2)) << std::endl;
618 
619  }
620 
621 
622  TCanvas* cResult = new TCanvas("cResult", "cResult");
623  hUncert_tot->GetXaxis()->SetRangeUser(1.0,10.0);
624  //hXsec_cv->Draw();
625  hUncert_tot->SetMarkerStyle(20);
626  hUncert_tot->Draw();
627  g->Draw("e2 same");
628  hXSec_systs[0]->Draw("hist same");
629  hXSec_systs[1]->Draw("hist same");
630  hXSec_systs[2]->Draw("hist same");
631  hXSec_systs[3]->Draw("hist same");
632  hXSec_systs[4]->Draw("hist same");
633  hXSec_systs[5]->Draw("hist same");
634  hXSec_systs[5]->Draw("hist same");
635  //for(uint j = 0; j < hXSec_systs.size(); j++){
636  //hXSec_systs[j]->Draw("hist same");
637  //}
638  hGenie_up->Draw("hist same");
639  hGenie_dwn->Draw("hist same");
640  hPPFX_up->Draw("hist same");
641  hPPFX_dwn->Draw("hist same");
642  hTrueSignal->Draw("hist same");
643  hUncert_tot->Draw("same");
644  TLegend* leg_total = ana::AutoPlaceLegend(0.3,0.3);
645  leg_total->AddEntry(hTrueSignal, "Truth", "l");
646  leg_total->AddEntry(hUncert_tot, "FakeData", "p");
647  leg_total->AddEntry(hXSec_systs[3], "Calibration -1 #sigma", "l");
648  leg_total->AddEntry(hXSec_systs[4], "Calibration +1 #sigma", "l");
649  leg_total->AddEntry(hXSec_systs[0], "Light Level -1 #sigma", "l");
650  leg_total->AddEntry(hXSec_systs[1], "Light Level +1 #sigma", "l");
651  leg_total->AddEntry(hXSec_systs[2], "Cherenkov Variation", "l");
652  leg_total->AddEntry(hXSec_systs[5], "Calibration Shape", "l");
653  leg_total->AddEntry(hPPFX_dwn, "Flux -1 #sigma", "l");
654  leg_total->AddEntry(hPPFX_up, "Flux +1 #sigma", "l");
655  leg_total->AddEntry(hGenie_dwn, "#nu-Interaction -1 #sigma", "l");
656  leg_total->AddEntry(hGenie_up, "#nu-Interaction +1 #sigma", "l");
657  leg_total->Draw();
658  sprintf(name, "%s%s_%s_%s_%s.png","AnalysisPlots/", dataName.c_str(),
659  "xsec", "results_all_systs",xsecResult.c_str());
660  cResult->SaveAs(name);
661 
662 
663  TGraph *gr_nue280=(TGraph*)fGenie->Get("nu_e_bar_C12/tot_cc");
664  TSpline3* spline_nue280=new TSpline3("mysplinenue280",gr_nue280);
665  TH1F *Xnue280=new TH1F("Xnue280","",100,0,10);
666 
667  for( int ibin=1; ibin<Xnue280->GetNbinsX()+1; ++ibin ){
668  double binctr=Xnue280->GetBinCenter(ibin);
669  double xnue280=spline_nue280->Eval(binctr)*pow(10,-38)/12.;
670  Xnue280->SetBinContent(ibin,xnue280);
671  Xnue280->SetBinError(ibin,0.);
672  }
673 
674  Xnue280->SetLineWidth(2);
675  Xnue280->SetLineColor(4);
676 
677  std::cout << "here" << std::endl;
678  TCanvas* c3 = new TCanvas("c3", "cResult");
679  hUncert_tot->GetXaxis()->SetRangeUser(1.0,10.0);
680  hUncert_tot->SetMarkerStyle(20);
681  hUncert_tot->Draw();
682  Xnue280->Draw("l same");
683  //g->Draw("e2 same");
684  hTrueSignal->Draw("hist same");
685  hUncert_tot->Draw("same");
686  //Xnue280->Draw("l same");
687  sprintf(name, "%s%s_%s_%s_%s.png","AnalysisPlots/", dataName.c_str(),
688  "xsec", "results_truth_compare",xsecResult.c_str());
689  c3->SaveAs(name);
690 
691 
692  //Calculate Relative Uncertainties
693  //CV
694  TH1F* hRelUnc_tot = (TH1F*)hUncert_tot->Clone(ana::UniqueName().c_str());
695  for(int i = 0; i <= (int)hRelUnc_tot->GetXaxis()->GetNbins();i++){
696  double uncert = hRelUnc_tot->GetBinError(i);
697  double total = hRelUnc_tot->GetBinContent(i);
698  hRelUnc_tot->SetBinContent(i,uncert/total);
699  }
700 
701  std::cout << "Here" << std::endl;
702  std::vector<TH1F*> hRelUncert;
703  for(int i = 0; i < (int)systs_up.size();i++){
704  TH1F* hHolder = (TH1F*)systs_up[i]->Clone(ana::UniqueName().c_str());
705  TH1F* hHolderLo = (TH1F*)systs_down[i]->Clone(ana::UniqueName().c_str());
706 
707  for(int ibin = 0; ibin <= hHolder->GetXaxis()->GetNbins(); ibin++){
708  double val_up = hHolder->GetBinContent(ibin);
709  double val_down = hHolderLo->GetBinContent(ibin);
710  double val_cv = hXsec_cv->GetBinContent(ibin);
711 
712  // double uncert = 0.5*( fabs(val_up - val_cv) + fabs(val_down - val_cv));
713  double uncert = std::max(fabs(val_up - val_cv), fabs(val_down - val_cv));
714  hHolder->SetBinContent(ibin,uncert/val_cv);
715  }
716  hRelUncert.push_back(hHolder);
717  }
718  std::cout << "Here 2" << std::endl;
719 
720  TH1F* hRelUnc_fit = (TH1F*)hXsec_cv->Clone(ana::UniqueName().c_str());
721  for(int i = 0; i <= (int)hRelUnc_tot->GetXaxis()->GetNbins();i++){
722  double uncert = hRelUnc_fit->GetBinError(i);
723  double total = hRelUnc_fit->GetBinContent(i);
724  hRelUnc_fit->SetBinContent(i,uncert/total);
725  }
726 
727  std::cout << "Here 3" << std::endl;
728 
729  //Same Order as vector used to generate TGraphAsmyErrors
730  hRelUncert[5]->SetLineColor(kGreen+2); //calibshape
731  hRelUncert[4]->SetLineColor(kOrange-3); //ckv
732  hRelUncert[0]->SetLineColor(kRed - 2); //genie
733  hRelUncert[1]->SetLineColor(kCyan-3); //ppfx
734  hRelUncert[2]->SetLineColor(kBlue-4); //light
735  hRelUncert[3]->SetLineColor(kMagenta+1); //calib
736  hRelUnc_fit->SetLineColor(kBlack);
737  hRelUnc_fit->SetLineStyle(7);
738  hRelUnc_tot->SetLineColor(kBlack);
739 
740  TCanvas* c2 = new TCanvas(ana::UniqueName().c_str(),"c2");
741  hRelUnc_tot->GetYaxis()->SetTitle("Relative Uncertainty");
742  hRelUnc_tot->SetTitle(" ");
743  hRelUnc_tot->GetYaxis()->SetRangeUser(0,1.0);
744  hRelUnc_tot->GetXaxis()->SetRangeUser(1,10);
745  hRelUnc_tot->Draw("hist");
746  hRelUncert[5]->Draw("hist same");
747  hRelUncert[4]->Draw("hist same");
748  hRelUncert[3]->Draw("hist same");
749  hRelUncert[2]->Draw("hist same");
750  hRelUncert[1]->Draw("hist same");
751  hRelUncert[0]->Draw("hist same");
752  hRelUnc_fit->Draw("hist same");
753  hRelUnc_tot->Draw("hist same");
754  TLegend* leg_systs = ana::AutoPlaceLegend(0.3,0.3);
755  leg_systs->AddEntry(hRelUnc_tot, "Total", "l");
756  leg_systs->AddEntry(hRelUnc_fit, "Fit", "l");
757  leg_systs->AddEntry(hRelUncert[3], "Calibration Normalization", "l");
758  leg_systs->AddEntry(hRelUncert[2], "Light Level", "l");
759  leg_systs->AddEntry(hRelUncert[1], "Flux", "l");
760  leg_systs->AddEntry(hRelUncert[0], "#nu-Interaction", "l");
761  leg_systs->AddEntry(hRelUncert[4], "Cherenkov Variation", "l");
762  leg_systs->AddEntry(hRelUncert[5], "Calibration Shape", "l");
763  leg_systs->Draw();
764  sprintf(name, "%s%s_%s_%s_%s.png","AnalysisPlots/", dataName.c_str(),
765  "xsec", "relative_uncertantiy",xsecResult.c_str());
766  c2->SaveAs(name);
767 
768 
769 
770 
771 
772 }
void Simulation()
Definition: tools.h:16
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
void MakePlots(TH1F *hFitted1D, std::vector< TH1F * > hSysts, TH1F *hTruth1D, std::string ytitle, std::string xtitle, std::string dataName, std::string xsecResult, std::string varName, bool calcCov=false)
const double pot
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const ana::Binning eelecbins
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:209
const SelDef chns[knumchns]
const std::string sGenie
const ana::Binning costhetabins
T sqrt(T number)
Definition: d0nt_math.hpp:156
void makeXSecPlots1D(std::string dataName="nominal", std::string xsecResult="energy")
constexpr T pow(T x)
Definition: pow.h:75
const std::string sSignalEst
const std::string sDir
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
c2
Definition: demo5.py:33
const double nnucleons
void PrintChiSq(TString str)
Definition: MakePlots.C:246
const Binning enubins
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:607
TH1F * DoUnfolding(TH1F *hMeasured, TH2F *hResponse, int iter)
Regular histogram.
Definition: Utilities.h:44
OStream cout
Definition: OStream.cxx:6
const std::string sAnalysis
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TSpline3 * spline_nue280
Definition: Xnue_genie.C:8
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal) const
c1
Definition: demo5.py:24
void MakePlots1D(TH1F *hData, TH1F *hFit, TH1F *hMC, std::string ytitle, std::string xtitle, std::string dataName, std::string xsecResult, std::string varName)
const Spectrum * LowerSigma(BandOptions opt=kBandFromNominal) const
std::string dataName
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
const double kAna2018FHCPOT
Definition: Exposures.h:207
const std::string sSysts
TPad * pad2
Definition: analysis.C:13
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
std::string name
Float_t e
Definition: plot.C:35
TH1F * Xnue280
Definition: Xnue_genie.C:14
TH3 * ToTH3(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, const Binning &binsz, ana::EBinType bintype)
Same as ToTH2, but with 3 dimensions.
Definition: Utilities.cxx:303
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
TGraphAsymmErrors * PlotSystErrorBand(TH1F *&nom, std::vector< TH1F * > &ups, std::vector< TH1F * > &dns, int col, int errCol, float headroom, bool newaxis)
TGraph * gr_nue280
Definition: Xnue_genie.C:7
TPad * pad1
Definition: analysis.C:13
unsigned int uint