Functions | Variables
makeXSecPlots2D.C File Reference
#include "NDAna/nuebarcc_inc/NuebarCCIncVars.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncCuts.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncExtra.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncTemplateFit.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncCrossSectionFunctions.h"
#include "RooUnfoldResponse.h"
#include "RooUnfoldBayes.h"
#include "Utilities/rootlogon.C"
#include "CAFAna/XSec/TargetCount.h"
#include "CAFAna/Analysis/Exposures.h"

Go to the source code of this file.

Functions

void MakePlots1D (TH2F *hData, TH2F *hFit, TH2F *hMC, std::string ytitle, std::string xtitle, std::string dataName, std::string xsecResult, std::string varName)
 
void MakePlots (TH2F *hFitted2D, std::vector< TH2F * > hSysts, TH2F *hTruth2D, std::string ytitle, std::string xtitle, std::string dataName, std::string xsecResult, std::string varName, bool calcCov=false)
 
void makeXSecPlots2D (std::string dataName="nominal", std::string xsecResult="double")
 

Variables

const double pot = kAna2018FHCPOT
 
const double nnucleons = 3.02225e+31
 
const std::string sDir = "/nova/ana/users/ddoyle/NuebarCrossSection_020319/"
 
const std::string sSignalEst = "CrossSection_TemplateFitResults_fakedata_development_stat.root"
 
const std::string sAnalysis = "SpectraNominal/fCrossSection_Nominal_FakeData_prongcvn_new.root"
 
const std::string sSysts = "SpectraSystematics/fCrossSection_EfficiencyFlux_Systs.root"
 

Function Documentation

void MakePlots ( TH2F *  hFitted2D,
std::vector< TH2F * >  hSysts,
TH2F *  hTruth2D,
std::string  ytitle,
std::string  xtitle,
std::string  dataName,
std::string  xsecResult,
std::string  varName,
bool  calcCov = false 
)

Definition at line 137 of file makeXSecPlots2D.C.

References ana::AutoPlaceLegend(), demo5::c1, demo5::c2, make_true_q0q3_plots::caption, om::cout, allTimeWatchdog::endl, MECModelEnuComparisons::i, make_syst_table_plots::ibin, makeTrainCVSamples::int, kBlue, kGreen, kOrange, kRed, MECModelEnuComparisons::leg, ana::ProjectionY(), string, produceInfoGainMatrix::total, and ana::UniqueName().

Referenced by makeXSecPlots2D().

141 {
142 
143  for(int xbin = 0; xbin <= hTruth2D->GetXaxis()->GetNbins(); xbin++){
144 
145  std::stringstream low_X;
146  low_X << std::fixed << std::setprecision(2) <<
147  hTruth2D->GetXaxis()->GetBinLowEdge(xbin);
148  std::stringstream hi_X;
149  hi_X << std::fixed << std::setprecision(2) <<
150  hTruth2D->GetXaxis()->GetBinUpEdge(xbin);
151 
152  std::string caption = low_X.str() + " #leq cos #theta_{e} < "+
153  hi_X.str();
154 
155  TH1F* hFitted =
156  (TH1F*)hFitted2D->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
157  TH1F* hTruth =
158  (TH1F*)hTruth2D->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
159 
160  std::cout << hSysts.size() << std::endl;
161 
162  hFitted->SetLineColor(kBlack);
163  hFitted->SetMarkerColor(kBlack);
164  hFitted->SetMarkerStyle(20);
165  hFitted->GetXaxis()->SetRangeUser(1.0,10.0);
166  hFitted->SetTitle(caption.c_str());
167  hTruth->SetLineColor(kRed);
168  hTruth->SetMarkerColor(kRed);
169 
170 
171  hFitted->GetYaxis()->SetTitle(ytitle.c_str());
172  hFitted->GetXaxis()->SetTitle(xtitle.c_str());
173 
174  std::vector<TH1F*> hRelUncert;
175  for(int i = 0; i < (int)hSysts.size();i++){
176  TH1F* hHolder =
177  (TH1F*)hSysts[i]->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
178  for(int ibin = 0; ibin <= hHolder->GetXaxis()->GetNbins(); ibin++){
179  double uncert = hHolder->GetBinError(ibin);
180  double total = hHolder->GetBinContent(ibin);
181 
182  if(total == 0)
183  hHolder->SetBinContent(ibin,0);
184  else
185  hHolder->SetBinContent(ibin,uncert/total);
186  }
187  hRelUncert.push_back(hHolder);
188  }
189 
190  TH1F* hHoldera = (TH1F*)hFitted->Clone(ana::UniqueName().c_str());
191  for(int ibin = 0; ibin <= hHoldera->GetXaxis()->GetNbins(); ibin++){
192  double uncert = hHoldera->GetBinError(ibin);
193  double total = hHoldera->GetBinContent(ibin);
194 
195  if(total == 0)
196  hHoldera->SetBinContent(ibin,0);
197  else
198  hHoldera->SetBinContent(ibin,uncert/total);
199  }
200 
201  hRelUncert.push_back(hHoldera);
202 
203  hRelUncert[0]->SetLineColor(kGreen+2); //calibshape
204  hRelUncert[1]->SetLineColor(kOrange-3); //ckv
205  hRelUncert[2]->SetLineColor(kRed - 2); //genie
206  hRelUncert[3]->SetLineColor(kCyan-3); //ppfx
207  hRelUncert[4]->SetLineColor(kBlue-4); //light
208  hRelUncert[5]->SetLineColor(kMagenta+1); //calib
209  hRelUncert[6]->SetLineColor(kBlack); //tot
210 
211  hSysts[0]->SetLineColor(kGreen+2); //calibshape
212  hSysts[1]->SetLineColor(kOrange-3); //ckv
213  hSysts[2]->SetLineColor(kRed - 2); //genie
214  hSysts[3]->SetLineColor(kCyan-3); //ppfx
215  hSysts[4]->SetLineColor(kBlue-4); //light
216  hSysts[5]->SetLineColor(kMagenta+1); //calib
217  hSysts[0]->SetMarkerColor(kGreen+2); //calibshape
218  hSysts[1]->SetMarkerColor(kOrange-3); //ckv
219  hSysts[2]->SetMarkerColor(kRed - 2); //genie
220  hSysts[3]->SetMarkerColor(kCyan-3); //ppfx
221  hSysts[4]->SetMarkerColor(kBlue-4); //light
222  hSysts[5]->SetMarkerColor(kMagenta+1); //calib
223 
224  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str(),"c1");
225  hFitted->SetMinimum(0);
226  hFitted->Draw();
227  hTruth->Draw("hist same");
228  TLegend* leg = ana::AutoPlaceLegend(0.3,0.3);
229  leg->AddEntry(hFitted, "Fitted Signal Estimate", "pl");
230  leg->AddEntry(hTruth,"Truth", "l");
231  leg->Draw();
232  char outname[50];
233  sprintf(outname, "%s%s_%s_%s_%i.png","AnalysisPlots/", dataName.c_str(),
234  varName.c_str(), xsecResult.c_str(),xbin);
235  c1->SaveAs(outname);
236  c1->Close();
237 
238 
239  TCanvas* c2 = new TCanvas(ana::UniqueName().c_str(),"c2");
240  hRelUncert[6]->GetYaxis()->SetTitle("Relative Uncertainty");
241  hRelUncert[6]->GetYaxis()->SetRangeUser(0,0.3);
242  hRelUncert[6]->Draw("hist");
243  hRelUncert[5]->Draw("hist same");
244  hRelUncert[4]->Draw("hist same");
245  hRelUncert[3]->Draw("hist same");
246  hRelUncert[2]->Draw("hist same");
247  hRelUncert[1]->Draw("hist same");
248  hRelUncert[0]->Draw("hist same");
249  hRelUncert[6]->Draw("hist same");
250  TLegend* leg_systs = ana::AutoPlaceLegend(0.3,0.3);
251  leg_systs->AddEntry(hRelUncert[6], "Total", "l");
252  leg_systs->AddEntry(hRelUncert[5], "Calibration Normalization", "l");
253  leg_systs->AddEntry(hRelUncert[4], "Light Level", "l");
254  leg_systs->AddEntry(hRelUncert[3], "Flux", "l");
255  leg_systs->AddEntry(hRelUncert[2], "#nu-Interaction", "l");
256  leg_systs->AddEntry(hRelUncert[1], "Cherenkov Variation", "l");
257  leg_systs->AddEntry(hRelUncert[0], "Calibration Shape", "l");
258  leg_systs->Draw();
259  sprintf(outname, "%s%s_%s_%s_%s_%i.png","AnalysisPlots/", dataName.c_str(),
260  varName.c_str(), "relative_uncertantiy",xsecResult.c_str(),xbin);
261  c2->SaveAs(outname);
262  }//loop over xbins
263 }
kOrange
enum BeamMode kOrange
kRed
enum BeamMode kRed
make_syst_table_plots.xtitle
string xtitle
Definition: make_syst_table_plots.py:204
demo5.c2
c2
Definition: demo5.py:33
make_syst_table_plots.ibin
ibin
Definition: make_syst_table_plots.py:246
make_true_q0q3_plots.caption
caption
Definition: make_true_q0q3_plots.py:404
om::cout
OStream cout
Definition: OStream.cxx:6
ana::ProjectionY
Eigen::ArrayXd ProjectionY(const Eigen::MatrixXd &mat)
Helper for WeightingVariable.
Definition: ReweightableSpectrum.cxx:105
demo5.c1
c1
Definition: demo5.py:24
dataName
std::string dataName
Definition: make_starplots.C:22
ana::AutoPlaceLegend
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:717
kGreen
enum BeamMode kGreen
kBlue
enum BeamMode kBlue
ana::UniqueName
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
string
enum BeamMode string
void MakePlots1D ( TH2F *  hData,
TH2F *  hFit,
TH2F *  hMC,
std::string  ytitle,
std::string  xtitle,
std::string  dataName,
std::string  xsecResult,
std::string  varName 
)

Definition at line 25 of file makeXSecPlots2D.C.

References ana::AutoPlaceLegend(), demo5::c1, make_true_q0q3_plots::caption, kBlue, kRed, MECModelEnuComparisons::leg, pad1, pad2, PrintChiSq(), Simulation(), string, art::to_string(), and ana::UniqueName().

Referenced by makeXSecPlots2D().

29 {
30 
31  for(int xbin = 0; xbin <= hData->GetXaxis()->GetNbins(); xbin++){
32 
33  std::stringstream low_X;
34  low_X << std::fixed << std::setprecision(2) <<
35  hData->GetXaxis()->GetBinLowEdge(xbin);
36  std::stringstream hi_X;
37  hi_X << std::fixed << std::setprecision(2) <<
38  hData->GetXaxis()->GetBinUpEdge(xbin);
39 
40  std::string caption = low_X.str() + " #leq cos #theta_{e} < "+
41  hi_X.str();
42 
43 
44  TH1F* hD = (TH1F*)hData->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
45  TH1F* hF = (TH1F*)hFit->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
46  TH1F* hN = (TH1F*)hMC->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
47 
48  hD->SetMarkerStyle(20);
49  hF->SetMarkerStyle(22);
50  hF->SetMarkerColor(kBlue);
51  hF->SetLineColor(kBlue);
52  hN->SetLineColor(kRed);
53  hF->SetTitle("");
54  hN->SetTitle("");
55  hD->GetYaxis()->SetTitle(ytitle.c_str());
56  hD->GetXaxis()->SetTitle(xtitle.c_str());
57 
58 
59  double chi_fit = 0;
60  int ndf_fit = 0, igood_fit = 0;
61  hD->Chi2TestX(hF,chi_fit, ndf_fit, igood_fit,"WW");
62  double chi_nom = 0;
63  int ndf_nom = 0, igood_nom = 0;
64  hD->Chi2TestX(hN,chi_nom, ndf_nom, igood_nom,"WW");
65 
66  std::string caption_chi =
67  "Fitted #chi^{2}: " + std::to_string(chi_fit) + "\n"
68  + "Nominal #chi^{2}: " + std::to_string(chi_nom);
69 
70 
71  TH1F* hratio_fit = (TH1F*)hF->Clone(ana::UniqueName().c_str());
72  TH1F* hratio_nom = (TH1F*)hN->Clone(ana::UniqueName().c_str());
73 
74  hD->GetXaxis()->SetRangeUser(1.0,10.0);
75  hF->GetXaxis()->SetRangeUser(1.0,10.0);
76  hN->GetXaxis()->SetRangeUser(1.0,10.0);
77  hratio_fit->GetXaxis()->SetRangeUser(1.0,10.0);
78  hratio_nom->GetXaxis()->SetRangeUser(1.0,10.0);
79 
80  hratio_fit->Divide(hD);
81  hratio_nom->Divide(hD);
82 
83  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str(), "c1");
84  c1->SetLeftMargin(0.10);
85  c1->SetRightMargin(0.10);
86  TPad *pad1 = new TPad(ana::UniqueName().c_str(), "pad1",0,0,1,1);
87  TPad *pad2 = new TPad(ana::UniqueName().c_str(), "pad2",0,0,1,1);
88  pad1->SetFillStyle(0);
89  pad2->SetFillStyle(0);
90  pad1->SetBottomMargin(0.30);
91  pad2->SetTopMargin(1-0.30);
92  pad2->SetGridx();
93  pad2->SetGridy();
94  c1->cd();
95  pad1->Draw();
96  pad1->cd();
97  hD->SetTitle(caption.c_str());
98  hD->GetXaxis()->SetLabelColor(kWhite);
99  hD->GetYaxis()->SetLabelSize(0.03);
100  hD->GetYaxis()->SetTitleSize(0.035);
101  hD->SetTitle(caption.c_str());
102  hD->Draw();
103  PrintChiSq(caption_chi.c_str());
104  hD->Draw("same");
105  hF->Draw("same");
106  hN->Draw("hist same");
107  TLegend* leg = ana::AutoPlaceLegend(0.3,0.3);
108  leg->AddEntry(hD, "Fake Data", "p");
109  leg->AddEntry(hN, "Nominal MC", "l");
110  leg->AddEntry(hF, "Fitted MC", "p");
111  leg->Draw();
112  Simulation();
113  gPad->RedrawAxis();
114  c1->cd();
115  pad2->Draw();
116  pad2->cd();
117  hratio_fit->GetYaxis()->SetTitle("MC/Data");
118  hratio_fit->GetXaxis()->SetTitle(xtitle.c_str());
119  hratio_fit->GetXaxis()->CenterTitle();
120  hratio_fit->GetYaxis()->SetLabelSize(0.03);
121  hratio_fit->GetYaxis()->SetTitleSize(0.035);
122  hratio_fit->GetYaxis()->SetRangeUser(0.4,1.5);
123  hratio_fit->Draw("same");
124  hratio_nom->Draw("hist same");
125  pad2->Draw("same");
126  gPad->RedrawAxis();
127  c1->cd();
128  c1->Update();
129  char outname[50];
130  sprintf(outname, "%s%s_%s_%i.png", "AnalysisPlots/",dataName.c_str(),
131  varName.c_str(), xbin);
132  c1->SaveAs(outname);
133  c1->Close();
134  } //loop over xbins
135 }
Simulation
void Simulation()
Definition: tools.h:16
kRed
enum BeamMode kRed
make_syst_table_plots.xtitle
string xtitle
Definition: make_syst_table_plots.py:204
make_true_q0q3_plots.caption
caption
Definition: make_true_q0q3_plots.py:404
PrintChiSq
void PrintChiSq(TString str)
Definition: MakePlots.C:246
demo5.c1
c1
Definition: demo5.py:24
dataName
std::string dataName
Definition: make_starplots.C:22
ana::AutoPlaceLegend
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:717
pad2
TPad * pad2
Definition: analysis.C:13
art::to_string
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
kBlue
enum BeamMode kBlue
ana::UniqueName
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
pad1
TPad * pad1
Definition: analysis.C:13
string
enum BeamMode string
void makeXSecPlots2D ( std::string  dataName = "nominal",
std::string  xsecResult = "double" 
)

Definition at line 264 of file makeXSecPlots2D.C.

References ana::AutoPlaceLegend(), demo5::c2, chisquared::c3, make_true_q0q3_plots::caption, ana::chns, plot_validation_datamc::Clone(), nuebarccinc::costhetabins, om::cout, dataName, DoUnfolding(), e, nuebarccinc::eelecbins, allTimeWatchdog::endl, ana::enubins, stan::math::fabs(), MECModelEnuComparisons::g, MECModelEnuComparisons::i, make_syst_table_plots::ibin, makeTrainCVSamples::int, Integral(), calib::j, ana::kBinContent, kBlue, kGreen, kOrange, ana::kPOT, kRed, ana::Spectrum::LoadFrom(), ana::GenieMultiverseSpectra::LowerSigma(), MakePlots(), MakePlots1D(), std::max(), ana::SelDef::name, nnucleons, PlotSystErrorBand(), pot, cet::pow(), ana::ProjectionX(), ana::ProjectionY(), sAnalysis, sDir, std::sqrt(), sSignalEst, sSysts, string, produceInfoGainMatrix::total, ana::Spectrum::ToTH1(), ana::ToTH2(), ana::ToTH3(), ana::UniqueName(), ana::GenieMultiverseSpectra::UpperSigma(), and febshutoff_auto::val.

266 {
267  TFile *fSignalEst = new TFile((sDir+sSignalEst).c_str(),"read");
268  TFile* fAnalysis = new TFile((sDir+sAnalysis).c_str(),"read");
269  TFile* fSysts = new TFile((sDir+sSysts).c_str(),"read");
270 
271  ////////////////////////////////////////////////////////////////////
272  ////////////////////// Compare Fit Results /////////////////////////
273  ////////////////////////////////////////////////////////////////////
274  char name[50];
275  sprintf(name, "mc_%s_%s_%s_postfit",dataName.c_str(),
276  "tot", chns[0].name.c_str());
277  Spectrum sFitTotal = *Spectrum::LoadFrom(fSignalEst, name);
278  TH3F* hFitTotal3D = (TH3F*)ana::ToTH3(sFitTotal,pot,kPOT,
279  costhetabins,eelecbins,
280  enubins,kBinContent);
281  TH2F* hFitTotal2D = (TH2F*)hFitTotal3D->Project3D("yx");
282  hFitTotal2D->SetName("hFitTotal2D");
283 
284  //Get Total Fake Data
285  sprintf(name, "mc_%s_%s_%s", dataName.c_str(), "analysis", "mc");
286  Spectrum sDataTotal = *Spectrum::LoadFrom(fAnalysis, name);
287  TH3F* hDataTotal3D = (TH3F*)ana::ToTH3(sDataTotal,pot,kPOT,
288  costhetabins,eelecbins,
289  enubins, kBinContent);
290  TH2F* hDataTotal2D = (TH2F*)hDataTotal3D->Project3D("yx");
291  hDataTotal2D->SetName("hDataTotal1D");
292 
293  //Get Total Nominal MC
294  sprintf(name, "mc_%s_%s_%s", "nominal", "analysis", "mc");
295  Spectrum sMCTotal = *Spectrum::LoadFrom(fAnalysis, name);
296  TH3F* hMCTotal3D = (TH3F*)ana::ToTH3(sMCTotal,pot,kPOT,
297  costhetabins,eelecbins,
298  enubins, kBinContent);
299  TH2F* hMCTotal2D = (TH2F*)hMCTotal3D->Project3D("yx");
300  hMCTotal2D->SetName("hMCTotal2D");
301 
302  MakePlots1D(hDataTotal2D,hFitTotal2D,hMCTotal2D,
303  "Events/8.09 #times 10^{20} POT",
304  "Reconstructed Electron Energy, E_{e} (GeV)",
305  dataName,xsecResult, "total_fit_comparison");
306 
307  ////////////////////////////////////////////////////////////////////
308  ////////////////////// Signal Estimate /////////////////////////////
309  ////////////////////////////////////////////////////////////////////
310  sprintf(name, "mc_%s_%s_signal_postfit", dataName.c_str(), "tot");
311  Spectrum sSignalEst = *Spectrum::LoadFrom(fSignalEst, name);
312  TH3F* hSignalEst3D = (TH3F*)ana::ToTH3(sSignalEst, pot, kPOT,
313  costhetabins,eelecbins,
314  enubins,kBinContent);
315  TH2F* hSignalEst = (TH2F*)hSignalEst3D->Project3D("yx");
316  hSignalEst->SetName("hSignalEst");
317 
318  //Get Individual Systs
319  std::vector<std::unique_ptr<Spectrum>> sSystsReco;
320  std::vector<TH3F*> hSystsReco3D;
321  std::vector<TH2F*> hSystsReco2D;
322  std::vector<std::string> sSystsName = {"calibshape", "ckv", "genie", "ppfx",
323  "light", "calib"};
324 
325  for(uint i = 0; i < sSystsName.size();i++){
326  sprintf(name, "mc_%s_%s_signal_postfit", dataName.c_str(),
327  sSystsName[i].c_str());
328  sSystsReco.push_back(Spectrum::LoadFrom(fSignalEst, name));
329  hSystsReco3D.push_back((TH3F*)ana::ToTH3(*sSystsReco[i],pot,kPOT,
330  costhetabins,eelecbins,
331  enubins,kBinContent));
332  sprintf(name, "mc_%s_%s_signal_postfit_1d", dataName.c_str(),
333  sSystsName[i].c_str());
334  hSystsReco2D.push_back((TH2F*)hSystsReco3D[i]->Project3D("yx"));
335  hSystsReco2D[i]->SetName(name);
336  }//loop over systs
337 
338 
339  //Get Fake Data Truth
340  sprintf(name, "mc_%s_analysis_%s", dataName.c_str(),chns[2].name.c_str()); //nuebar
341  Spectrum sTrueRecoSignal =
342  *Spectrum::LoadFrom(fAnalysis, name);
343  TH3F* hTrueRecoSignal3D = (TH3F*)ana::ToTH3(sTrueRecoSignal,
344  pot,kPOT,
345  costhetabins,eelecbins,
346  enubins,kBinContent);
347  TH2F* hTrueRecoSignal =
348  (TH2F*)hTrueRecoSignal3D->Project3D("yx");
349  hTrueRecoSignal->SetName("hTrueRecoSignal");
350 
351  //Get Nominal MC
352  sprintf(name, "mc_nominal_analysis_%s", chns[2].name.c_str());
353  Spectrum sTrueRecoSignal_nofit =
354  *Spectrum::LoadFrom(fAnalysis, name);
355  TH3F* hTrueRecoSignal_nofit3D = (TH3F*)ana::ToTH3(sTrueRecoSignal_nofit,
356  pot,kPOT,
357  costhetabins,eelecbins,
358  enubins,kBinContent);
359  TH2F* hTrueRecoSignal_nominal =
360  (TH2F*)hTrueRecoSignal_nofit3D->Project3D("yx");
361  hTrueRecoSignal_nominal->SetName("hTrueRecoSignal_nominal");
362 
363  MakePlots(hSignalEst, hSystsReco2D, hTrueRecoSignal,
364  "Signal Events/8.09 #times 10^{20}",
365  "Reconstructed Electron Energy, E_{e} (GeV)",
366  dataName, xsecResult,"reco_signal_estimate");
367 
368 
369  ////////////////////////////////////////////////////////////////////
370  ////////////////////// Efficiency //////////////////////////////////
371  ////////////////////////////////////////////////////////////////////
372  sprintf(name, "%s_%s_%s_%s", "mc", "nominal", "efficiency", "num_2d");
373  Spectrum eff_num_nom = *Spectrum::LoadFrom(fSysts, name);
374  sprintf(name, "%s_%s_%s_%s", "mc", "nominal", "efficiency", "denom_2d");
375  Spectrum eff_denom_nom = *Spectrum::LoadFrom(fSysts, name);
376 
377  std::vector<std::unique_ptr<Spectrum>> vSystsEfficiencyNum;
378  std::vector<std::unique_ptr<Spectrum>> vSystsEfficiencyDenom;
379  std::vector<std::string> dataset_labels = {"lightdown",
380  "lightup","ckv", "calibneg",
381  "calibpos","calibshape",
382  "genie","ppfx"};
383 
384  for(uint idataset = 0; idataset < dataset_labels.size(); idataset++){
385  if(dataset_labels[idataset] == "genie" ||
386  dataset_labels[idataset] == "ppfx"){
387 
388  for(uint i = 0; i < 100; i++){
389  sprintf(name, "%s_%s_%i_%s_%s", "mc",
390  dataset_labels[idataset].c_str(),i,
391  "efficiency","num_2d");
392  vSystsEfficiencyNum.push_back(Spectrum::LoadFrom
393  (fSysts->GetDirectory(name)));
394  sprintf(name, "%s_%s_%i_%s_%s", "mc",
395  dataset_labels[idataset].c_str(),i,
396  "efficiency","denom_2d");
397  vSystsEfficiencyDenom.push_back(Spectrum::LoadFrom
398  (fSysts->GetDirectory(name)));
399  }
400  }
401  else{
402  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[idataset].c_str(),
403  "efficiency","num_2d");
404  vSystsEfficiencyNum.push_back(Spectrum::LoadFrom
405  (fSysts->GetDirectory(name)));
406  sprintf(name, "%s_%s_%s_%s", "mc", dataset_labels[idataset].c_str(),
407  "efficiency","denom_2d");
408  vSystsEfficiencyDenom.push_back(Spectrum::LoadFrom
409  (fSysts->GetDirectory(name)));
410  }
411  }
412 
413  TH2F* cv_eff =(TH2F*)ana::ToTH2(eff_num_nom,pot,kPOT,
414  costhetabins,eelecbins,kBinContent);
415  cv_eff->Divide((TH2F*)ana::ToTH2(eff_denom_nom,pot,kPOT,
416  costhetabins,eelecbins,kBinContent));
417 
418  std::vector<TH2F*> vSystsEff;
419  for(uint i = 0; i < vSystsEfficiencyNum.size();i++){
420  TH2F* hHolder = (TH2F*)ana::ToTH2(*vSystsEfficiencyNum[i],pot,kPOT,
421  costhetabins,eelecbins,kBinContent);
422  hHolder->SetName(ana::UniqueName().c_str());
423  TH2F* hHolder2 = (TH2F*)ana::ToTH2(*vSystsEfficiencyDenom[i],pot,kPOT,
424  costhetabins,eelecbins,kBinContent);
425  hHolder2->SetName(ana::UniqueName().c_str());
426  hHolder->Divide(hHolder2);
427  vSystsEff.push_back(hHolder);
428  }
429 
430  ////////////////////////////////////////////////////////////////////
431  ////////////////////// Flux ///////////////////////////////////////
432  ////////////////////////////////////////////////////////////////////
433  sprintf(name, "%s_%s_%s", "mc", "nominal", "flux");
434  Spectrum sflux_nom = *Spectrum::LoadFrom(fSysts, name);
435  TH1F* hFlux = (TH1F*)sflux_nom.ToTH1(pot);
436  hFlux->Scale(1e-4); //nu/cm^2 from nu/m^2
437 
438  std::vector<double> hAvgFluxVal;
439  hAvgFluxVal.push_back(hFlux->Integral());
440 
441  std::vector<std::unique_ptr<Spectrum>> vsSystsFlux;
442  std::vector<TH1F*> vSystsFlux;
443  for(uint i = 0; i < 100; i++){
444  sprintf(name, "%s_%s_%i_%s", "mc","ppfx",i,"flux");
445  vsSystsFlux.push_back(Spectrum::LoadFrom
446  (fSysts->GetDirectory(name)));
447  vSystsFlux.push_back((TH1F*)vsSystsFlux[i]->ToTH1(pot));
448  vSystsFlux[i]->Scale(1e-4); //nu/cm^2 from nu/m^2
449  hAvgFluxVal.push_back(vSystsFlux[i]->Integral());
450  }
451 
452  ////////////////////////////////////////////////////////////////////
453  ////////////////////// Unfolding////////////////////////////////////
454  ////////////////////////////////////////////////////////////////////
455  sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), "unfolding", "2d");
456  TH2F* hResponseMatrix =
457  (TH2F*)(Spectrum::LoadFrom(fAnalysis, name)))->ToTH2(pot;
458 
459  //Need to convert TH2 into TH1 for unfolding
460  int numBins = hSignalEst->GetXaxis()->GetNbins() *
461  hSignalEst->GetYaxis()->GetNbins();
462 
463  std::vector<TH1F*> hForUnfoldHists;
464  for(int j = 0; j <=(int)hSystsReco2D.size(); j++){
465  TH1F* hHolder = new TH1F(ana::UniqueName().c_str(),"",
466  numBins,0,numBins);
467  for(int i = 0; i <= hHolder->GetNbinsX();++i){
468  if(j == (int)hSystsReco2D.size()){
469  const int ix = i/ hSignalEst->GetYaxis()->GetNbins();
470  const int iy = i % hSignalEst->GetYaxis()->GetNbins();
471  const double val = hSignalEst->GetBinContent(ix+1, iy+1);
472  const double err = hSignalEst->GetBinError (ix+1, iy+1);
473  hHolder->SetBinContent(i+1,val);
474  hHolder->SetBinError(i+1,err);
475  }
476  else{
477  const int ix = i/ hSystsReco2D[j]->GetYaxis()->GetNbins();
478  const int iy = i % hSystsReco2D[j]->GetYaxis()->GetNbins();
479  const double val = hSystsReco2D[j]->GetBinContent(ix+1, iy+1);
480  const double err = hSystsReco2D[j]->GetBinError (ix+1, iy+1);
481  hHolder->SetBinContent(i+1,val);
482  hHolder->SetBinError(i+1,err);
483  }
484  }//loop over all kinematic bins
485  hForUnfoldHists.push_back(hHolder);
486  }//loop over all systematics + nominal
487  //Note last entry of the vector is the nominal histogram
488 
489  //TH1F* hUnfolded_SignalEst = DoUnfolding(hForUnfoldHists.back(),
490  // hResponseMatrix,8);
491  std::vector<TH1F*> hUnfoldedSysts;
492  for(int i = 0; i < (int)hForUnfoldHists.size(); i++){
493  // hUnfoldedSysts.push_back((TH1F*)DoUnfolding(hForUnfoldHists[i],
494  // hResponseMatrix,8));
495  hUnfoldedSysts.push_back((TH1F*)DoUnfolding(hForUnfoldHists[i],
496  hResponseMatrix,1)); // try with one iteration to see if we get input back
497  }
498 
499  //Now take unfolded hists and convert back to TH2s
500  std::vector<TH2F*> hUnfolded2D;
501  for(int j = 0; j < (int) hUnfoldedSysts.size(); j++){
502  TH2F* hHolder = (TH2F*)hSignalEst->Clone(ana::UniqueName().c_str());
503  for(int i = 0; i < hUnfoldedSysts[0]->GetNbinsX(); ++i){
504  const double val = hUnfoldedSysts[j]->GetBinContent(i+1);
505  const double err = hUnfoldedSysts[j]->GetBinError(i+1);
506  const int ix = i/ hSignalEst->GetYaxis()->GetNbins();
507  const int iy = i % hSignalEst->GetYaxis()->GetNbins();
508 
509  hHolder->SetBinContent(ix+1, iy+1, val);
510  hHolder->SetBinError (ix+1, iy+1, err);
511  }
512  hUnfolded2D.push_back(hHolder);
513  }
514 
515  TH2F* hUnfolded_SignalEst =
516  (TH2F*)hUnfolded2D.back()->Clone("hUnfoldedSignal");
517 
518 
519  ////////////////////////////////////////////////////////////////////
520  ////////////////////// Calculate XSec //////////////////////////////
521  ////////////////////////////////////////////////////////////////////
522  std::cout << "Calculate XSec" << std::endl;
523  TH2F* hXsec_cv_2d = (TH2F*)hUnfolded_SignalEst->Clone(ana::UniqueName().c_str());
524  hXsec_cv_2d->SetTitle("Cross Section");
525 
526  hXsec_cv_2d->Scale(1,"width");
527  hXsec_cv_2d->Divide(cv_eff);
528  hXsec_cv_2d->Scale(1./hAvgFluxVal[0]);
529  hXsec_cv_2d->Scale(1./nnucleons);
530  //hXsec_cv_2d->Scale(1e+38);
531 
532  std::cout << "Here " << std::endl;
533  std::vector<TH2F*> hXSec_systs;
534  for(int i = 0; i < (int) vSystsEff.size();i++){
535  TH2F* hHolder =
536  (TH2F*)hUnfolded_SignalEst->Clone(ana::UniqueName().c_str());
537 
538  hHolder->Scale(1,"width");
539  hHolder->Scale(1./nnucleons);
540  hHolder->Divide(vSystsEff[i]);
541 
542  if(i > 105){
543  hHolder->Scale(1./hAvgFluxVal[i-106]);
544  }
545  else{
546  hHolder->Scale(1./hAvgFluxVal[0]);
547  }
548  //hHolder->Scale(1e+38);
549  hXSec_systs.push_back(hHolder);
550  }
551 
552  //Get Systs Ready
553  //Order: {"lightdown","lightup","ckv", "calibneg","calibpos",
554  // "calibshape","genie","ppfx"};
555  std::cout << "Here 2" << std::endl;
556  hXSec_systs[0]->SetLineColor(kRed+3);
557  hXSec_systs[0]->SetLineStyle(7);
558  hXSec_systs[1]->SetLineColor(kRed-3);
559  hXSec_systs[1]->SetLineStyle(7);
560  hXSec_systs[2]->SetLineColor(kGreen);
561  hXSec_systs[3]->SetLineColor(kBlue-3);
562  hXSec_systs[4]->SetLineColor(kBlue+3);
563  hXSec_systs[5]->SetLineColor(kGreen);
564  hXSec_systs[5]->SetLineStyle(7);
565 
566  //Now For the multiverse stuff.
567  //Have to convert back to TH1s->Spectra->Multiverse->TH1->TH2
568  std::cout << "Here 3" << std::endl;
569  std::vector<TH1F*> hForGenieMultiverse;
570  for(int j = 6; j < 100+6; j++){
571  TH1F* hHolder = new TH1F(ana::UniqueName().c_str(),"",
572  numBins,0,numBins);
573  for(int i = 0; i <= hHolder->GetNbinsX();++i){
574  const int ix = i/ hXSec_systs[j]->GetYaxis()->GetNbins();
575  const int iy = i % hXSec_systs[j]->GetYaxis()->GetNbins();
576  const double val = hXSec_systs[j]->GetBinContent(ix+1, iy+1);
577  const double err = hXSec_systs[j]->GetBinError (ix+1, iy+1);
578  hHolder->SetBinContent(i+1,val);
579  hHolder->SetBinError(i+1,err);
580  }//loop over all kinematic bins
581  hForGenieMultiverse.push_back(hHolder);
582  }//loop over all genie multiverse xsecs
583 
584  std::cout << "Here 4" << std::endl;
585  std::vector<TH1F*> hForPPFXMultiverse;
586  for(int j = 106; j < 200+6; j++){
587  TH1F* hHolder = new TH1F(ana::UniqueName().c_str(),"",
588  numBins,0,numBins);
589  for(int i = 0; i <= hHolder->GetNbinsX();++i){
590  const int ix = i/ hXSec_systs[j]->GetYaxis()->GetNbins();
591  const int iy = i % hXSec_systs[j]->GetYaxis()->GetNbins();
592  const double val = hXSec_systs[j]->GetBinContent(ix+1, iy+1);
593  const double err = hXSec_systs[j]->GetBinError (ix+1, iy+1);
594  hHolder->SetBinContent(i+1,val);
595  hHolder->SetBinError(i+1,err);
596  }//loop over all kinematic bins
597  hForPPFXMultiverse.push_back(hHolder);
598  }//loop over all ppfx multiverse xsecs
599 
600  //Convert to spectrum
601  std::cout << "Here 5" << std::endl;
602  std::vector<Spectrum*> spects_xsec_genie;
603  std::vector<std::unique_ptr<Spectrum>> vspects_xsec_genie;
604  for(int i = 0; i < (int)hForGenieMultiverse.size(); i++){
605  Spectrum* spect = new Spectrum(hForGenieMultiverse[i],pot,0);
606  spects_xsec_genie.push_back(spect);
607  vspects_xsec_genie.push_back((std::unique_ptr<Spectrum>)
608  spects_xsec_genie[i]);
609  }
610 
611  std::cout << "Here 6" << std::endl;
612  std::vector<Spectrum*> spects_xsec_ppfx;
613  std::vector<std::unique_ptr<Spectrum>> vspects_xsec_ppfx;
614  for(int i = 0; i < (int)hForPPFXMultiverse.size(); i++){
615  Spectrum* spect = new Spectrum(hForPPFXMultiverse[i],pot,0);
616  spects_xsec_ppfx.push_back(spect);
617  vspects_xsec_ppfx.push_back((std::unique_ptr<Spectrum>)
618  spects_xsec_ppfx[i]);
619  }
620 
621  //Convert to Multiverse
622  GenieMultiverseSpectra genie_xsec_syst =
623  GenieMultiverseSpectra(100,vspects_xsec_genie,true);
624  GenieMultiverseSpectra ppfx_xsec_syst =
625  GenieMultiverseSpectra(100,vspects_xsec_ppfx,true);
626 
627  //Get +/- 1sigma curves (TH1)
628  TH1F* hGenie_up = (TH1F*)(genie_xsec_syst.UpperSigma())->ToTH1(pot);
629  TH1F* hGenie_dwn = (TH1F*)(genie_xsec_syst.LowerSigma())->ToTH1(pot);
630  TH1F* hPPFX_up = (TH1F*)(ppfx_xsec_syst.UpperSigma())->ToTH1(pot);
631  TH1F* hPPFX_dwn = (TH1F*)(ppfx_xsec_syst.LowerSigma())->ToTH1(pot);
632 
633  hGenie_up->SetLineColor(kOrange+3);
634  hGenie_dwn->SetLineColor(kOrange-3);
635  hPPFX_up->SetLineColor(kOrange+3);
636  hPPFX_dwn->SetLineColor(kOrange-3);
637  hPPFX_up->SetLineStyle(7);
638  hPPFX_dwn->SetLineStyle(7);
639 
640  //Convert Back to TH2s
641  TH2F* hGenie_up_2d =
642  (TH2F*)hSignalEst->Clone(ana::UniqueName().c_str());
643  TH2F* hGenie_down_2d =
644  (TH2F*)hSignalEst->Clone(ana::UniqueName().c_str());
645  TH2F* hPPFX_up_2d =
646  (TH2F*)hSignalEst->Clone(ana::UniqueName().c_str());
647  TH2F* hPPFX_down_2d =
648  (TH2F*)hSignalEst->Clone(ana::UniqueName().c_str());
649  for(int i = 0; i < hGenie_up->GetNbinsX(); ++i){
650  const double val_gu = hGenie_up->GetBinContent(i+1);
651  const double err_gu = hGenie_up->GetBinError(i+1);
652  const double val_gd = hGenie_dwn->GetBinContent(i+1);
653  const double err_gd = hGenie_dwn->GetBinError(i+1);
654  const double val_pu = hPPFX_up->GetBinContent(i+1);
655  const double err_pu = hPPFX_up->GetBinError(i+1);
656  const double val_pd = hPPFX_dwn->GetBinContent(i+1);
657  const double err_pd = hPPFX_dwn->GetBinError(i+1);
658  const int ix = i/ hSignalEst->GetYaxis()->GetNbins();
659  const int iy = i % hSignalEst->GetYaxis()->GetNbins();
660 
661  hGenie_up_2d->SetBinContent(ix+1, iy+1, val_gu);
662  hGenie_up_2d->SetBinError (ix+1, iy+1, err_gu);
663  hGenie_down_2d->SetBinContent(ix+1, iy+1, val_gd);
664  hGenie_down_2d->SetBinError (ix+1, iy+1, err_gd);
665  hPPFX_up_2d->SetBinContent(ix+1, iy+1, val_pu);
666  hPPFX_up_2d->SetBinError (ix+1, iy+1, err_pu);
667  hPPFX_down_2d->SetBinContent(ix+1, iy+1, val_pd);
668  hPPFX_down_2d->SetBinError (ix+1, iy+1, err_pd);
669  }
670 
671  std::vector<TH2F*>systs_up_2d =
672  {(TH2F*)hGenie_up_2d->Clone(ana::UniqueName().c_str()),
673  (TH2F*)hPPFX_up_2d->Clone(ana::UniqueName().c_str()),
674  (TH2F*)hXSec_systs[1]->Clone(ana::UniqueName().c_str()),
675  (TH2F*)hXSec_systs[4]->Clone(ana::UniqueName().c_str()),
676  (TH2F*)hXSec_systs[2]->Clone(ana::UniqueName().c_str()),
677  (TH2F*)hXSec_systs[5]->Clone(ana::UniqueName().c_str())
678  };
679 
680 
681  std::vector<TH2F*> systs_down_2d =
682  {(TH2F*)hGenie_down_2d->Clone(ana::UniqueName().c_str()),
683  (TH2F*)hPPFX_down_2d->Clone(ana::UniqueName().c_str()),
684  (TH2F*)hXSec_systs[0]->Clone(ana::UniqueName().c_str()),
685  (TH2F*)hXSec_systs[3]->Clone(ana::UniqueName().c_str()),
686  (TH2F*)hXsec_cv_2d->Clone(ana::UniqueName().c_str()),
687  (TH2F*)hXsec_cv_2d->Clone(ana::UniqueName().c_str())};
688 
689  TH2F* hCVHolder_2d = (TH2F*)hXsec_cv_2d->Clone(ana::UniqueName().c_str());
690 
691  //Get True Cross Section From Fake Data
692  sprintf(name, "%s_%s_%s_%s","mc",dataName.c_str(),"efficiency", "denom_2d");
693  Spectrum sTrueSignal =
694  *Spectrum::LoadFrom(fAnalysis, name);
695  TH2F* hTrueSignal_2d = (TH2F*)ana::ToTH2(sTrueSignal,pot, kPOT,
696  costhetabins,eelecbins,kBinContent);
697  hTrueSignal_2d->Scale(1./hAvgFluxVal[0]);
698  hTrueSignal_2d->Scale(1./nnucleons);
699  hTrueSignal_2d->Scale(1, "width");
700  //hTrueSignal_2d->Scale(1e38);
701  hTrueSignal_2d->SetLineColor(kRed);
702 
703 
704  //Finally Project Results down to a single axis
705  for(int xbin = 0; xbin <= hXsec_cv_2d->GetXaxis()->GetNbins();xbin++){
706  std::stringstream low_X;
707  low_X << std::fixed << std::setprecision(2) <<
708  hXsec_cv_2d->GetXaxis()->GetBinLowEdge(xbin);
709  std::stringstream hi_X;
710  hi_X << std::fixed << std::setprecision(2) <<
711  hXsec_cv_2d->GetXaxis()->GetBinUpEdge(xbin);
712 
713  std::string caption = low_X.str() + " #leq cos #theta_{e} < "+
714  hi_X.str();
715 
716  TH1F* hXsec_cv =
717  (TH1F*)hXsec_cv_2d->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
718  hXsec_cv->GetXaxis()->SetTitle("Electron Energy, E_{e} (GeV)");
719  TH1F* hTrueSignal =
720  (TH1F*)hTrueSignal_2d->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
721  TH1F* hCVHolder =
722  (TH1F*)hXsec_cv_2d->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
723 
724  std::vector<TH1F*>systs_up;
725  std::vector<TH1F*>systs_down;
726 
727  for(int i = 0; i < (int)systs_up_2d.size(); i++){
728  systs_up.push_back((TH1F*)systs_up_2d[i]->
729  ProjectionY(ana::UniqueName().c_str(),xbin,xbin));
730  systs_down.push_back((TH1F*)systs_down_2d[i]->
731  ProjectionY(ana::UniqueName().c_str(),xbin,xbin));
732  }
733 
734  TH1F* hUnfolded_SignalEst_1d =
735  (TH1F*)hUnfolded_SignalEst->ProjectionY(ana::UniqueName().c_str(),
736  xbin,xbin);
737 
738  TGraphAsymmErrors* g =
739  PlotSystErrorBand(hCVHolder,
740  systs_up,systs_down,
741  -1,-1,1.3,false);
742 
743 
744  //Propagate Uncertainties
745  //CV
746  hXsec_cv->SetMarkerColor(kBlack);
747  hXsec_cv->SetLineColor(kBlack);
748  TH1F* hUncert_tot = (TH1F*)hXsec_cv->Clone(ana::UniqueName().c_str());
749  for(int i = 1; i <= (int)hUncert_tot->GetXaxis()->GetNbins();i++){
750  double uncert = hUnfolded_SignalEst_1d->GetBinError(i);
751  double total = hUncert_tot->GetBinContent(i);
752  double selected_tot = hUnfolded_SignalEst_1d->GetBinContent(i);
753  double err_hi = g->GetErrorYhigh(i);
754  double err_lo = g->GetErrorYlow(i);
755 
756  //double err_fluxeff = (err_hi)/total;
757  double err_fluxeff = (err_hi > err_lo)? err_hi/total : err_lo/total; // test relative uncertainties
758 
759 
760  double tot_err = total*sqrt(pow(uncert/selected_tot,2) +
761  pow(err_fluxeff,2));
762  if(total == 0) tot_err = 0;
763 
764  hUncert_tot->SetBinError(i,tot_err);
765 
766  std::cout << i << "\t" << total << "\t" <<
767  uncert << "\t" << selected_tot << "\t" <<
768  err_fluxeff << "\t" <<
769  sqrt(pow(uncert/selected_tot,2) +
770  pow(err_fluxeff,2)) << std::endl;
771 
772  }
773 
774  hUncert_tot->GetYaxis()->SetTitle("#frac{d^{2} #sigma}{dcos#theta dE_{e}} #left(#frac{10^{-38} cm^{2}}{nucleon GeV}#right)");
775  hUncert_tot->SetTitle(caption.c_str());
776 
777  std::vector<TH2F*>systs_up_2d =
778  {(TH2F*)hGenie_up_2d->Clone(ana::UniqueName().c_str()),
779  (TH2F*)hPPFX_up_2d->Clone(ana::UniqueName().c_str()),
780  (TH2F*)hXSec_systs[1]->Clone(ana::UniqueName().c_str()),
781  (TH2F*)hXSec_systs[4]->Clone(ana::UniqueName().c_str()),
782  (TH2F*)hXSec_systs[2]->Clone(ana::UniqueName().c_str()),
783  (TH2F*)hXSec_systs[5]->Clone(ana::UniqueName().c_str())
784  };
785 
786 
787  std::vector<TH2F*> systs_down_2d =
788  {(TH2F*)hGenie_down_2d->Clone(ana::UniqueName().c_str()),
789  (TH2F*)hPPFX_down_2d->Clone(ana::UniqueName().c_str()),
790  (TH2F*)hXSec_systs[0]->Clone(ana::UniqueName().c_str()),
791  (TH2F*)hXSec_systs[3]->Clone(ana::UniqueName().c_str()),
792  (TH2F*)hXsec_cv_2d->Clone(ana::UniqueName().c_str()),
793  (TH2F*)hXsec_cv_2d->Clone(ana::UniqueName().c_str())};
794 
795 
796 
797  TCanvas* cResult = new TCanvas(ana::UniqueName().c_str(), "cResult");
798  cResult->SetLeftMargin(0.15);
799  hUncert_tot->GetXaxis()->SetRangeUser(1.0,10.0);
800  hUncert_tot->SetMarkerStyle(20);
801  hUncert_tot->Draw();
802  g->Draw("e2 same");
803  systs_down[2]->Draw("hist same");
804  systs_up[2]->Draw("hist same");
805  systs_up[4]->Draw("hist same");
806  systs_down[3]->Draw("hist same");
807  systs_up[3]->Draw("hist same");
808  systs_up[5]->Draw("hist same");
809  systs_up[0]->Draw("hist same");
810  systs_down[0]->Draw("hist same");
811  systs_up[1]->Draw("hist same");
812  systs_down[1]->Draw("hist same");
813  hTrueSignal->Draw("hist same");
814  hUncert_tot->Draw("same");
815  TLegend* leg_total = ana::AutoPlaceLegend(0.3,0.3);
816  leg_total->AddEntry(hTrueSignal, "Truth", "l");
817  leg_total->AddEntry(hUncert_tot, "FakeData", "p");
818  leg_total->AddEntry(systs_down[3], "Calibration -1 #sigma", "l");
819  leg_total->AddEntry(systs_up[3], "Calibration +1 #sigma", "l");
820  leg_total->AddEntry(systs_down[2], "Light Level -1 #sigma", "l");
821  leg_total->AddEntry(systs_up[2], "Light Level +1 #sigma", "l");
822  leg_total->AddEntry(systs_up[4], "Cherenkov Variation", "l");
823  leg_total->AddEntry(systs_up[5], "Calibration Shape", "l");
824  leg_total->AddEntry(systs_down[1], "Flux -1 #sigma", "l");
825  leg_total->AddEntry(systs_up[1], "Flux +1 #sigma", "l");
826  leg_total->AddEntry(systs_down[0], "#nu-Interaction -1 #sigma", "l");
827  leg_total->AddEntry(systs_up[0], "#nu-Interaction +1 #sigma", "l");
828  leg_total->Draw();
829  sprintf(name, "%s%s_%s_%s_%s_%i.png","AnalysisPlots/", dataName.c_str(),
830  "xsec", "results_all_systs",xsecResult.c_str(),xbin);
831  cResult->SaveAs(name);
832 
833  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str(), "cResult");
834  c3->SetLeftMargin(0.15);
835  hUncert_tot->GetXaxis()->SetRangeUser(1.0,10.0);
836  hUncert_tot->SetMarkerStyle(20);
837  hUncert_tot->Draw();
838  hTrueSignal->Draw("hist same");
839  hUncert_tot->Draw("same");
840  sprintf(name, "%s%s_%s_%s_%s_%i.png","AnalysisPlots/", dataName.c_str(),
841  "xsec", "results_truth_compare",xsecResult.c_str(),xbin);
842  c3->SaveAs(name);
843 
844  //Calculate Relative Uncertainties
845  //CV
846  TH1F* hRelUnc_tot = (TH1F*)hUncert_tot->Clone(ana::UniqueName().c_str());
847  for(int i = 0; i <= (int)hRelUnc_tot->GetXaxis()->GetNbins();i++){
848  double uncert = hRelUnc_tot->GetBinError(i);
849  double total = hRelUnc_tot->GetBinContent(i);
850  if(total == 0) hRelUnc_tot->SetBinContent(i,0);
851  else hRelUnc_tot->SetBinContent(i,uncert/total);
852  }
853 
854  std::cout << "Here" << std::endl;
855  std::vector<TH1F*> hRelUncert;
856  for(int i = 0; i < (int)systs_up.size();i++){
857  TH1F* hHolder = (TH1F*)systs_up[i]->Clone(ana::UniqueName().c_str());
858  TH1F* hHolderLo =
859  (TH1F*)systs_down[i]->Clone(ana::UniqueName().c_str());
860 
861  for(int ibin = 0; ibin <= hHolder->GetXaxis()->GetNbins(); ibin++){
862  double val_up = hHolder->GetBinContent(ibin);
863  double val_down = hHolderLo->GetBinContent(ibin);
864  double val_cv = hXsec_cv->GetBinContent(ibin);
865 
866  //double uncert = 0.5*( fabs(val_up - val_cv) +
867  // fabs(val_down - val_cv));
868  //double uncert = fabs(val_up - val_cv); // test rel uncertainties
869  double uncert = std::max(fabs(val_up - val_cv), fabs(val_down - val_cv));
870 
871  if(val_cv == 0) hHolder->SetBinContent(ibin,0);
872  else hHolder->SetBinContent(ibin,uncert/val_cv);
873  }
874  hRelUncert.push_back(hHolder);
875  }
876  std::cout << "Here 2" << std::endl;
877 
878  TH1F* hRelUnc_fit = (TH1F*)hXsec_cv->Clone(ana::UniqueName().c_str());
879  for(int i = 0; i <= (int)hRelUnc_tot->GetXaxis()->GetNbins();i++){
880  double uncert = hRelUnc_fit->GetBinError(i);
881  double total = hRelUnc_fit->GetBinContent(i);
882  if(total == 0) hRelUnc_fit->SetBinContent(i,0);
883  else hRelUnc_fit->SetBinContent(i,uncert/total);
884  }
885 
886  std::cout << "Here 3" << std::endl;
887 
888  //Same Order as vector used to generate TGraphAsmyErrors
889  hRelUncert[5]->SetLineColor(kGreen+2); //calibshape
890  hRelUncert[4]->SetLineColor(kOrange-3); //ckv
891  hRelUncert[0]->SetLineColor(kRed - 2); //genie
892  hRelUncert[1]->SetLineColor(kCyan-3); //ppfx
893  hRelUncert[2]->SetLineColor(kBlue-4); //light
894  hRelUncert[3]->SetLineColor(kMagenta+1); //calib
895  hRelUnc_fit->SetLineColor(kBlack);
896  hRelUnc_fit->SetLineStyle(7);
897  hRelUnc_tot->SetLineColor(kBlack);
898 
899  TCanvas* c2 = new TCanvas(ana::UniqueName().c_str(),"c2");
900  hRelUnc_tot->GetYaxis()->SetTitle("Relative Uncertainty");
901  hRelUnc_tot->SetTitle(" ");
902  hRelUnc_tot->GetYaxis()->SetRangeUser(0,1.0);
903  hRelUnc_tot->GetXaxis()->SetRangeUser(1,10);
904  hRelUnc_tot->Draw("hist");
905  hRelUncert[5]->Draw("hist same");
906  hRelUncert[4]->Draw("hist same");
907  hRelUncert[3]->Draw("hist same");
908  hRelUncert[2]->Draw("hist same");
909  hRelUncert[1]->Draw("hist same");
910  hRelUncert[0]->Draw("hist same");
911  hRelUnc_fit->Draw("hist same");
912  hRelUnc_tot->Draw("hist same");
913  TLegend* leg_systs = ana::AutoPlaceLegend(0.3,0.3);
914  leg_systs->AddEntry(hRelUnc_tot, "Total", "l");
915  leg_systs->AddEntry(hRelUnc_fit, "Fit", "l");
916  leg_systs->AddEntry(hRelUncert[3], "Calibration Normalization", "l");
917  leg_systs->AddEntry(hRelUncert[2], "Light Level", "l");
918  leg_systs->AddEntry(hRelUncert[1], "Flux", "l");
919  leg_systs->AddEntry(hRelUncert[0], "#nu-Interaction", "l");
920  leg_systs->AddEntry(hRelUncert[4], "Cherenkov Variation", "l");
921  leg_systs->AddEntry(hRelUncert[5], "Calibration Shape", "l");
922  leg_systs->Draw();
923  sprintf(name, "%s%s_%s_%s_%s_%i.png","AnalysisPlots/", dataName.c_str(),
924  "xsec", "relative_uncertantiy",xsecResult.c_str(),xbin);
925  c2->SaveAs(name);
926  }//loop over xbins
927 
928  ////////////////////////////////////////////////////////////////////////
929  ////////////////Single Differentail Results/////////////////////////////
930  ////////////////////////////////////////////////////////////////////////
931  {
932  TH1F* hXsec_cv =
933  (TH1F*)hXsec_cv_2d->ProjectionY(ana::UniqueName().c_str());
934  TH1F* hTrueSignal =
935  (TH1F*)hTrueSignal_2d->ProjectionY(ana::UniqueName().c_str());
936  TH1F* hCVHolder =
937  (TH1F*)hXsec_cv_2d->ProjectionY(ana::UniqueName().c_str());
938 
939  std::vector<TH1F*>systs_up;
940  std::vector<TH1F*>systs_down;
941 
942  for(int i = 0; i < (int)systs_up_2d.size(); i++){
943  systs_up.push_back((TH1F*)systs_up_2d[i]->
944  ProjectionY(ana::UniqueName().c_str()));
945  systs_down.push_back((TH1F*)systs_down_2d[i]->
946  ProjectionY(ana::UniqueName().c_str()));
947  }
948 
949  TH1F* hUnfolded_SignalEst_1d =
950  (TH1F*)hUnfolded_SignalEst->ProjectionY(ana::UniqueName().c_str());
951 
952  TGraphAsymmErrors* g =
953  PlotSystErrorBand(hCVHolder,
954  systs_up,systs_down,
955  -1,-1,1.3,false);
956 
957 
958  //Propagate Uncertainties
959  //CV
960  hXsec_cv->SetMarkerColor(kBlack);
961  hXsec_cv->SetLineColor(kBlack);
962  TH1F* hUncert_tot = (TH1F*)hXsec_cv->Clone(ana::UniqueName().c_str());
963  for(int i = 1; i <= (int)hUncert_tot->GetXaxis()->GetNbins();i++){
964  double uncert = hUnfolded_SignalEst_1d->GetBinError(i);
965  double total = hUncert_tot->GetBinContent(i);
966  double selected_tot = hUnfolded_SignalEst_1d->GetBinContent(i);
967  double err_hi = g->GetErrorYhigh(i);
968  double err_lo = g->GetErrorYlow(i);
969 
970  //double err_fluxeff = (err_hi)/total;
971  double err_fluxeff = (err_hi > err_lo)? err_hi/total : err_lo/total;
972 
973 
974  double tot_err = total*sqrt(pow(uncert/selected_tot,2) +
975  pow(err_fluxeff,2));
976  if(total == 0) tot_err = 0;
977 
978  hUncert_tot->SetBinError(i,tot_err);
979 
980  std::cout << i << "\t" << total << "\t" <<
981  uncert << "\t" << selected_tot << "\t" <<
982  err_fluxeff << "\t" <<
983  sqrt(pow(uncert/selected_tot,2) +
984  pow(err_fluxeff,2)) << std::endl;
985 
986  }
987 
988  hUncert_tot->GetYaxis()->SetTitle("#frac{d #sigma}{dE_{e}} #left(#frac{10^{-38} cm^{2}}{nucleon GeV}#right)");
989  hUncert_tot->GetXaxis()->SetTitle("Electron Energy, E_{#e} (GeV)");
990  hUncert_tot->SetTitle("");
991 
992  TCanvas* cResult = new TCanvas(ana::UniqueName().c_str(), "cResult");
993  cResult->SetLeftMargin(0.15);
994  hUncert_tot->GetXaxis()->SetRangeUser(1.0,10.0);
995  hUncert_tot->SetMarkerStyle(20);
996  hUncert_tot->Draw();
997  g->Draw("e2 same");
998  systs_down[2]->Draw("hist same");
999  systs_up[2]->Draw("hist same");
1000  systs_up[4]->Draw("hist same");
1001  systs_down[3]->Draw("hist same");
1002  systs_up[3]->Draw("hist same");
1003  systs_up[5]->Draw("hist same");
1004  systs_up[0]->Draw("hist same");
1005  systs_down[0]->Draw("hist same");
1006  systs_up[1]->Draw("hist same");
1007  systs_down[1]->Draw("hist same");
1008  hTrueSignal->Draw("hist same");
1009  hUncert_tot->Draw("same");
1010  TLegend* leg_total = ana::AutoPlaceLegend(0.3,0.3);
1011  leg_total->AddEntry(hTrueSignal, "Truth", "l");
1012  leg_total->AddEntry(hUncert_tot, "FakeData", "p");
1013  leg_total->AddEntry(systs_down[3], "Calibration -1 #sigma", "l");
1014  leg_total->AddEntry(systs_up[3], "Calibration +1 #sigma", "l");
1015  leg_total->AddEntry(systs_down[2], "Light Level -1 #sigma", "l");
1016  leg_total->AddEntry(systs_up[2], "Light Level +1 #sigma", "l");
1017  leg_total->AddEntry(systs_up[4], "Cherenkov Variation", "l");
1018  leg_total->AddEntry(systs_up[5], "Calibration Shape", "l");
1019  leg_total->AddEntry(systs_down[1], "Flux -1 #sigma", "l");
1020  leg_total->AddEntry(systs_up[1], "Flux +1 #sigma", "l");
1021  leg_total->AddEntry(systs_down[0], "#nu-Interaction -1 #sigma", "l");
1022  leg_total->AddEntry(systs_up[0], "#nu-Interaction +1 #sigma", "l");
1023  leg_total->Draw();
1024  sprintf(name, "%s%s_%s_%s_%s_%s.png","AnalysisPlots/", dataName.c_str(),
1025  "xsec", "results_all_systs",xsecResult.c_str(),"elecE");
1026  cResult->SaveAs(name);
1027 
1028  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str(), "cResult");
1029  c3->SetLeftMargin(0.15);
1030  hUncert_tot->GetXaxis()->SetRangeUser(1.0,10.0);
1031  hUncert_tot->SetMarkerStyle(20);
1032  hUncert_tot->Draw();
1033  hTrueSignal->Draw("hist same");
1034  hUncert_tot->Draw("same");
1035  sprintf(name, "%s%s_%s_%s_%s_%s.png","AnalysisPlots/", dataName.c_str(),
1036  "xsec", "results_truth_compare",xsecResult.c_str(),"elecE");
1037  c3->SaveAs(name);
1038 
1039  //Calculate Relative Uncertainties
1040  //CV
1041  TH1F* hRelUnc_tot = (TH1F*)hUncert_tot->Clone(ana::UniqueName().c_str());
1042  for(int i = 0; i <= (int)hRelUnc_tot->GetXaxis()->GetNbins();i++){
1043  double uncert = hRelUnc_tot->GetBinError(i);
1044  double total = hRelUnc_tot->GetBinContent(i);
1045  if(total == 0) hRelUnc_tot->SetBinContent(i,0);
1046  else hRelUnc_tot->SetBinContent(i,uncert/total);
1047  }
1048 
1049  std::cout << "Here" << std::endl;
1050  std::vector<TH1F*> hRelUncert;
1051  for(int i = 0; i < (int)systs_up.size();i++){
1052  TH1F* hHolder = (TH1F*)systs_up[i]->Clone(ana::UniqueName().c_str());
1053  TH1F* hHolderLo =
1054  (TH1F*)systs_down[i]->Clone(ana::UniqueName().c_str());
1055 
1056  for(int ibin = 0; ibin <= hHolder->GetXaxis()->GetNbins(); ibin++){
1057  double val_up = hHolder->GetBinContent(ibin);
1058  double val_down = hHolderLo->GetBinContent(ibin);
1059  double val_cv = hXsec_cv->GetBinContent(ibin);
1060 
1061  // double uncert = 0.5*( fabs(val_up - val_cv) +
1062  //fabs(val_down - val_cv));
1063  double uncert = std::max(fabs(val_up - val_cv), fabs(val_down - val_cv));
1064  if(val_cv == 0) hHolder->SetBinContent(ibin,0);
1065  else hHolder->SetBinContent(ibin,uncert/val_cv);
1066  }
1067  hRelUncert.push_back(hHolder);
1068  }
1069  std::cout << "Here 2" << std::endl;
1070 
1071  TH1F* hRelUnc_fit = (TH1F*)hXsec_cv->Clone(ana::UniqueName().c_str());
1072  for(int i = 0; i <= (int)hRelUnc_tot->GetXaxis()->GetNbins();i++){
1073  double uncert = hRelUnc_fit->GetBinError(i);
1074  double total = hRelUnc_fit->GetBinContent(i);
1075  if(total == 0) hRelUnc_fit->SetBinContent(i,0);
1076  else hRelUnc_fit->SetBinContent(i,uncert/total);
1077  }
1078 
1079  std::cout << "Here 3" << std::endl;
1080 
1081  //Same Order as vector used to generate TGraphAsmyErrors
1082  hRelUncert[5]->SetLineColor(kGreen+2); //calibshape
1083  hRelUncert[4]->SetLineColor(kOrange-3); //ckv
1084  hRelUncert[0]->SetLineColor(kRed - 2); //genie
1085  hRelUncert[1]->SetLineColor(kCyan-3); //ppfx
1086  hRelUncert[2]->SetLineColor(kBlue-4); //light
1087  hRelUncert[3]->SetLineColor(kMagenta+1); //calib
1088  hRelUnc_fit->SetLineColor(kBlack);
1089  hRelUnc_fit->SetLineStyle(7);
1090  hRelUnc_tot->SetLineColor(kBlack);
1091 
1092  TCanvas* c2 = new TCanvas(ana::UniqueName().c_str(),"c2");
1093  hRelUnc_tot->GetYaxis()->SetTitle("Relative Uncertainty");
1094  hRelUnc_tot->SetTitle(" ");
1095  hRelUnc_tot->GetYaxis()->SetRangeUser(0,1.0);
1096  hRelUnc_tot->GetXaxis()->SetRangeUser(1,10);
1097  hRelUnc_tot->Draw("hist");
1098  hRelUncert[5]->Draw("hist same");
1099  hRelUncert[4]->Draw("hist same");
1100  hRelUncert[3]->Draw("hist same");
1101  hRelUncert[2]->Draw("hist same");
1102  hRelUncert[1]->Draw("hist same");
1103  hRelUncert[0]->Draw("hist same");
1104  hRelUnc_fit->Draw("hist same");
1105  hRelUnc_tot->Draw("hist same");
1106  TLegend* leg_systs = ana::AutoPlaceLegend(0.3,0.3);
1107  leg_systs->AddEntry(hRelUnc_tot, "Total", "l");
1108  leg_systs->AddEntry(hRelUnc_fit, "Fit", "l");
1109  leg_systs->AddEntry(hRelUncert[3], "Calibration Normalization", "l");
1110  leg_systs->AddEntry(hRelUncert[2], "Light Level", "l");
1111  leg_systs->AddEntry(hRelUncert[1], "Flux", "l");
1112  leg_systs->AddEntry(hRelUncert[0], "#nu-Interaction", "l");
1113  leg_systs->AddEntry(hRelUncert[4], "Cherenkov Variation", "l");
1114  leg_systs->AddEntry(hRelUncert[5], "Calibration Shape", "l");
1115  leg_systs->Draw();
1116  sprintf(name, "%s%s_%s_%s_%s_%s.png","AnalysisPlots/", dataName.c_str(),
1117  "xsec", "relative_uncertantiy",xsecResult.c_str(),"elecE");
1118  c2->SaveAs(name);
1119  }
1120 
1121  {
1122  TH1F* hXsec_cv =
1123  (TH1F*)hXsec_cv_2d->ProjectionX(ana::UniqueName().c_str());
1124  TH1F* hTrueSignal =
1125  (TH1F*)hTrueSignal_2d->ProjectionX(ana::UniqueName().c_str());
1126  TH1F* hCVHolder =
1127  (TH1F*)hXsec_cv_2d->ProjectionX(ana::UniqueName().c_str());
1128 
1129  std::vector<TH1F*>systs_up;
1130  std::vector<TH1F*>systs_down;
1131 
1132  for(int i = 0; i < (int)systs_up_2d.size(); i++){
1133  systs_up.push_back((TH1F*)systs_up_2d[i]->
1134  ProjectionX(ana::UniqueName().c_str()));
1135  systs_down.push_back((TH1F*)systs_down_2d[i]->
1136  ProjectionX(ana::UniqueName().c_str()));
1137  }
1138 
1139  TH1F* hUnfolded_SignalEst_1d =
1140  (TH1F*)hUnfolded_SignalEst->ProjectionX(ana::UniqueName().c_str());
1141 
1142  TGraphAsymmErrors* g =
1143  PlotSystErrorBand(hCVHolder,
1144  systs_up,systs_down,
1145  -1,-1,1.3,false);
1146 
1147 
1148  //Propagate Uncertainties
1149  //CV
1150  hXsec_cv->SetMarkerColor(kBlack);
1151  hXsec_cv->SetLineColor(kBlack);
1152  TH1F* hUncert_tot = (TH1F*)hXsec_cv->Clone(ana::UniqueName().c_str());
1153  for(int i = 1; i <= (int)hUncert_tot->GetXaxis()->GetNbins();i++){
1154  double uncert = hUnfolded_SignalEst_1d->GetBinError(i);
1155  double total = hUncert_tot->GetBinContent(i);
1156  double selected_tot = hUnfolded_SignalEst_1d->GetBinContent(i);
1157  double err_hi = g->GetErrorYhigh(i);
1158  double err_lo = g->GetErrorYlow(i);
1159 
1160  //double err_fluxeff = (err_hi)/total;
1161  double err_fluxeff = (err_hi > err_lo) ? err_hi/total : err_lo/total;
1162 
1163 
1164  double tot_err = total*sqrt(pow(uncert/selected_tot,2) +
1165  pow(err_fluxeff,2));
1166  if(total == 0) tot_err = 0;
1167 
1168  hUncert_tot->SetBinError(i,tot_err);
1169 
1170  std::cout << i << "\t" << total << "\t" <<
1171  uncert << "\t" << selected_tot << "\t" <<
1172  err_fluxeff << "\t" <<
1173  sqrt(pow(uncert/selected_tot,2) +
1174  pow(err_fluxeff,2)) << std::endl;
1175 
1176  }
1177 
1178  hUncert_tot->GetYaxis()->SetTitle("#frac{d #sigma}{dcos #theta} #left(#frac{10^{-38} cm^{2}}{nucleon}#right)");
1179  hUncert_tot->GetXaxis()->SetTitle("cos #theta_{e}");
1180  hUncert_tot->SetTitle("");
1181 
1182  TCanvas* cResult = new TCanvas(ana::UniqueName().c_str(), "cResult");
1183  cResult->SetLeftMargin(0.15);
1184  hUncert_tot->GetXaxis()->SetRangeUser(0.75,1.0);
1185  hUncert_tot->SetMarkerStyle(20);
1186  hUncert_tot->Draw();
1187  g->Draw("e2 same");
1188  systs_down[2]->Draw("hist same");
1189  systs_up[2]->Draw("hist same");
1190  systs_up[4]->Draw("hist same");
1191  systs_down[3]->Draw("hist same");
1192  systs_up[3]->Draw("hist same");
1193  systs_up[5]->Draw("hist same");
1194  systs_up[0]->Draw("hist same");
1195  systs_down[0]->Draw("hist same");
1196  systs_up[1]->Draw("hist same");
1197  systs_down[1]->Draw("hist same");
1198  hTrueSignal->Draw("hist same");
1199  hUncert_tot->Draw("same");
1200  TLegend* leg_total = ana::AutoPlaceLegend(0.3,0.3);
1201  leg_total->AddEntry(hTrueSignal, "Truth", "l");
1202  leg_total->AddEntry(hUncert_tot, "FakeData", "p");
1203  leg_total->AddEntry(systs_down[3], "Calibration -1 #sigma", "l");
1204  leg_total->AddEntry(systs_up[3], "Calibration +1 #sigma", "l");
1205  leg_total->AddEntry(systs_down[2], "Light Level -1 #sigma", "l");
1206  leg_total->AddEntry(systs_up[2], "Light Level +1 #sigma", "l");
1207  leg_total->AddEntry(systs_up[4], "Cherenkov Variation", "l");
1208  leg_total->AddEntry(systs_up[5], "Calibration Shape", "l");
1209  leg_total->AddEntry(systs_down[1], "Flux -1 #sigma", "l");
1210  leg_total->AddEntry(systs_up[1], "Flux +1 #sigma", "l");
1211  leg_total->AddEntry(systs_down[0], "#nu-Interaction -1 #sigma", "l");
1212  leg_total->AddEntry(systs_up[0], "#nu-Interaction +1 #sigma", "l");
1213  leg_total->Draw();
1214  sprintf(name, "%s%s_%s_%s_%s_%s.png","AnalysisPlots/", dataName.c_str(),
1215  "xsec", "results_all_systs",xsecResult.c_str(),"costheta");
1216  cResult->SaveAs(name);
1217 
1218  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str(), "cResult");
1219  c3->SetLeftMargin(0.15);
1220  hUncert_tot->GetXaxis()->SetRangeUser(0.75,1.0);
1221  hUncert_tot->SetMarkerStyle(20);
1222  hUncert_tot->Draw();
1223  hTrueSignal->Draw("hist same");
1224  hUncert_tot->Draw("same");
1225  sprintf(name, "%s%s_%s_%s_%s_%s.png","AnalysisPlots/", dataName.c_str(),
1226  "xsec", "results_truth_compare",xsecResult.c_str(),"costheta");
1227  c3->SaveAs(name);
1228 
1229  //Calculate Relative Uncertainties
1230  //CV
1231  TH1F* hRelUnc_tot = (TH1F*)hUncert_tot->Clone(ana::UniqueName().c_str());
1232  for(int i = 0; i <= (int)hRelUnc_tot->GetXaxis()->GetNbins();i++){
1233  double uncert = hRelUnc_tot->GetBinError(i);
1234  double total = hRelUnc_tot->GetBinContent(i);
1235  if(total == 0) hRelUnc_tot->SetBinContent(i,0);
1236  else hRelUnc_tot->SetBinContent(i,uncert/total);
1237  }
1238 
1239  std::cout << "Here" << std::endl;
1240  std::vector<TH1F*> hRelUncert;
1241  for(int i = 0; i < (int)systs_up.size();i++){
1242  TH1F* hHolder = (TH1F*)systs_up[i]->Clone(ana::UniqueName().c_str());
1243  TH1F* hHolderLo =
1244  (TH1F*)systs_down[i]->Clone(ana::UniqueName().c_str());
1245 
1246  for(int ibin = 0; ibin <= hHolder->GetXaxis()->GetNbins(); ibin++){
1247  double val_up = hHolder->GetBinContent(ibin);
1248  double val_down = hHolderLo->GetBinContent(ibin);
1249  double val_cv = hXsec_cv->GetBinContent(ibin);
1250 
1251  // double uncert = 0.5*( fabs(val_up - val_cv) +
1252  // fabs(val_down - val_cv));
1253  double uncert = std::max(fabs(val_up - val_cv), fabs(val_down - val_cv));
1254  if(val_cv == 0) hHolder->SetBinContent(ibin,0);
1255  else hHolder->SetBinContent(ibin,uncert/val_cv);
1256  }
1257  hRelUncert.push_back(hHolder);
1258  }
1259  std::cout << "Here 2" << std::endl;
1260 
1261  TH1F* hRelUnc_fit = (TH1F*)hXsec_cv->Clone(ana::UniqueName().c_str());
1262  for(int i = 0; i <= (int)hRelUnc_tot->GetXaxis()->GetNbins();i++){
1263  double uncert = hRelUnc_fit->GetBinError(i);
1264  double total = hRelUnc_fit->GetBinContent(i);
1265  if(total == 0) hRelUnc_fit->SetBinContent(i,0);
1266  else hRelUnc_fit->SetBinContent(i,uncert/total);
1267  }
1268 
1269  std::cout << "Here 3" << std::endl;
1270 
1271  //Same Order as vector used to generate TGraphAsmyErrors
1272  hRelUncert[5]->SetLineColor(kGreen+2); //calibshape
1273  hRelUncert[4]->SetLineColor(kOrange-3); //ckv
1274  hRelUncert[0]->SetLineColor(kRed - 2); //genie
1275  hRelUncert[1]->SetLineColor(kCyan-3); //ppfx
1276  hRelUncert[2]->SetLineColor(kBlue-4); //light
1277  hRelUncert[3]->SetLineColor(kMagenta+1); //calib
1278  hRelUnc_fit->SetLineColor(kBlack);
1279  hRelUnc_fit->SetLineStyle(7);
1280  hRelUnc_tot->SetLineColor(kBlack);
1281 
1282  TCanvas* c2 = new TCanvas(ana::UniqueName().c_str(),"c2");
1283  hRelUnc_tot->GetYaxis()->SetTitle("Relative Uncertainty");
1284  hRelUnc_tot->SetTitle(" ");
1285  hRelUnc_tot->GetYaxis()->SetRangeUser(0,1.0);
1286  hRelUnc_tot->GetXaxis()->SetRangeUser(0.75,1);
1287  hRelUnc_tot->Draw("hist");
1288  hRelUncert[5]->Draw("hist same");
1289  hRelUncert[4]->Draw("hist same");
1290  hRelUncert[3]->Draw("hist same");
1291  hRelUncert[2]->Draw("hist same");
1292  hRelUncert[1]->Draw("hist same");
1293  hRelUncert[0]->Draw("hist same");
1294  hRelUnc_fit->Draw("hist same");
1295  hRelUnc_tot->Draw("hist same");
1296  TLegend* leg_systs = ana::AutoPlaceLegend(0.3,0.3);
1297  leg_systs->AddEntry(hRelUnc_tot, "Total", "l");
1298  leg_systs->AddEntry(hRelUnc_fit, "Fit", "l");
1299  leg_systs->AddEntry(hRelUncert[3], "Calibration Normalization", "l");
1300  leg_systs->AddEntry(hRelUncert[2], "Light Level", "l");
1301  leg_systs->AddEntry(hRelUncert[1], "Flux", "l");
1302  leg_systs->AddEntry(hRelUncert[0], "#nu-Interaction", "l");
1303  leg_systs->AddEntry(hRelUncert[4], "Cherenkov Variation", "l");
1304  leg_systs->AddEntry(hRelUncert[5], "Calibration Shape", "l");
1305  leg_systs->Draw();
1306  sprintf(name, "%s%s_%s_%s_%s_%s.png","AnalysisPlots/", dataName.c_str(),
1307  "xsec", "relative_uncertantiy",xsecResult.c_str(),"costheta");
1308  c2->SaveAs(name);
1309  }
1310 
1311 
1312 
1313 }
std::max
T max(const caf::Proxy< T > &a, T b)
Definition: BasicTypesProxy.h:416
name
const XML_Char * name
Definition: expat.h:151
kOrange
enum BeamMode kOrange
pot
const double pot
Definition: makeXSecPlots2D.C:17
kRed
enum BeamMode kRed
MakePlots
void MakePlots(TH2F *hFitted2D, std::vector< TH2F * > hSysts, TH2F *hTruth2D, std::string ytitle, std::string xtitle, std::string dataName, std::string xsecResult, std::string varName, bool calcCov=false)
Definition: makeXSecPlots2D.C:137
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
ana::ProjectionX
Eigen::ArrayXd ProjectionX(const Eigen::MatrixXd &mat)
Helper for Unweighted.
Definition: ReweightableSpectrum.cxx:99
nuebarccinc::eelecbins
const ana::Binning eelecbins
Definition: NuebarCCIncBins.h:19
ana::Spectrum::ToTH1
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
ana::chns
const SelDef chns[knumchns]
Definition: NCPi0SemiIncCuts.h:309
Integral
double Integral(const Spectrum &s, const double pot, int cvnregion)
Definition: syst_table_from_preds.C:46
nuebarccinc::costhetabins
const ana::Binning costhetabins
Definition: NuebarCCIncBins.h:23
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
plot_validation_datamc.Clone
def Clone(obj)
Definition: plot_validation_datamc.py:15
cet::pow
constexpr T pow(T x)
Definition: pow.h:75
sSignalEst
const std::string sSignalEst
Definition: makeXSecPlots2D.C:21
sDir
const std::string sDir
Definition: makeXSecPlots2D.C:20
ana::LoadFrom
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
demo5.c2
c2
Definition: demo5.py:33
make_syst_table_plots.ibin
ibin
Definition: make_syst_table_plots.py:246
make_true_q0q3_plots.caption
caption
Definition: make_true_q0q3_plots.py:404
nnucleons
const double nnucleons
Definition: makeXSecPlots2D.C:18
DoUnfolding
TH1F * DoUnfolding(TH1F *hMeasured, TH2F *hResponse, int iter)
Definition: NueCCIncCrossSectionFunctions.h:150
calib::j
const double j
Definition: BetheBloch.cxx:29
cmf::Spectrum
std::vector< float > Spectrum
Definition: Constants.h:728
ana::kBinContent
Regular histogram.
Definition: UtilsExt.h:30
om::cout
OStream cout
Definition: OStream.cxx:6
sAnalysis
const std::string sAnalysis
Definition: makeXSecPlots2D.C:22
nuebarccinc::enubins
const ana::Binning enubins
Definition: NuebarCCIncBins.h:14
ana::ProjectionY
Eigen::ArrayXd ProjectionY(const Eigen::MatrixXd &mat)
Helper for WeightingVariable.
Definition: ReweightableSpectrum.cxx:105
MakePlots1D
void MakePlots1D(TH2F *hData, TH2F *hFit, TH2F *hMC, std::string ytitle, std::string xtitle, std::string dataName, std::string xsecResult, std::string varName)
Definition: makeXSecPlots2D.C:25
ana::GenieMultiverseSpectra::UpperSigma
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal) const
Definition: GenieMultiverseSyst.cxx:388
ana::GenieMultiverseSpectra::LowerSigma
const Spectrum * LowerSigma(BandOptions opt=kBandFromNominal) const
Definition: GenieMultiverseSyst.cxx:648
dataName
std::string dataName
Definition: make_starplots.C:22
ana::AutoPlaceLegend
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:717
sSysts
const std::string sSysts
Definition: makeXSecPlots2D.C:23
ana::SelDef::name
std::string name
Definition: NCPi0SemiIncCuts.h:294
e
Float_t e
Definition: plot.C:35
kGreen
enum BeamMode kGreen
kBlue
enum BeamMode kBlue
ana::ToTH3
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: UtilsExt.cxx:162
ana::UniqueName
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
ana::ToTH2
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
For use with Var2D.
Definition: UtilsExt.cxx:115
PlotSystErrorBand
TGraphAsymmErrors * PlotSystErrorBand(TH1F *&nom, std::vector< TH1F * > &ups, std::vector< TH1F * > &dns, int col, int errCol, float headroom, bool newaxis)
Definition: NuebarCCIncTemplateFit.h:373
febshutoff_auto.val
val
Definition: febshutoff_auto.py:277
string
enum BeamMode string
uint
unsigned int uint
Definition: DAQChannelMapBaseOLD.h:22

Variable Documentation

const double nnucleons = 3.02225e+31

Definition at line 18 of file makeXSecPlots2D.C.

Referenced by makeXSecPlots2D().

const double pot = kAna2018FHCPOT

Definition at line 17 of file makeXSecPlots2D.C.

Referenced by makeXSecPlots2D().

const std::string sAnalysis = "SpectraNominal/fCrossSection_Nominal_FakeData_prongcvn_new.root"

Definition at line 22 of file makeXSecPlots2D.C.

Referenced by makeXSecPlots2D().

const std::string sDir = "/nova/ana/users/ddoyle/NuebarCrossSection_020319/"

Definition at line 20 of file makeXSecPlots2D.C.

Referenced by makeXSecPlots2D().

const std::string sSignalEst = "CrossSection_TemplateFitResults_fakedata_development_stat.root"

Definition at line 21 of file makeXSecPlots2D.C.

Referenced by makeXSecPlots2D().

const std::string sSysts = "SpectraSystematics/fCrossSection_EfficiencyFlux_Systs.root"

Definition at line 23 of file makeXSecPlots2D.C.

Referenced by makeXSecPlots2D().

Generated on Wed Feb 24 2021 14:15:34 for NOvA by   doxygen 1.8.11