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");
253 sprintf(name,
"mc_%s_%s_%s_postfit",
dataName.c_str(),
259 TH1F* hFitTotal1D = (TH1F*)hFitTotal3D->ProjectionZ(
"hFitTotal1D");
263 sprintf(name,
"mc_%s_%s_%s",
dataName.c_str(),
"analysis",
"mc");
268 TH1F* hDataTotal1D = (TH1F*)hDataTotal3D->ProjectionZ(
"hDataTotal1D");
271 sprintf(name,
"mc_%s_%s_%s",
"nominal",
"analysis",
"mc");
276 TH1F* hMCTotal1D = (TH1F*)hMCTotal3D->ProjectionZ(
"hMCTotal1D");
279 "Events/8.09 #times 10^{20} POT",
280 "Reconstructed Neutrino Energy, E_{#nu} (GeV)",
281 dataName, xsecResult,
"total_fit_comparison");
287 sprintf(name,
"mc_%s_%s_signal_postfit",
dataName.c_str(),
"tot");
292 TH1F* hSignalEst = (TH1F*)hSignalEst3D->ProjectionZ(
"hSignalEst");
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",
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());
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));
317 TH3F* hTrueRecoSignal3D = (TH3F*)
ana::ToTH3(sTrueRecoSignal,
321 TH1F* hTrueRecoSignal =
322 (TH1F*)hTrueRecoSignal3D->ProjectionZ(
"hTrueRecoSignal");
325 sprintf(name,
"mc_nominal_analysis_%s",
chns[2].name.c_str());
328 TH3F* hTrueRecoSignal_nofit3D = (TH3F*)
ana::ToTH3(sTrueRecoSignal_nofit,
332 TH1F* hTrueRecoSignal_nominal =
333 (TH1F*)hTrueRecoSignal_nofit3D->ProjectionZ(
"hTrueRecoSignal_nominal");
345 sprintf(name,
"%s_%s_%s_%s",
"mc",
"nominal",
"efficiency",
"num_1d");
347 sprintf(name,
"%s_%s_%s_%s",
"mc",
"nominal",
"efficiency",
"denom_1d");
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",
357 for(
uint idataset = 0; idataset < dataset_labels.size(); idataset++){
358 if(dataset_labels[idataset] ==
"genie" ||
359 dataset_labels[idataset] ==
"ppfx"){
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");
366 (fSysts->GetDirectory(name)));
367 sprintf(name,
"%s_%s_%i_%s_%s",
"mc",
368 dataset_labels[idataset].c_str(),i,
369 "efficiency",
"denom_1d");
371 (fSysts->GetDirectory(name)));
375 sprintf(name,
"%s_%s_%s_%s",
"mc", dataset_labels[idataset].c_str(),
376 "efficiency",
"num_1d");
378 (fSysts->GetDirectory(name)));
379 sprintf(name,
"%s_%s_%s_%s",
"mc", dataset_labels[idataset].c_str(),
380 "efficiency",
"denom_1d");
382 (fSysts->GetDirectory(name)));
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);
393 TH1F* hHolder2 = (TH1F*)vSystsEfficiencyDenom[i]->ToTH1(
pot);
395 hHolder->Divide(hHolder2);
396 vSystsEff.push_back(hHolder);
399 TH1F* hLightDwnEfficiency =
401 TH1F* hLightUpEfficiency =
403 TH1F* hCkvEfficiency =
405 TH1F* hCalibDwnEfficiency =
407 TH1F* hCalibUpEfficiency =
409 TH1F* hCalibShapeEfficiency =
412 std::vector<Spectrum*> spects_genie;
413 std::vector<std::unique_ptr<Spectrum>> genie_hists;
414 for(
int i = 6; i < 100+6; i++){
416 spects_genie.push_back(spect);
417 genie_hists.push_back((std::unique_ptr<Spectrum>)spects_genie[i-6]);
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++){
423 spects_ppfx.push_back(spect);
424 ppfx_hists.push_back((std::unique_ptr<Spectrum>)spects_ppfx[i-106]);
435 sprintf(name,
"%s_%s_%s",
"mc",
"nominal",
"flux");
437 TH1F* hFlux = (TH1F*)sflux_nom.
ToTH1(
pot);
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");
444 (fSysts->GetDirectory(name)));
445 vSystsFlux.push_back((TH1F*)vsSystsFlux[i]->ToTH1(
pot));
446 vSystsFlux[
i]->Scale(1
e-4);
452 sprintf(name,
"%s_%s_%s_%s",
"mc",
dataName.c_str(),
"unfolding",
"1d");
453 TH2F* hResponseMatrix =
457 TH1F* hUnfolded_SignalEst =
DoUnfolding(hSignalEst, hResponseMatrix,1);
458 std::vector<TH1F*> hUnfoldedSysts;
459 for(
int i = 0; i < (
int)hSystsReco1D.size(); i++){
462 hUnfoldedSysts.push_back((TH1F*)
DoUnfolding(hSystsReco1D[i],
469 TH1F* hXsec_cv = (TH1F*)hUnfolded_SignalEst->Clone(
ana::UniqueName().c_str());
470 hXsec_cv->SetTitle(
"Cross Section");
473 hXsec_cv->Divide(cv_eff);
474 hXsec_cv->Divide(hFlux);
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();
481 std::vector<TH1F*>hXSec_systs;
482 for(
int i = 0; i < (
int)vSystsEff.size(); i++){
487 hHolder->Divide(vSystsEff[i]);
491 hHolder->Divide(vSystsFlux[i-106]);
494 hHolder->Divide(hFlux);
495 hHolder->SetLineColor(
kBlue);
496 hXSec_systs.push_back(hHolder);
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);
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++){
517 spects_xsec_genie.push_back(spect);
518 vspects_xsec_genie.push_back((std::unique_ptr<Spectrum>)
519 spects_xsec_genie[i-6]);
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++){
526 spects_xsec_ppfx.push_back(spect);
527 vspects_xsec_ppfx.push_back((std::unique_ptr<Spectrum>)
528 spects_xsec_ppfx[i-106]);
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);
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);
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);
560 std::vector<TH1F*> systs_up =
568 std::vector<TH1F*> systs_down =
578 TGraphAsymmErrors*
g =
584 sprintf(name,
"%s_%s_%s_%s",
"mc",
dataName.c_str(),
"efficiency",
"denom_1d");
587 hTrueSignal->Divide(hFlux);
589 hTrueSignal->SetLineColor(
kRed);
593 hXsec_cv->SetMarkerColor(kBlack);
594 hXsec_cv->SetLineColor(kBlack);
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);
606 double err_fluxeff = (err_hi > err_lo)? err_hi/total : err_lo/total;
609 double tot_err = total*
sqrt(
pow(uncert/selected_tot,2) +
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) +
622 TCanvas* cResult =
new TCanvas(
"cResult",
"cResult");
623 hUncert_tot->GetXaxis()->SetRangeUser(1.0,10.0);
625 hUncert_tot->SetMarkerStyle(20);
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");
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");
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");
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);
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);
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.);
674 Xnue280->SetLineWidth(2);
675 Xnue280->SetLineColor(4);
678 TCanvas*
c3 =
new TCanvas(
"c3",
"cResult");
679 hUncert_tot->GetXaxis()->SetRangeUser(1.0,10.0);
680 hUncert_tot->SetMarkerStyle(20);
682 Xnue280->Draw(
"l same");
684 hTrueSignal->Draw(
"hist same");
685 hUncert_tot->Draw(
"same");
687 sprintf(name,
"%s%s_%s_%s_%s.png",
"AnalysisPlots/",
dataName.c_str(),
688 "xsec",
"results_truth_compare",xsecResult.c_str());
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);
702 std::vector<TH1F*> hRelUncert;
703 for(
int i = 0; i < (
int)systs_up.size();i++){
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);
714 hHolder->SetBinContent(
ibin,uncert/val_cv);
716 hRelUncert.push_back(hHolder);
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);
730 hRelUncert[5]->SetLineColor(
kGreen+2);
731 hRelUncert[4]->SetLineColor(
kOrange-3);
732 hRelUncert[0]->SetLineColor(
kRed - 2);
733 hRelUncert[1]->SetLineColor(kCyan-3);
734 hRelUncert[2]->SetLineColor(
kBlue-4);
735 hRelUncert[3]->SetLineColor(kMagenta+1);
736 hRelUnc_fit->SetLineColor(kBlack);
737 hRelUnc_fit->SetLineStyle(7);
738 hRelUnc_tot->SetLineColor(kBlack);
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");
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");
764 sprintf(name,
"%s%s_%s_%s_%s.png",
"AnalysisPlots/",
dataName.c_str(),
765 "xsec",
"relative_uncertantiy",xsecResult.c_str());
T max(const caf::Proxy< T > &a, T b)
fvar< T > fabs(const fvar< T > &x)
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.
const SelDef chns[knumchns]
const ana::Binning costhetabins
const std::string sSignalEst
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Representation of a spectrum in any variable, with associated POT.
TH1F * DoUnfolding(TH1F *hMeasured, TH2F *hResponse, int iter)
std::vector< float > Spectrum
const std::string sAnalysis
const ana::Binning enubins
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal) const
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
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
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.
std::string UniqueName()
Return a different string each time, for creating histograms.
TGraphAsymmErrors * PlotSystErrorBand(TH1F *&nom, std::vector< TH1F * > &ups, std::vector< TH1F * > &dns, int col, int errCol, float headroom, bool newaxis)