19 TDirectory *
tmp = gDirectory;
20 TCanvas *
c =
new TCanvas();
21 c->SetLeftMargin(0.15);
22 c->SetRightMargin(0.01);
27 hist->Draw(draw_option.c_str());
28 hist->SetTitle(title.c_str());
29 c->Print(name.c_str());
36 std::vector<std::string>
labels,
41 TDirectory *
tmp = gDirectory;
42 TCanvas *
c =
new TCanvas();
43 c->SetLeftMargin(0.15);
44 c->SetRightMargin(0.01);
52 TLegend *
leg =
new TLegend(0.4,0.7,0.7,0.9);
54 double max_val = -1e9;
55 for(
auto ihist = 0
u; ihist < hists.size(); ihist++) {
56 if(hists[ihist] == NULL)
continue;
61 for(
auto ihist = 0
u; ihist < hists.size(); ihist++) {
63 if(hists[ihist] == NULL)
continue;
66 hists[ihist]->SetLineColor(colors[ihist]);
67 hists[ihist]->SetTitle(title.c_str());
69 leg->AddEntry(hists[ihist], labels[ihist].c_str(),
"l");
74 hists[ihist]->GetYaxis()->SetRangeUser(0, max_val * 1.2);
76 hists[ihist]->GetYaxis()->SetRangeUser(0, ymax);
77 hists[ihist]->Draw(draw_option.c_str());
80 hists[ihist]->Draw((draw_option +
" same").c_str());
86 c->Print(name.c_str());
97 TH1 * h_nominal = nominal->
ToTH1(nominal->
POT());
98 PlotDebug(universes, h_nominal, title, name);
107 TDirectory *
tmp = gDirectory;
108 TCanvas *
c =
new TCanvas();
109 c->SetLeftMargin(0.15);
110 c->SetRightMargin(0.01);
113 double max_val = -1e9;
114 for(
auto ihist = 0
u; ihist < universes.size(); ihist++) {
117 h_nominal->GetYaxis()->SetRangeUser(0, max_val * 1.2);
118 h_nominal->SetTitle(title.c_str());
119 h_nominal->Draw(
"hist");
121 for(
auto ihist = 0
u; ihist < universes.size(); ihist++) {
122 universes[ihist]->Draw(
"hist same");
125 h_nominal->Draw(
"hist same");
129 c->Print(name.c_str());
136 const Cut * selection_cut)
180 for(
auto isyst = 0
u; isyst <
fSystDefs.size(); isyst++) {
181 for(
auto iload = 0
u; iload <
fSystDefs[isyst]->LoadersUp().size(); iload++)
182 fSystDefs[isyst]->LoadersUp()[iload]->SetSpillCut(spillcut);
184 for(
auto iload = 0
u; iload <
fSystDefs[isyst]->LoadersDown().size(); iload++)
185 fSystDefs[isyst]->LoadersDown()[iload]->SetSpillCut(spillcut);
195 std::vector<Var> mv_weights)
234 std::vector<SystShifts> mv_shifts) {
268 for(
auto isyst = 0
u; isyst <
fSystDefs.size(); isyst++)
270 for(
auto isyst = 0
u; isyst <
fMVSystDefs.size(); isyst++)
279 TDirectory *
tmp = gDirectory;
281 dir = dir->mkdir(name.c_str());
284 TObjString(
"CutOptimization").Write(
"type");
286 TDirectory * nom_dir = dir->mkdir(
"Nominal");
298 nsysts.Write(
"nsysts");
303 for(
auto isyst = 0
u; isyst <
fSystDefs.size(); isyst++) {
318 nsysts.Write(
"nmultiverse");
342 std::unique_ptr<CutOptimization>
345 dir = dir->GetDirectory(name.c_str());
354 std::vector<SystematicDef *> syst_defs;
355 std::map<SystematicDef *, UpDownPair<Spectrum> > syst_signal;
356 std::map<SystematicDef *, UpDownPair<Spectrum> > syst_bkgd;
357 std::map<SystematicDef *, UpDownPair<Spectrum> > syst_selected_signal;
358 std::map<SystematicDef *, UpDownPair<Spectrum> > syst_selected_bkgd;
359 std::map<SystematicDef *, UpDownPair<Spectrum> > syst_selected;
361 std::vector<SystematicDef *> mv_defs;
362 std::map<SystematicDef *, Multiverse * > mv_signal;
363 std::map<SystematicDef *, Multiverse * > mv_bkgd;
364 std::map<SystematicDef *, Multiverse * > mv_selected_signal;
365 std::map<SystematicDef *, Multiverse * > mv_selected_bkgd;
366 std::map<SystematicDef *, Multiverse * > mv_selected;
368 TDirectory * nom_dir = dir->GetDirectory(
"Nominal");
376 int nsysts = (*vsysts)[0];
377 for(
int isyst = 0; isyst <
nsysts; isyst++) {
380 TDirectory * syst_dir = dir->GetDirectory(
TString::Format(
"syst_%s", syst_defs.back()->GetName().c_str()));
383 syst_selected_signal[syst_defs.back()] =
LoadFromUpDownSpectra(syst_dir->GetDirectory(
"fSelectedSignal" ));
384 syst_selected_bkgd [syst_defs.back()] =
LoadFromUpDownSpectra(syst_dir->GetDirectory(
"fSelectedBkgd" ));
389 int nmultiverse = (*vmultiverse)[0];
390 for(
int imv = 0; imv < nmultiverse; imv++) {
393 TDirectory * mv_dir = dir->GetDirectory(
TString::Format(
"mv_%s", mv_defs.back()->GetName().c_str()));
396 mv_selected_signal[mv_defs.back()] =
Multiverse::LoadFrom(mv_dir,
"fMVSelectedSignal" ).release();
403 return std::make_unique<CutOptimization> (nominal_signal,
405 nominal_selected_signal,
406 nominal_selected_bkgd,
411 syst_selected_signal,
428 std::vector<SystematicDef*> syst_defs,
434 std::vector<SystematicDef*> mv_defs,
435 std::map<SystematicDef*, Multiverse * > mv_signal,
436 std::map<SystematicDef*, Multiverse * > mv_bkgd,
437 std::map<SystematicDef*, Multiverse * > mv_selected_signal,
438 std::map<SystematicDef*, Multiverse * > mv_selected_bkgd,
439 std::map<SystematicDef*, Multiverse * > mv_selected)
472 for(
auto isyst = 0
u; isyst < rhs.
fSystDefs.size(); isyst++) {
480 for(
auto imv = 0
u; imv < rhs.
fMVSystDefs.size(); imv++) {
505 for(
auto isyst = 0
u; isyst < rhs.fSystDefs.size(); isyst++) {
506 this->
fSystDefs.push_back(std::move(rhs.fSystDefs[isyst]));
513 for(
auto imv = 0
u; imv < rhs.fMVSystDefs.size(); imv++) {
514 this->
fMVSystDefs.push_back(std::move(rhs.fMVSystDefs[imv]));
515 this->
fMVSignal [this->
fMVSystDefs.back()] = std::move(rhs.fMVSignal .at(rhs.fMVSystDefs[imv]));
516 this->
fMVBkgd [this->
fMVSystDefs.back()] = std::move(rhs.fMVBkgd .at(rhs.fMVSystDefs[imv]));
528 if(
this == &rhs)
return *
this;
540 for(
auto isyst = 0
u; isyst < rhs.
fSystDefs.size(); isyst++) {
548 for(
auto imv = 0
u; imv < rhs.
fMVSystDefs.size(); imv++) {
574 for(
auto isyst = 0
u; isyst < rhs.fSystDefs.size(); isyst++) {
575 this->
fSystDefs.push_back(std::move(rhs.fSystDefs[isyst]));
582 for(
auto imv = 0
u; imv < rhs.fMVSystDefs.size(); imv++) {
583 this->
fMVSystDefs.push_back(std::move(rhs.fMVSystDefs[imv]));
584 this->
fMVSignal [this->
fMVSystDefs.back()] = std::move(rhs.fMVSignal .at(rhs.fMVSystDefs[imv]));
585 this->
fMVBkgd [this->
fMVSystDefs.back()] = std::move(rhs.fMVBkgd .at(rhs.fMVSystDefs[imv]));
597 auto ndims = tmp->GetDimension();
598 auto NX = tmp->GetNbinsX();
599 auto NY = tmp->GetNbinsY();
600 auto NZ = tmp->GetNbinsZ();
610 assert(ndims == 1 &&
"I don't know how to handle this histogram for kIntegratedInsideBounds");
616 auto array = tmp->GetXaxis()->GetXbins()->GetArray();
619 NX, tmp->GetXaxis()->GetXbins()->GetArray(),
620 NX, tmp->GetXaxis()->GetXbins()->GetArray());
623 auto minx = tmp->GetXaxis()->GetBinLowEdge(1);
624 auto maxx = tmp->GetXaxis()->GetBinLowEdge(NX+1);
630 result->GetYaxis()->SetTitle(
TString::Format(
"%s %s", tmp->GetXaxis()->GetTitle(),
"(Upper)"));
631 result->GetXaxis()->SetTitle(
TString::Format(
"%s %s", tmp->GetXaxis()->GetTitle(),
"(Lower)"));
633 for(
auto ibinlower = 1; ibinlower <= NX; ibinlower++) {
634 for(
auto ibinupper = ibinlower; ibinupper <=NX; ibinupper++) {
636 double integral = tmp->IntegralAndError(ibinlower, ibinupper, error);
637 result->SetBinContent(ibinlower, ibinupper, integral);
638 result->SetBinContent(ibinlower, ibinupper, integral);
647 for(
auto ibinx = 1; ibinx <= NX; ibinx++) {
651 double integral = tmp->IntegralAndError(x_a, x_b, error);
652 result->SetBinContent(ibinx, integral);
653 result->SetBinError (ibinx, error );
661 TH2 * result_2d =
static_cast<TH2*
> (
result);
662 const TH2 * tmp_2d =
static_cast<const TH2*
>(
tmp );
663 for(
auto ibinx = 1; ibinx <= NX; ibinx++) {
664 for(
auto ibiny = 1; ibiny <= NY; ibiny++) {
671 double integral = tmp_2d->IntegralAndError(x_a, x_b,
673 result_2d->SetBinContent(ibinx, ibiny, integral);
674 result_2d->SetBinError (ibinx, ibiny, error );
677 result =
static_cast<TH1*
>(result_2d);
683 TH3 * result_3d =
static_cast<TH3*
> (
result);
684 const TH3 * tmp_3d =
static_cast<const TH3*
>(
tmp );
685 for(
auto ibinx = 1; ibinx <= NX; ibinx++) {
686 for(
auto ibiny = 1; ibiny <= NY; ibiny++) {
687 for(
auto ibinz = 1; ibinz <= NZ; ibinz++) {
696 double integral = tmp_3d->IntegralAndError(x_a, x_b,
699 result_3d->SetBinContent(ibinx, ibiny, ibinz, integral);
700 result_3d->SetBinError (ibinx, ibiny, ibinz, error );
704 result =
static_cast<TH1*
>(result_3d);
721 TH1 *
ret = (TH1*) h_up->Clone();
722 for(
auto ibinx = 1; ibinx <= ret->GetNbinsX(); ibinx++) {
723 for(
auto ibiny = 1; ibiny <= ret->GetNbinsY(); ibiny++) {
724 for(
auto ibinz = 1; ibinz <= ret->GetNbinsZ(); ibinz++) {
725 double up = h_up ->GetBinContent(ibinx, ibiny, ibinz);
726 double down = h_down ->GetBinContent(ibinx, ibiny, ibinz);
727 double nominal = h_nominal ->GetBinContent(ibinx, ibiny, ibinz);
729 ret->SetBinContent(ibinx, ibiny, ibinz, abs_shift + nominal);
745 if(spec == NULL)
return NULL;
748 if(POT < 0) POT = spec->
POT();
754 tmp = spec->
ToTH1(POT);
757 tmp = spec->
ToTH2(POT);
760 tmp = spec->
ToTH3(POT);
763 std::cout <<
"Spectrum appears to be empty. Something went terribly wrong here" <<
std::endl;
768 ret = (TH1*) tmp->Clone();
798 auto NX = denom->GetNbinsX();
799 auto NY = denom->GetNbinsY();
800 auto NZ = denom->GetNbinsZ();
803 integral = method ==
kIntegratedUp ? denom->GetBinContent(1) : denom->GetBinContent(NX, NY, NZ);
810 integral = denom->GetBinContent(1, NY, 1);
813 eff->Scale(1/integral);
815 for(
auto ibinx = 1; ibinx <= NX; ibinx++) {
816 for(
auto ibiny = 1; ibiny <= NY; ibiny++) {
817 for(
auto ibinz = 1; ibinz <= NZ; ibinz++) {
818 double efficiency = eff->GetBinContent(ibinx, ibiny, ibinz);
819 double binom_error = 1/integral *
std::sqrt( efficiency * (1 - efficiency / integral) );
820 eff->SetBinError(ibinx, ibiny, ibinz, binom_error);
827 eff->Divide(eff, denom, 1, 1,
"B");
830 eff->SetMaximum(1.2*eff->GetMaximum());
854 [method](TH1 * univ_num, TH1 * univ_denom)
860 eff->
Divide(*denom,
true);
871 TH1 *
result = (TH1*) univ->Clone();
892 TH1 * hNominalEfficiency =
Efficiency(hNominalSelectedSignal, hNominalSignal, method);
896 std::map<SystematicDef*, Multiverse *> fMVEfficiency;
900 fMVBkgd [*mv]->Scale(mc_pot / mc_pot);
916 std::vector<UpDownPair<TH1 > > hSystSignal;
917 std::vector<UpDownPair<TH1 > > hSystBkgd;
918 std::vector<UpDownPair<TH1 > > hSystSelectedSignal;
919 std::vector<UpDownPair<TH1 > > hSystSelectedBkgd;
920 std::vector<UpDownPair<TH1 > > hSystSelected;
921 std::vector<UpDownPair<TH1 > > hSystEfficiency;
923 std::vector<SystematicDef*>
systs;
924 for(
auto isyst = 0
u; isyst <
fSystDefs.size(); isyst++) {
926 systs.push_back(syst);
930 hSystBkgd .push_back({
ToHist(
fBkgd .at(syst).up, method, mc_pot ),
ToHist(
fBkgd .at(syst).down, method, mc_pot )});
933 hSystEfficiency .push_back(
Efficiency(hSystSelectedSignal.back(), hSystSignal.back(), method));
945 hSystEfficiency .push_back(
ToUpDownHist(fMVEfficiency [mv], hNominalEfficiency ));
948 TH1 * hFracUncertaintySelected = (TH1*) hNominalSelected ->
Clone(
UniqueName().c_str());
949 TH1 * hFracUncertaintyBkgdStat = (TH1*) hNominalSelectedBkgd ->
Clone(
UniqueName().c_str());
950 TH1 * hFracUncertaintyBkgdSyst = (TH1*) hNominalBkgd ->
Clone(
UniqueName().c_str());
951 TH1 * hFracUncertaintyEfficiency = (TH1*) hNominalSignal ->
Clone(
UniqueName().c_str());
952 TH1 * hFracUncertaintyXSec = (TH1*) hNominalSignal ->
Clone(
UniqueName().c_str());
954 hFracUncertaintySelected ->Scale(0);
955 hFracUncertaintyBkgdStat ->Scale(0);
956 hFracUncertaintyBkgdSyst ->Scale(0);
957 hFracUncertaintyEfficiency ->Scale(0);
958 hFracUncertaintyXSec ->Scale(0);
960 std::vector<TH1*> vabs_uncert_eff_syst (systs.size());
961 std::vector<TH1*> vabs_uncert_bkgd_syst (systs.size());
962 std::vector<TH1*> vfrac_uncert_eff_syst (systs.size());
963 std::vector<TH1*> vfrac_uncert_bkgd_syst(systs.size());
964 TH1* sel_minus_bkgd = (TH1*) hNominalSelected->Clone(
UniqueName().c_str());
965 sel_minus_bkgd->Scale(0);
966 for(
auto isyst = 0
u; isyst < systs.size(); isyst++) {
967 vabs_uncert_eff_syst [isyst] = (TH1*) hNominalSelected->Clone(
UniqueName().c_str());
968 vabs_uncert_bkgd_syst [isyst] = (TH1*) hNominalSelected->Clone(
UniqueName().c_str());
969 vfrac_uncert_eff_syst [isyst] = (TH1*) hNominalSelected->Clone(
UniqueName().c_str());
970 vfrac_uncert_bkgd_syst[isyst] = (TH1*) hNominalSelected->Clone(
UniqueName().c_str());
972 if(hNominalSelected->GetDimension() == 1) {
973 vfrac_uncert_bkgd_syst[isyst] ->GetYaxis()->SetTitle(
"#frac{s#upoint#delta N^{syst}_{bkgd}}{N_{sel} - s#upoint N_{bkgd}}");
974 vfrac_uncert_eff_syst [isyst] ->GetYaxis()->SetTitle(
"#frac{#delta#varepsilon}{#varepsilon}");
975 vabs_uncert_bkgd_syst [isyst] ->GetYaxis()->SetTitle(
"s#upoint#delta N^{syst}_{bkgd}");
976 vabs_uncert_eff_syst [isyst] ->GetYaxis()->SetTitle(
"#delta#varepsilon");
979 else if(hNominalSelected->GetDimension() == 2) {
980 vfrac_uncert_bkgd_syst[isyst] ->GetZaxis()->SetTitle(
"#frac{s#upoint#delta N^{syst}_{bkgd}}{N_{sel} - s#upoint N_{bkgd}}");
981 vfrac_uncert_eff_syst [isyst] ->GetZaxis()->SetTitle(
"#frac{#delta#varepsilon}{#varepsilon}");
982 vabs_uncert_bkgd_syst [isyst] ->GetZaxis()->SetTitle(
"s#upoint#delta N^{syst}_{bkgd}");
983 vabs_uncert_eff_syst [isyst] ->GetZaxis()->SetTitle(
"#delta#varepsilon");
985 vabs_uncert_eff_syst [isyst] ->Scale(0);
986 vabs_uncert_bkgd_syst [isyst] ->Scale(0);
987 vfrac_uncert_eff_syst [isyst] ->Scale(0);
988 vfrac_uncert_bkgd_syst[isyst] ->Scale(0);
991 auto NX = hNominalSignal->GetNbinsX();
992 auto NY = hNominalSignal->GetNbinsY();
993 auto NZ = hNominalSignal->GetNbinsZ();
994 for(
auto ibinx = 1; ibinx <= NX; ibinx++) {
995 for(
auto ibiny = 1; ibiny <= NY; ibiny++) {
996 for(
auto ibinz = 1; ibinz <= NZ; ibinz++) {
997 double abs_uncert_bkgd_syst = 0;
998 double abs_uncert_eff = 0;
1000 double nom_sel = hNominalSelected ->GetBinContent(ibinx, ibiny, ibinz);
1001 double nom_bkgd = hNominalSelectedBkgd ->GetBinContent(ibinx, ibiny, ibinz);
1002 double nom_eff = hNominalEfficiency ->GetBinContent(ibinx, ibiny, ibinz);
1003 double nom_signal = nom_sel - s * nom_bkgd;
1005 sel_minus_bkgd->SetBinContent(ibinx, ibiny, ibinz, nom_signal);
1006 for(
auto isyst = 0
u; isyst < systs.size(); isyst++) {
1007 double syst_bkgd = 0;
1008 double syst_eff = 0;
1010 if(systs[isyst]->IsTwoSided()) {
1012 double up_bkgd = hSystSelectedBkgd[isyst].up ->GetBinContent(ibinx, ibiny, ibinz);
1013 double down_bkgd = hSystSelectedBkgd[isyst].down->GetBinContent(ibinx, ibiny, ibinz);
1015 double up_eff = hSystEfficiency [isyst].up ->GetBinContent(ibinx, ibiny, ibinz);
1016 double down_eff = hSystEfficiency [isyst].down->GetBinContent(ibinx, ibiny, ibinz);
1022 syst_bkgd =
std::abs(hSystSelectedBkgd[isyst].up ->GetBinContent(ibinx, ibiny, ibinz) - nom_bkgd);
1023 syst_eff =
std::abs(hSystEfficiency [isyst].up ->GetBinContent(ibinx, ibiny, ibinz) - nom_eff );
1027 syst_bkgd = s * syst_bkgd;
1029 abs_uncert_bkgd_syst +=
std::pow(syst_bkgd, 2);
1030 abs_uncert_eff +=
std::pow(syst_eff , 2);
1032 vabs_uncert_eff_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_eff );
1033 vabs_uncert_bkgd_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_bkgd);
1035 vfrac_uncert_bkgd_syst[isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_bkgd / nom_signal);
1037 vfrac_uncert_bkgd_syst[isyst]->SetBinContent(ibinx, ibiny, ibinz, 100);
1039 vfrac_uncert_eff_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_eff / nom_eff);
1041 vfrac_uncert_eff_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, 100);
1045 double frac_uncert_sel = 100;
1046 double frac_uncert_bkgd_stat = 100;
1047 double frac_uncert_bkgd_syst = 100;
1048 double frac_uncert_eff = 100;
1050 frac_uncert_sel =
std::sqrt(nom_sel) / nom_signal;
1051 frac_uncert_bkgd_stat = s *
std::sqrt(nom_bkgd) / nom_signal;
1052 frac_uncert_bkgd_syst =
std::sqrt(abs_uncert_bkgd_syst) / nom_signal;
1055 frac_uncert_eff =
std::sqrt(abs_uncert_eff) / nom_eff;
1059 std::pow(frac_uncert_bkgd_stat, 2) +
1060 std::pow(frac_uncert_bkgd_syst, 2) +
1065 hFracUncertaintyXSec ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_sigma );
1066 hFracUncertaintySelected ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_sel );
1067 hFracUncertaintyBkgdStat ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_bkgd_stat);
1068 hFracUncertaintyBkgdSyst ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_bkgd_syst);
1069 hFracUncertaintyEfficiency ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_eff );
1078 TFile * fdebug =
new TFile((plot_dump +
"/debug.root").c_str(),
"recreate");
1080 hNominalSelectedBkgd->Write(
"hNominalSelectedBkgd");
1081 hNominalSelected->Write(
"hNominalSelected");
1082 hNominalSelectedSignal->Write(
"hNominalSelectedSignal");
1083 hNominalSignal->Write(
"hNominalSignal");
1084 hNominalBkgd->Write(
"hNominalBkgd");
1085 hNominalEfficiency->Write(
"hNominalEfficiency");
1087 hFracUncertaintyXSec ->Write(
"hFracUncertaintyXSec" );
1088 hFracUncertaintySelected ->Write(
"hFracUncertaintySelected" );
1089 hFracUncertaintyBkgdStat ->Write(
"hFracUncertaintyBkgdStat" );
1090 hFracUncertaintyBkgdSyst ->Write(
"hFracUncertaintyBkgdSyst" );
1091 hFracUncertaintyEfficiency ->Write(
"hFracUncertaintyEfficiency" );
1094 for(
auto isyst = 0
u; isyst < systs.size(); isyst++ ) {
1095 hSystSelectedBkgd [isyst].up->Write((systs[isyst]->
GetName() +
"_hSystSelectedBkgd" +
"_UP").c_str());
1096 hSystSelected [isyst].up->Write((systs[isyst]->
GetName() +
"_hSystSelected" +
"_UP").c_str());
1097 hSystSelectedSignal[isyst].up->Write((systs[isyst]->
GetName() +
"_hSystSelectedSignal" +
"_UP").c_str());
1098 hSystBkgd [isyst].up->Write((systs[isyst]->
GetName() +
"_hSystBkgd" +
"_UP").c_str());
1099 hSystSignal [isyst].up->Write((systs[isyst]->
GetName() +
"_hSystSignal" +
"_UP").c_str());
1100 hSystEfficiency [isyst].up->Write((systs[isyst]->
GetName() +
"_hSystEfficiency" +
"_UP").c_str());
1102 if(hSystSelectedBkgd[isyst].down != NULL) {
1103 hSystSelectedBkgd [isyst].down->Write((systs[isyst]->
GetName() +
"_hSystSelectedBkgd" +
"_DOWN").c_str());
1104 hSystSelected [isyst].down->Write((systs[isyst]->
GetName() +
"_hSystSelected" +
"_DOWN").c_str());
1105 hSystSelectedSignal[isyst].down->Write((systs[isyst]->
GetName() +
"_hSystSelectedSignal" +
"_DOWN").c_str());
1106 hSystBkgd [isyst].down->Write((systs[isyst]->
GetName() +
"_hSystBkgd" +
"_DOWN").c_str());
1107 hSystSignal [isyst].down->Write((systs[isyst]->
GetName() +
"_hSystSignal" +
"_DOWN").c_str());
1108 hSystEfficiency [isyst].down->Write((systs[isyst]->
GetName() +
"_hSystEfficiency" +
"_DOWN").c_str());
1111 vabs_uncert_eff_syst [isyst]->Write((systs[isyst]->
GetName() +
"_abs_uncert_eff_syst" ).c_str());
1112 vabs_uncert_bkgd_syst [isyst]->Write((systs[isyst]->
GetName() +
"_abs_uncert_bkgd_syst" ).c_str());
1113 vfrac_uncert_eff_syst [isyst]->Write((systs[isyst]->
GetName() +
"_frac_uncert_eff_syst" ).c_str());
1114 vfrac_uncert_bkgd_syst[isyst]->Write((systs[isyst]->
GetName() +
"_frac_uncert_bkgd_syst").c_str());
1160 fMVSignal [syst]->SaveTo(fdebug,
"MV_" + syst->GetName() +
"_Signal" );
1161 fMVBkgd [syst]->SaveTo(fdebug,
"MV_" + syst->GetName() +
"_Bkgd" );
1162 fMVSelected [syst]->SaveTo(fdebug,
"MV_" + syst->GetName() +
"_Selected" );
1163 fMVSelectedSignal[syst]->SaveTo(fdebug,
"MV_" + syst->GetName() +
"_SelectedSignal");
1164 fMVSelectedBkgd [syst]->SaveTo(fdebug,
"MV_" + syst->GetName() +
"_SelectedBkgd" );
1165 fMVEfficiency [syst]->SaveTo(fdebug,
"MV_" + syst->GetName() +
"_Efficiency" );
1198 hFracUncertaintyXSec ->SetTitle(
"Fractional Uncertainty on Cross-section" );
1199 hFracUncertaintySelected ->SetTitle(
"Fractional Stat. Uncertainty on Selection" );
1200 hFracUncertaintyBkgdStat ->SetTitle(
"Fractional Stat. Uncertainty on Background");
1201 hFracUncertaintyBkgdSyst ->SetTitle(
"Fractional Syst. Uncertainty on Background");
1202 hFracUncertaintyEfficiency ->SetTitle(
"Fractional Syst. Uncertainty on Efficiency");
1206 hFracUncertaintyXSec ->GetYaxis()->SetRangeUser(0, maxy);
1207 hFracUncertaintySelected ->GetYaxis()->SetRangeUser(0, maxy);
1208 hFracUncertaintyBkgdStat ->GetYaxis()->SetRangeUser(0, maxy);
1209 hFracUncertaintyBkgdSyst ->GetYaxis()->SetRangeUser(0, maxy);
1210 hFracUncertaintyEfficiency ->GetYaxis()->SetRangeUser(0, maxy);
1212 hFracUncertaintyXSec ->GetYaxis()->SetTitle(
"#frac{#delta#sigma}{#sigma}" );
1213 hFracUncertaintySelected ->GetYaxis()->SetTitle(
"#frac{#sqrt{N_{sel}}}{N_{sel}-s#upoint N_{bkg}}" );
1214 hFracUncertaintyBkgdStat ->GetYaxis()->SetTitle(
"#frac{s#upoint#sqrt{N_{bkg}}}{N_{sel}-s#upoint N_{bkg}}");
1215 hFracUncertaintyBkgdSyst ->GetYaxis()->SetTitle(
"#frac{#sqrt{s#upoint#deltaN^{syst}_{bkg}}}{N_{sel}-s#upoint N_{bkg}}");
1216 hFracUncertaintyEfficiency ->GetYaxis()->SetTitle(
"#frac{#delta#varepsilon}{#varepsilon}" );
1218 hFracUncertaintyXSec ->GetYaxis()->CenterTitle();
1219 hFracUncertaintySelected ->GetYaxis()->CenterTitle();
1220 hFracUncertaintyBkgdStat ->GetYaxis()->CenterTitle();
1221 hFracUncertaintyBkgdSyst ->GetYaxis()->CenterTitle();
1222 hFracUncertaintyEfficiency ->GetYaxis()->CenterTitle();
1224 hFracUncertaintyXSec ->SetLineColor(
kRed);
1225 hFracUncertaintySelected ->SetLineColor(
kRed);
1226 hFracUncertaintyBkgdStat ->SetLineColor(
kRed);
1227 hFracUncertaintyBkgdSyst ->SetLineColor(
kRed);
1228 hFracUncertaintyEfficiency ->SetLineColor(
kRed);
1231 TCanvas *
canvas =
new TCanvas();
1232 canvas->SetLeftMargin(0.15);
1233 canvas->SetRightMargin(0.01);
1236 hFracUncertaintyXSec ->Draw(
"hist");
1240 hFracUncertaintySelected ->Draw(
"hist");
1244 hFracUncertaintyBkgdStat ->Draw(
"hist");
1245 canvas ->Print(
TString::Format(
"%s/optimize_%s_frac_uncert_bkgd_stat.pdf",
1248 hFracUncertaintyBkgdSyst ->Draw(
"hist");
1249 canvas ->Print(
TString::Format(
"%s/optimize_%s_frac_uncert_bkgd_syst.pdf",
1252 hFracUncertaintyEfficiency ->Draw(
"hist");
1260 hFracUncertaintyXSec ->GetZaxis()->SetRangeUser(0, maxy);
1261 hFracUncertaintySelected ->GetZaxis()->SetRangeUser(0, maxy);
1262 hFracUncertaintyBkgdStat ->GetZaxis()->SetRangeUser(0, maxy);
1263 hFracUncertaintyBkgdSyst ->GetZaxis()->SetRangeUser(0, maxy);
1264 hFracUncertaintyEfficiency ->GetZaxis()->SetRangeUser(0, maxy);
1266 hFracUncertaintyXSec ->GetZaxis()->SetTitle(
"#frac{#delta#sigma}{#sigma}" );
1267 hFracUncertaintySelected ->GetZaxis()->SetTitle(
"#frac{#sqrt{N_{sel}}}{N_{sel}-N_{bkg}}" );
1268 hFracUncertaintyBkgdStat ->GetZaxis()->SetTitle(
"#frac{#sqrt{N_{bkg}}}{N_{sel}-N_{bkg}}" );
1269 hFracUncertaintyBkgdSyst ->GetZaxis()->SetTitle(
"#frac{#deltaN^{syst}_{bkg}}{N_{sel}-N_{bkg}}");
1270 hFracUncertaintyEfficiency ->GetZaxis()->SetTitle(
"#frac{#delta#varepsilon}{#varepsilon}" );
1272 hFracUncertaintyXSec ->GetYaxis()->CenterTitle();
1273 hFracUncertaintySelected ->GetYaxis()->CenterTitle();
1274 hFracUncertaintyBkgdStat ->GetYaxis()->CenterTitle();
1275 hFracUncertaintyBkgdSyst ->GetYaxis()->CenterTitle();
1276 hFracUncertaintyEfficiency ->GetYaxis()->CenterTitle();
1279 TCanvas *
canvas =
new TCanvas();
1283 hFracUncertaintyXSec ->Draw(
"colz");
1284 canvas ->Print(
TString::Format(
"%s/optimize_%s_%s_frac_uncert_xsec.pdf",
1288 hFracUncertaintySelected ->Draw(
"colz");
1289 canvas ->Print(
TString::Format(
"%s/optimize_%s_%s_frac_uncert_sel.pdf",
1293 hFracUncertaintyBkgdStat ->Draw(
"colz");
1294 canvas ->Print(
TString::Format(
"%s/optimize_%s_%s_frac_uncert_bkgd_stat.pdf",
1298 hFracUncertaintyBkgdSyst ->Draw(
"colz");
1299 canvas ->Print(
TString::Format(
"%s/optimize_%s_%s_frac_uncert_bkgd_syst.pdf",
1303 hFracUncertaintyEfficiency ->Draw(
"colz");
1304 canvas ->Print(
TString::Format(
"%s/optimize_%s_%s_frac_uncert_eff.pdf",
1316 hFracUncertaintyXSec ->Write(
"frac_uncert_xsec" );
1317 hFracUncertaintySelected ->Write(
"frac_uncert_sel" );
1318 hFracUncertaintyBkgdStat ->Write(
"frac_uncert_bkgd_stat");
1319 hFracUncertaintyBkgdSyst ->Write(
"frac_uncert_bkgd_syst");
1320 hFracUncertaintyEfficiency ->Write(
"frac_uncert_eff" );
1327 TH1 *
ret = (TH1*) syst.
up->Clone();
1328 int NX = syst.
up->GetNbinsX();
1329 int NY = syst.
up->GetNbinsY();
1330 int NZ = syst.
up->GetNbinsZ();
1331 for(
auto ibinx = 1; ibinx <= NX; ibinx++) {
1332 for(
auto ibiny = 1; ibiny <= NY; ibiny++) {
1333 for(
auto ibinz = 1; ibinz <= NZ; ibinz++) {
1334 double nom = nominal->GetBinContent(ibinx, ibiny, ibinz);
1336 double up = syst.
up->GetBinContent(ibinx, ibiny, ibinz);
1338 if(syst.
down) down = syst.
down->GetBinContent(ibinx, ibiny, ibinz);
1341 ret->SetBinContent(ibinx, ibiny, ibinz, error);
T max(const caf::Proxy< T > &a, T b)
UpDownPair< TH1 > ToUpDownHist(Multiverse *mv, const TH1 *h_nominal)
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
_HistAxis< Var > HistAxis
Cuts and Vars for the 2020 FD DiF Study.
const SystShifts * Shifts() const
Spectrum * GetPlusOneSigmaShift(const Spectrum *)
static std::unique_ptr< GenericSystematicDef< SRType > > LoadFrom(TDirectory *dir, const std::string &name)
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.
std::map< SystematicDef *, UpDownPair< Spectrum > > fSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedBkgd
std::vector< SpectrumLoaderBase * > Loaders() const
Spectrum * BuildSpectrumDown(const _Cut< SRType > cut)
const _Var< SRType > * Weight() const
CutOptimization(const HistAxis *hist_axis, const Cut *signal_cut, const Cut *selection_cut)
const HistAxis * fOptiHistAxis
static std::unique_ptr< CutOptimization > LoadFrom(TDirectory *dir, const std::string &name)
std::vector< double > Spectrum
std::map< SystematicDef *, Multiverse * > fMVBkgd
static TH1 * ToHist(const Spectrum *spec, OptimizeMethod_t method, double POT=-1)
std::map< SystematicDef *, UpDownPair< Spectrum > > fBkgd
TCanvas * canvas(const char *nm, const char *ti, const char *a)
GenericSystematicDef< caf::SRProxy > SystematicDef
Spectrum * GetMinusOneSigmaShift(const Spectrum *)
void DefineNominal(SystematicDef *nominal)
UpDownPair< Spectrum > CopyUpDownSpectrum(const UpDownPair< Spectrum > ©)
Representation of a spectrum in any variable, with associated POT.
CutOptimization & operator=(const CutOptimization &rhs)
static void PlotDebug(TH1 *hist, std::string draw_option, std::string title, std::string name)
std::string GetName(int i)
Spectrum * BuildSpectrumUp(const _Cut< SRType > cut)
Spectrum * fNominalSignal
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
std::map< SystematicDef *, Multiverse * > fMVSelectedBkgd
void SaveTo(TDirectory *dir, const std::string &name) const
void DefineMultiverseSystematic(SystematicDef *syst, std::vector< Var > mv_weights)
static TH1 * Efficiency(TH1 *num, TH1 *denom, OptimizeMethod_t method)
const Cut * fSelectionCut
Spectrum * fNominalSelectedSignal
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
void SetSpillCuts(const SpillCut &)
void SaveTo(TDirectory *dir, const std::string &name) const
unsigned int NDimensions() const
static void Integrate(TH1 *&result, const TH1 *tmp, OptimizeMethod_t method)
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelectedSignal
static std::unique_ptr< Multiverse > LoadFrom(TDirectory *dir, const std::string &name)
std::vector< double > POT
void Optimize(OptimizeMethod_t method, FOM_t fom, double data_pot, TDirectory *save_dir, std::string plot_dump=".", bool debug=false)
void OptimizedSigmaOverSigma(OptimizeMethod_t method, double data_pot, TDirectory *save_dir, std::string plot_dump=".", bool debug=false)
Spectrum * fNominalSelected
TH3 * ToTH3(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 3D to obtain TH3.
void SaveToUpDownSpectra(const UpDownPair< Spectrum > &save, TDirectory *dir)
void DefineSystematic(SystematicDef *syst)
std::map< SystematicDef *, Multiverse * > fMVSelectedSignal
std::map< SystematicDef *, UpDownPair< Spectrum > > fSelected
assert(nhit_max >=nhit_nbins)
std::vector< SystematicDef * > fSystDefs
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
UpDownPair< Spectrum > LoadFromUpDownSpectra(TDirectory *dir)
TH1 * AbsUncertaintySquared(const UpDownPair< TH1 > syst, const TH1 *nominal) const
Spectrum * BuildSpectrum(const _Cut< SRType > cut)
const std::vector< std::string > & GetLabels() const
std::map< SystematicDef *, Multiverse * > fMVSelected
void Transform(const Multiverse &rhs, const std::function< ProductFunc_t > &transform)
std::string UniqueName()
Return a different string each time, for creating histograms.
std::vector< SystematicDef * > fMVSystDefs
std::vector< SpectrumLoaderBase * > LoadersUp() const
SystematicDef * fNominalSystDef
void Divide(Multiverse &, bool=false)
Spectrum * fNominalSelectedBkgd
std::map< SystematicDef *, Multiverse * > fMVSignal
TH1 * AbsUncertainty(const UpDownPair< TH1 > syst, const TH1 *nominal) const