22 #include "TGraphAsymmErrors.h" 33 std::vector<ColumnDef>
cols;
36 const bool beamLab =
beam ==
"fhc" ?
true :
false;
38 std::vector<std::vector<const ISyst *>>
systs;
39 std::vector<std::string> extrapStr;
40 std::vector<std::string> extrapLab;
45 extrapStr.push_back(
"noextrap_noPt");
46 extrapLab.push_back(
"Not Extrapolated");
74 extrapStr.push_back(
"extrap");
75 extrapLab.push_back(
"Extrapolated");
96 extrapStr.push_back(
"combo");
97 extrapLab.push_back(
"Extrapolated");
116 extrapStr.push_back(
"extrap");
117 extrapLab.push_back(
"Extrapolated");
138 extrapStr.push_back(
"prop");
139 extrapLab.push_back(
"Extrapolated");
142 gSystem->MakeDirectory(
"plots" );
144 const unsigned int cosmicsIdx = cols.size()-1;
150 const double dCP = .82;
151 const double ssth23 = .568;
153 const double dmsq32 = 2.406e-3;
156 calcosc->
SetdCP ( dCP*(TMath::Pi()) );
161 std::vector<std::map<std::string,std::pair<double,double>>> tots, sigs, bkgs;
162 std::vector<std::pair<double,double>> totQuad, sigQuad, bkgQuad;
164 for(
unsigned int ext = 0;
ext < extrapStr.size(); ++
ext){
165 std::vector< std::vector<Spectrum> > noms;
166 std::vector<const IPrediction*>
pred;
167 std::vector< std::pair <Spectrum*, double> >
cosmics;
177 for (
size_t pr=0;
pr<pred.size(); ++
pr) {
178 std::vector<Spectrum> Tempnom;
180 if(
col.name ==
"Cosmics")
182 Tempnom.push_back(*cosmics[
pr].first);
187 Tempnom.push_back(pred[
pr]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign));
191 noms.push_back(Tempnom);
198 std::vector<Spectrum> Tempnom;
200 if(
col.name ==
"Cosmics")
202 Spectrum sumCosmics = *cosmics[0].first + *cosmics[1].first + *cosmics[2].first + *cosmics[3].first;
203 Tempnom.push_back(sumCosmics);
208 Spectrum sumBeam = pred[0]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign) +
209 pred[1]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign) +
210 pred[2]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign) +
211 pred[3]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign);
212 Tempnom.push_back(sumBeam);
216 noms.push_back(Tempnom);
227 unsigned int predSize = pred.size();
230 if(!
IsNuMu && quantsComp) predSize += 3;
232 for (
unsigned int prIdx=0; prIdx<predSize; ++prIdx) {
233 unsigned int pr = prIdx;
234 if(!
IsNuMu && quantsComp) pr = 0;
236 TH1* hNomTot = noms[
pr][0].ToTH1(pot);
238 hNomTot->Add(noms[pr][cosmicsIdx].ToTH1(livetime,
kLivetime));
240 TH1* hNomSig = noms[
pr][1].ToTH1(pot);
245 noms[
pr][2] -= noms[
pr][1];
247 noms[
pr][3] -= noms[
pr][1];
250 TH1* hNomBkg = noms[
pr][2].ToTH1(pot);
251 hNomBkg->Add(noms[pr][cosmicsIdx].ToTH1(livetime,
kLivetime));
253 std::map<std::string, TH1*> hUpsTot, hUpsSig, hUpsBkg;
254 std::map<std::string, TH1*> hDnsTot, hDnsSig, hDnsBkg;
267 sSamp +=
" All Quartiles";
277 std::cout<<
"\\multicolumn{19}{c}{\\bf{Systematic Summary for " << sSamp <<
" (\\%)}} \\\\"<<
std::endl;
280 if(prIdx==1) sBin =
"Low PID";
281 else if(prIdx==2) sBin =
"High PID";
282 else if(prIdx==3) sBin =
"Peripheral";
288 std::cout<<
"\\multicolumn{21}{c}{\\bf{Systematic Summary for " << sSamp <<
" (\\%)}} \\\\"<<
std::endl;
296 std::pair<double,double> quads[cols.size()];
297 for(
unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx) quads[colIdx] = {0,0};
299 std::map<std::string, std::pair<double,double>> barsTot, barsSig, barsBkg;
305 std::vector<Spectrum> shiftsUp, shiftsDn;
308 if(!
IsNuMu || pr != pred.size())
310 if(
col.name !=
"Cosmics")
314 shiftsUp.push_back(pred[pr]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign));
315 shiftsDn.push_back(pred[pr]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign));
319 shiftsUp.push_back(pred[pr]->PredictComponentSyst(calcosc,
SystShifts(syst, +1),
col.flav,
col.curr,
col.sign));
320 shiftsDn.push_back(pred[pr]->PredictComponentSyst(calcosc,
SystShifts(syst, -1),
col.flav,
col.curr,
col.sign));
325 shiftsUp.push_back(*cosmics[pr].first);
326 shiftsDn.push_back(*cosmics[pr].first);
338 if(
col.name !=
"Cosmics")
342 shiftsUp.push_back(pred[0]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign)+
343 pred[1]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign)+
344 pred[2]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign)+
345 pred[3]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign));
346 shiftsDn.push_back(pred[0]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign)+
347 pred[1]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign)+
348 pred[2]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign)+
349 pred[3]->PredictComponent(calcosc,
col.flav,
col.curr,
col.sign));
353 shiftsUp.push_back(pred[0]->PredictComponentSyst(calcosc,
SystShifts(syst, +1),
col.flav,
col.curr,
col.sign)+
354 pred[1]->PredictComponentSyst(calcosc,
SystShifts(syst, +1),
col.flav,
col.curr,
col.sign)+
355 pred[2]->PredictComponentSyst(calcosc,
SystShifts(syst, +1),
col.flav,
col.curr,
col.sign)+
356 pred[3]->PredictComponentSyst(calcosc,
SystShifts(syst, +1),
col.flav,
col.curr,
col.sign));
357 shiftsDn.push_back(pred[0]->PredictComponentSyst(calcosc,
SystShifts(syst, -1),
col.flav,
col.curr,
col.sign)+
358 pred[1]->PredictComponentSyst(calcosc,
SystShifts(syst, -1),
col.flav,
col.curr,
col.sign)+
359 pred[2]->PredictComponentSyst(calcosc,
SystShifts(syst, -1),
col.flav,
col.curr,
col.sign)+
360 pred[3]->PredictComponentSyst(calcosc,
SystShifts(syst, -1),
col.flav,
col.curr,
col.sign));
377 shiftsUp.push_back(cosmicsUp);
378 shiftsDn.push_back(cosmicsDn);
382 shiftsUp.push_back(*cosmics[0].first+*cosmics[1].first+*cosmics[2].first+*cosmics[3].first);
383 shiftsDn.push_back(*cosmics[0].first+*cosmics[1].first+*cosmics[2].first+*cosmics[3].first);
428 shiftsUp[2] -= shiftsUp[1];
429 shiftsDn[2] -= shiftsDn[1];
432 shiftsUp[3] -= shiftsUp[1];
433 shiftsDn[3] -= shiftsDn[1];
435 for(
unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx){
437 const Spectrum& cosNom = noms[
pr][cosmicsIdx];
439 double ratioUp, errUp, ratioDn, errDn;
440 if(colIdx == 0 || colIdx == 2)
443 errUp = 100.0*(ratioUp-1.0);
446 errDn = 100.0*(ratioDn-1.0);
448 else if(colIdx == cosmicsIdx)
451 errUp = 100.0*(ratioUp-1.0);
454 errDn = 100.0*(ratioDn-1.0);
459 errUp = 100.0*(ratioUp-1.0);
462 errDn = 100.0*(ratioDn-1.0);
466 if(!
IsNuMu && quantsComp && prIdx != 0)
468 std::pair<int,int> Int = {1,9};
469 if(prIdx==2) Int = {10,18};
470 if(prIdx==3) Int = {19,hNomTot->GetNbinsX()};
472 if(colIdx == 0 || colIdx == 2)
474 ratioUp = 1.0*( shiftsUp[colIdx].ToTH1(pot)->Integral(Int.first,Int.second) + shiftsUp[cosmicsIdx].ToTH1(livetime,
kLivetime)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.
ToTH1(pot)->Integral(Int.first,Int.second) + cosNom.
ToTH1(livetime,
kLivetime)->Integral(Int.first,Int.second) ) );
475 errUp = 100.0*(ratioUp-1.0);
477 ratioDn = 1.0*( shiftsDn[colIdx].ToTH1(pot)->Integral(Int.first,Int.second) + shiftsDn[cosmicsIdx].ToTH1(livetime,
kLivetime)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.
ToTH1(pot)->Integral(Int.first,Int.second) + cosNom.
ToTH1(livetime,
kLivetime)->Integral(Int.first,Int.second) ) );
478 errDn = 100.0*(ratioDn-1.0);
480 else if(colIdx == cosmicsIdx)
482 ratioUp = 1.0*( shiftsUp[colIdx].ToTH1(livetime,
kLivetime)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.
ToTH1(livetime,
kLivetime)->Integral(Int.first,Int.second) ) );
483 errUp = 100.0*(ratioUp-1.0);
485 ratioDn = 1.0*( shiftsDn[colIdx].ToTH1(livetime,
kLivetime)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.
ToTH1(livetime,
kLivetime)->Integral(Int.first,Int.second) ) );
486 errDn = 100.0*(ratioDn-1.0);
490 ratioUp = 1.0*( shiftsUp[colIdx].ToTH1(pot)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.
ToTH1(pot)->Integral(Int.first,Int.second) ) );
491 errUp = 100.0*(ratioUp-1.0);
493 ratioDn = 1.0*( shiftsDn[colIdx].ToTH1(pot)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.
ToTH1(pot)->Integral(Int.first,Int.second) ) );
494 errDn = 100.0*(ratioDn-1.0);
500 if(max < 0.01) max = 0;
501 if(min > -0.01) min = 0;
503 quads[colIdx].first += min*
min;
504 quads[colIdx].second += max*
max;
519 std::cout <<
" & \\cellcolor{" << minstr <<
"}" << (min == 0 ?
"-" :
"") << min <<
" & \\cellcolor{"<< maxstr <<
"}+" << max;
523 if(cols[colIdx].
name ==
"$\\nu_e$ Signal" ||
524 cols[colIdx].
name ==
"$\\bar\\nu_e$ Signal" ||
525 cols[colIdx].
name ==
"$\\nu_\\mu$ Signal" ||
526 cols[colIdx].
name ==
"$\\bar\\nu_\\mu$ Signal")
529 if(cols[colIdx].
name ==
"Total Pred")
532 if(cols[colIdx].
name ==
"Total Bkg")
542 for(
unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx)
std::cout <<
" & -" <<
sqrt(quads[colIdx].first) <<
" & +" <<
sqrt(quads[colIdx].
second);
558 std::cout<<
"\\multicolumn{7}{c}{\\bf{Summed Systematic for " << sSamp <<
" (\\%)}} \\\\"<<
std::endl;
560 for(
int colIdx = 0; colIdx < 3; ++colIdx)
std::cout <<
" & \\multicolumn{2}{c}{" << cols[colIdx].
name <<
"}";
575 for(
int colIdx = 0; colIdx < 3; ++colIdx)
std::cout <<
" & -" <<
sqrt(quads[colIdx].first) <<
" & +" <<
sqrt(quads[colIdx].second);
584 tots.push_back(barsTot);
585 sigs.push_back(barsSig);
586 bkgs.push_back(barsBkg);
588 for(
unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx){
589 quads[colIdx].first =
sqrt(quads[colIdx].first);
590 quads[colIdx].second =
sqrt(quads[colIdx].second);
593 totQuad.push_back(quads[0]);
594 sigQuad.push_back(quads[1]);
595 bkgQuad.push_back(quads[2]);
601 if(pr==pred.size()) sQuant =
"_AllQuartiles";
605 sQuant =
"_AllPrediction";
606 if(prIdx==1) sQuant =
"_lowPID";
607 else if(prIdx==2) sQuant =
"_highPID";
608 else if(prIdx==3) sQuant =
"_peripheral";
613 sampLab =
"#nu_{e} Selection";
614 if(
beam ==
"rhc") sampLab =
"#bar{#nu}_{e} Selection";
619 sampLab =
"#nu_{#mu} Selection "+quart;
620 if(
beam ==
"rhc") sampLab =
"#bar{#nu}_{#mu} Selection "+quart;
624 if(prIdx==1) sampLab +=
" Low PID";
625 if(prIdx==2) sampLab +=
" High PID";
626 if(prIdx==3) sampLab +=
" Peripheral";
630 TH1* hD = (TH1*)hNomTot->Clone();
631 TH1* hNomTotRatio = (TH1*)hNomTot->Clone();
632 hNomTotRatio->Divide(hNomTot);
633 std::vector<TH1*> tempUp, tempDn;
634 for(
auto it: hUpsTot) tempUp.push_back(
it.second);
635 for(
auto it: hDnsTot) tempDn.push_back(
it.second);
637 std::vector<TH1*> hUpsTotRatios =
MakeRatios(hNomTot, tempUp);
638 std::vector<TH1*> hDnsTotRatios =
MakeRatios(hNomTot, tempDn);
640 std::map<std::string, std::vector<TH1*>> hUpsTotType =
SortSystHists(hUpsTot);
641 std::map<std::string, std::vector<TH1*>> hDnsTotType =
SortSystHists(hDnsTot);
643 std::vector<int>
types = {0,1,2,4,5};
644 if(
IsNuMu) types = {1,2,4,5,6};
645 TLegend *lTypes =
new TLegend(.1,.1,.9,.9);
646 lTypes->SetFillStyle(0);
648 TCanvas *
c =
new TCanvas(TString(
"c")+TString(sQuant),
"c");
651 lTypes->AddEntry(hD,
"Total syst. uncertainty",
"l");
659 TH1* hDummy = (TH1*)hNomTot->Clone();
661 lTypes->AddEntry(hDummy,
systTypes[j].c_str(),
"l");
667 gPad->Print((
"plots/ratioband_"+sNuMu+
"_"+extrapStr[ext]+
"_"+sumStr+
"_"+
beam+sQuant+
"_tot.pdf").c_str());
669 TCanvas *cLeg =
new TCanvas(TString(
"cLeg")+TString(sQuant),
"cLeg");
672 gPad->Print((
"plots/ratioband_"+sNuMu+
"_"+extrapStr[ext]+
"_"+sumStr+
"_"+
beam+sQuant+
"_legend.pdf").c_str());
675 TCanvas *cTot =
new TCanvas(TString(
"cTot")+TString(sQuant),
"cTot");
676 ErrorBarChart(barsTot,{totErr,totErr},
"Total Prediction Uncertainty (%)");
677 gPad->SetLeftMargin(.35);
681 gPad->Print((
"plots/syst_summary_"+sNuMu+
"_"+extrapStr[ext]+
"_"+sumStr+
"_"+
beam+sQuant+
"_tot.pdf").c_str());
683 TCanvas *cTot2 =
new TCanvas(TString(
"cTot2")+TString(sQuant),
"cTot2");
684 ErrorBarChart(barsTot,{0,0},
"Total Prediction Uncertainty (%)");
685 gPad->SetLeftMargin(.35);
689 gPad->Print((
"plots/syst_summary_"+sNuMu+
"_"+extrapStr[ext]+
"_"+sumStr+
"_"+
beam+sQuant+
"_tot_nostat.pdf").c_str());
691 TCanvas *cSig =
new TCanvas(TString(
"cSig")+TString(sQuant),
"cSig");
692 ErrorBarChart(barsSig,{sigErr,sigErr},
"Signal Uncertainty (%)");
693 gPad->SetLeftMargin(.35);
697 gPad->Print((
"plots/syst_summary_"+sNuMu+
"_"+extrapStr[ext]+
"_"+sumStr+
"_"+
beam+sQuant+
"_sig.pdf").c_str());
699 TCanvas *cSig2 =
new TCanvas(TString(
"cSig2")+TString(sQuant),
"cSig2");
701 gPad->SetLeftMargin(.35);
705 gPad->Print((
"plots/syst_summary_"+sNuMu+
"_"+extrapStr[ext]+
"_"+sumStr+
"_"+
beam+sQuant+
"_sig_nostat.pdf").c_str());
707 TCanvas *cBkg =
new TCanvas(TString(
"cBkg")+TString(sQuant),
"cBkg");
708 ErrorBarChart(barsBkg,{bkgErr,bkgErr},
"Background Uncertainty (%)");
709 gPad->SetLeftMargin(.35);
713 gPad->Print((
"plots/syst_summary_"+sNuMu+
"_"+extrapStr[ext]+
"_"+sumStr+
"_"+
beam+sQuant+
"_bkg.pdf").c_str());
715 TCanvas *cBkg2 =
new TCanvas(TString(
"cBkg2")+TString(sQuant),
"cBkg2");
717 gPad->SetLeftMargin(.35);
721 gPad->Print((
"plots/syst_summary_"+sNuMu+
"_"+extrapStr[ext]+
"_"+sumStr+
"_"+
beam+sQuant+
"_bkg_nostat.pdf").c_str());
732 unsigned int nQuants = 4;
734 for(
unsigned int i = 0;
i < extrapStr.size(); ++
i)
736 std::vector<std::string> quantLab = {
"Quartile 1",
"Quartile 2",
"Quartile 3",
"Quartile 4"};
738 std::vector<std::string>
labels;
739 std::vector<std::vector<std::pair<double,double>>> fTot;
740 std::vector<std::vector<std::pair<double,double>>>
fSig;
741 std::vector<std::vector<std::pair<double,double>>> fBkg;
742 std::vector<std::pair<double,double>> fTotQuad;
743 std::vector<std::pair<double,double>> fSigQuad;
744 std::vector<std::pair<double,double>> fBkgQuad;
746 for(
unsigned int q = 0; q < nQuants; ++q)
748 std::vector<std::pair<double,double>> tempTot;
749 std::vector<std::pair<double,double>> tempSig;
750 std::vector<std::pair<double,double>> tempBkg;
751 int idx =
i*(nQuants+1) + q;
759 tempTot.push_back(tots[idx][
systTypes[j]]);
760 tempSig.push_back(sigs[idx][
systTypes[j]]);
761 tempBkg.push_back(bkgs[idx][
systTypes[j]]);
766 for(
const ISyst* syst: systs[0]){
775 fTot.push_back(tempTot);
776 fSig.push_back(tempSig);
777 fBkg.push_back(tempBkg);
778 fTotQuad.push_back(totQuad[idx]);
779 fSigQuad.push_back(sigQuad[idx]);
780 fBkgQuad.push_back(bkgQuad[idx]);
783 sampLab =
"#nu_{#mu} Selection";
784 if(
beam ==
"rhc") sampLab =
"#bar{#nu}_{#mu} Selection";
788 myBarChart(
"Total Prediction Uncertainty (%)", labels,
789 fTot, fTotQuad, quantLab, colors);
790 gPad->SetLeftMargin(.35);
795 gPad->Print((
"plots/syst_quant_comp_"+sNuMu+
"_"+sumStr+
"_"+
beam+
"_tot_"+extrapStr[
i]+
".pdf").c_str());
799 fSig, fSigQuad, quantLab, colors);
800 gPad->SetLeftMargin(.35);
805 gPad->Print((
"plots/syst_quant_comp_"+sNuMu+
"_"+sumStr+
"_"+
beam+
"_sig_"+extrapStr[
i]+
".pdf").c_str());
808 myBarChart(
"Background Uncertainty (%)", labels,
809 fBkg, fBkgQuad, quantLab, colors);
810 gPad->SetLeftMargin(.35);
815 gPad->Print((
"plots/syst_quant_comp_"+sNuMu+
"_"+sumStr+
"_"+
beam+
"_bkg_"+extrapStr[
i]+
".pdf").c_str());
822 unsigned int nBins = 3;
824 for(
unsigned int i = 0;
i < extrapStr.size(); ++
i)
826 std::vector<std::string> quantLab = {
"Low PID",
"High PID",
"Peripheral"};
828 std::vector<std::string>
labels;
829 std::vector<std::vector<std::pair<double,double>>> fTot;
830 std::vector<std::vector<std::pair<double,double>>>
fSig;
831 std::vector<std::vector<std::pair<double,double>>> fBkg;
832 std::vector<std::pair<double,double>> fTotQuad;
833 std::vector<std::pair<double,double>> fSigQuad;
834 std::vector<std::pair<double,double>> fBkgQuad;
836 for(
unsigned int b = 1;
b <=
nBins; ++
b)
838 std::vector<std::pair<double,double>> tempTot;
839 std::vector<std::pair<double,double>> tempSig;
840 std::vector<std::pair<double,double>> tempBkg;
841 int idx =
i*(nBins+1) +
b;
849 tempTot.push_back(tots[idx][
systTypes[j]]);
850 tempSig.push_back(sigs[idx][
systTypes[j]]);
851 tempBkg.push_back(bkgs[idx][
systTypes[j]]);
856 for(
const ISyst* syst: systs[0]){
865 fTot.push_back(tempTot);
866 fSig.push_back(tempSig);
867 fBkg.push_back(tempBkg);
868 fTotQuad.push_back(totQuad[idx]);
869 fSigQuad.push_back(sigQuad[idx]);
870 fBkgQuad.push_back(bkgQuad[idx]);
873 sampLab =
"#nu_{e} Selection";
874 if(
beam ==
"rhc") sampLab =
"#bar{#nu}_{e} Selection";
878 myBarChart(
"Total Prediction Uncertainty (%)", labels,
879 fTot, fTotQuad, quantLab, colors);
880 gPad->SetLeftMargin(.35);
885 gPad->Print((
"plots/syst_bins_comp_"+sNuMu+
"_"+sumStr+
"_"+
beam+
"_tot_"+extrapStr[
i]+
".pdf").c_str());
889 fSig, fSigQuad, quantLab, colors);
890 gPad->SetLeftMargin(.35);
895 gPad->Print((
"plots/syst_bins_comp_"+sNuMu+
"_"+sumStr+
"_"+
beam+
"_sig_"+extrapStr[
i]+
".pdf").c_str());
898 myBarChart(
"Background Uncertainty (%)", labels,
899 fBkg, fBkgQuad, quantLab, colors);
900 gPad->SetLeftMargin(.35);
905 gPad->Print((
"plots/syst_bins_comp_"+sNuMu+
"_"+sumStr+
"_"+
beam+
"_bkg_"+extrapStr[
i]+
".pdf").c_str());
911 unsigned int nQuants = 1;
912 if(!
IsNuMu && quantsComp) nQuants = 4;
915 for(
unsigned int q = 0; q < nQuants; ++q)
918 std::vector<std::string>
labels;
919 std::vector<std::vector<std::pair<double,double>>> fTot;
920 std::vector<std::vector<std::pair<double,double>>>
fSig;
921 std::vector<std::vector<std::pair<double,double>>> fBkg;
922 std::vector<std::pair<double,double>> fTotQuad;
923 std::vector<std::pair<double,double>> fSigQuad;
924 std::vector<std::pair<double,double>> fBkgQuad;
926 for(
unsigned int i = 0;
i < extrapStr.size(); ++
i)
928 std::vector<std::pair<double,double>> tempTot;
929 std::vector<std::pair<double,double>> tempSig;
930 std::vector<std::pair<double,double>> tempBkg;
931 int idx =
i*nQuants + q;
939 tempTot.push_back(tots[idx][
systTypes[j]]);
940 tempSig.push_back(sigs[idx][
systTypes[j]]);
941 tempBkg.push_back(bkgs[idx][
systTypes[j]]);
947 for(
const ISyst* syst: systs[1]){
949 if (
i==0 &&
FixSystName(syst->ShortName()).Contains(
"LeptonAngle") )
951 tempTot.push_back({.0,.0});
952 tempSig.push_back({.0,.0});
953 tempBkg.push_back({.0,.0});
955 else if (
i==0 &&
FixSystName(syst->ShortName()).Contains(
"accept") )
959 tempTot.push_back(tots[idx][ (
std::string)
"accept signalkin FHC 2020" ]);
960 tempSig.push_back(sigs[idx][ (
std::string)
"accept signalkin FHC 2020" ]);
961 tempBkg.push_back(bkgs[idx][ (
std::string)
"accept signalkin FHC 2020" ]);
965 tempTot.push_back(tots[idx][ (
std::string)
"accept signalkin RHC 2020" ]);
966 tempSig.push_back(sigs[idx][ (
std::string)
"accept signalkin RHC 2020" ]);
967 tempBkg.push_back(bkgs[idx][ (
std::string)
"accept signalkin RHC 2020" ]);
979 fTot.push_back(tempTot);
980 fSig.push_back(tempSig);
981 fBkg.push_back(tempBkg);
982 fTotQuad.push_back(totQuad[idx]);
983 fSigQuad.push_back(sigQuad[idx]);
984 fBkgQuad.push_back(bkgQuad[idx]);
990 if(q==4) sQuant =
"_AllQuartiles";
993 sQuant =
"_AllPrediction";
994 if(q==1) sQuant =
"_lowPID";
995 if(q==2) sQuant =
"_highPID";
996 if(q==3) sQuant =
"_peripheral";
999 sampLab =
"#nu_{e} Selection";
1000 if(
beam ==
"rhc") sampLab =
"#bar{#nu}_{e} Selection";
1005 sampLab =
"#nu_{#mu} Selection "+quart;
1006 if(
beam ==
"rhc") sampLab =
"#bar{#nu}_{#mu} Selection "+quart;
1010 if(q==1) sampLab +=
" Low PID";
1011 if(q==2) sampLab +=
" High PID";
1012 if(q==3) sampLab +=
" Peripheral";
1016 myBarChart(
"Total Prediction Uncertainty (%)", labels,
1017 fTot, fTotQuad, extrapLab, {-1});
1018 gPad->SetLeftMargin(.35);
1023 gPad->Print((
"plots/syst_extrap_comp_"+sNuMu+
"_"+sumStr+
"_"+
beam+
"_tot"+sQuant+
".pdf").c_str());
1027 fSig, fSigQuad, extrapLab, {-1});
1028 gPad->SetLeftMargin(.35);
1033 gPad->Print((
"plots/syst_extrap_comp_"+sNuMu+
"_"+sumStr+
"_"+
beam+
"_sig"+sQuant+
".pdf").c_str());
1036 myBarChart(
"Background Uncertainty (%)", labels,
1037 fBkg, fBkgQuad, extrapLab, {-1});
1038 gPad->SetLeftMargin(.35);
1043 gPad->Print((
"plots/syst_extrap_comp_"+sNuMu+
"_"+sumStr+
"_"+
beam+
"_bkg"+sQuant+
".pdf").c_str());
T max(const caf::Proxy< T > &a, T b)
void systematics_summary_from_pred_interp(std::string beam="rhc", std::string extrapStr="prop", bool sum=true)
Cuts and Vars for the 2020 FD DiF Study.
void drawBeamLabel(bool isFHC, double x=0.35)
std::map< std::string, std::vector< TH1 * > > SortSystHists(std::map< std::string, TH1 * > shifts)
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.
Simple record of shifts applied to systematic parameters.
std::string systTypes[kNumSystTypes]
virtual void SetTh13(const T &th13)=0
std::map< std::string, std::pair< double, double > > SumSysts(std::map< std::string, std::pair< double, double >> bars)
Spectrum * ShiftedCosmics(Spectrum *cosmic, double sigma)
std::vector< const IPrediction * > GetNumuPredictions2020(const int nq=4, std::string decomp="noPt", osc::IOscCalc *calc=DefaultOscCalc(), bool useSysts=true, std::string beam="fhc", bool isFakeData=false, bool GetFromUPS=false, bool minimizeMemory=true, bool NERSC=false)
const Color_t kTotalMCErrorBandColor
const CosmicBkgScaleSyst kCosmicBkgScaleSyst
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Encapsulate code to systematically shift a caf::SRProxy.
virtual void SetDmsq32(const T &dmsq32)=0
Representation of a spectrum in any variable, with associated POT.
std::pair< Spectrum *, double > GetNueCosmics2020(std::string beam, bool GetFromUPS=false, bool NERSC=false)
void documentHeader(bool IsNuMu)
Charged-current interactions.
Interactions of both types.
std::vector< TH1 * > MakeRatios(TH1 *nom, std::vector< TH1 * > shifts)
const IPrediction * GetNuePrediction2020(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=true, bool NERSC=false)
const double kAna2020FHCLivetime
const double kAna2020FHCPOT
const double kAna2020RHCPOT
std::vector< std::pair< Spectrum *, double > > GetNumuCosmics2020(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
static float min(const float a, const float b, const float c)
void drawSampleLabel(std::string s, double x=.915, double y=.9)
void MyPlotWithSystErrorBand(TH1 *&nom, std::vector< TH1 * > &ups, std::vector< TH1 * > &dns, int col, int errCol, float linewidth, bool newaxis)
double Integral(const Spectrum &s, double pot)
const double kAna2020RHCLivetime
virtual void SetTh23(const T &th23)=0
TString FixSystName(TString mystring)
Neutral-current interactions.
std::vector< const ISyst * > get3FlavorAna2020AllSysts(const EAnaType2020 ana, const bool smallgenie, const BeamType2020 beam, const bool isFit, const bool ptExtrap)
void ErrorBarChart(const std::map< std::string, std::pair< double, double >> &systErrors, const std::pair< double, double > &statErr, const std::string &label)
Make a simple plot of relative size of different errors.
Both neutrinos and antineutrinos.
int systCols[kNumSystTypes]
T min(const caf::Proxy< T > &a, T b)
static void Add(TH3D *h, const int bx, const int by, const int bz, const double w)
All neutrinos, any flavor.
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
std::string to_string(ModuleType const mt)
const DummyRockScaleSyst kRockScaleSyst
virtual void SetdCP(const T &dCP)=0