8 #include "Utilities/func/MathUtil.h" 11 #include "TFeldmanCousins.h" 13 #include "TGraphAsymmErrors.h" 17 #include "THLimitsFinder.h" 22 #include "TPaveText.h" 41 TH1* hData = data.
ToTH1(data.
POT(), kBlack, kSolid,
kPOT, bintype);
43 hData->SetMarkerStyle(kFullCircle);
46 if(hMC->GetMaximum() > hData->GetMaximum()+
sqrt(hData->GetMaximum())){
49 hData->Draw(
"ep same");
54 hMC->Draw(
"hist same");
55 hData->Draw(
"ep same");
83 hData->SetMarkerStyle(kFullCircle);
85 hMC->Scale(hData->Integral()/hMC->Integral());
88 if(hMC->GetMaximum() > hData->GetMaximum()+
sqrt(hData->GetMaximum())){
91 hData->Draw(
"ep same");
96 hMC->Draw(
"hist same");
97 hData->Draw(
"ep same");
141 TH1* hData = data.
ToTH1(data.
POT());
142 hData->SetMarkerStyle(kFullCircle);
145 if(hMC->GetMaximum() > hData->GetMaximum()+
sqrt(hData->GetMaximum())){
148 hData->Draw(
"ep same");
153 hMC->Draw(
"hist same");
154 hData->Draw(
"ep same");
157 hNC->Draw(
"hist same");
158 hCC->Draw(
"hist same");
159 hBeam->Draw(
"hist same");
170 double miny,
double maxy)
178 double miny,
double maxy)
182 h->GetYaxis()->SetTitle(
"Data / MC");
183 h->GetYaxis()->SetRangeUser(miny, maxy);
186 TLine*
one =
new TLine(h->GetXaxis()->GetXmin(), 1,
187 h->GetXaxis()->GetXmax(), 1);
188 one->SetLineWidth(2);
189 one->SetLineStyle(7);
199 double miny,
double maxy)
207 double miny,
double maxy)
219 double miny,
double maxy)
221 Ratio fitRatio(fit, expected);
222 Ratio dataRatio(data, expected);
224 TH1*
h = fitRatio.
ToTH1();
225 h->GetYaxis()->SetTitle(
"Ratio to expectation");
226 h->GetYaxis()->SetRangeUser(miny, maxy);
230 h = dataRatio.
ToTH1();
231 h->SetMarkerStyle(kFullCircle);
234 TLine*
one =
new TLine(h->GetXaxis()->GetXmin(), 1,
235 h->GetXaxis()->GetXmax(), 1);
236 one->SetLineWidth(2);
237 one->SetLineStyle(7);
246 std::vector<const ISyst*>
systs,
247 std::vector<TH1*> &hUps,
248 std::vector<TH1*> &hDns,
254 for (
const ISyst* syst: systs) {
273 std::vector<const ISyst*>
systs,
275 std::vector<TH1*> &hUps,
276 std::vector<TH1*> &hDns,
282 for (
const ISyst* syst: systs) {
305 const std::vector<const ISyst*>&
systs,
308 int col,
int errCol,
float headroom,
316 std::vector<Spectrum> ups, dns;
318 for(
const ISyst* syst: systs){
334 const std::vector<Spectrum>& upShifts,
335 const std::vector<Spectrum>& downShifts,
336 double pot,
int col,
int errCol,
337 float headroom,
bool newaxis,
342 TH1* nom = nominal.
ToTH1(pot, kBlack, kSolid,
kPOT, bintype);
344 std::vector<TH1*> ups, dns;
346 for(
const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot, kBlack, kSolid,
kPOT, bintype));
347 for(
const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot, kBlack, kSolid,
kPOT, bintype));
351 for(TH1* up: ups)
delete up;
352 for(TH1* dn: dns)
delete dn;
360 std::vector<TH1*>& ups,
361 std::vector<TH1*>& dns,
363 float headroom,
bool newaxis,
double alpha)
369 else if(errCol == -1) errCol = col-7;
371 nom->SetLineColor(col);
372 nom->GetXaxis()->CenterTitle();
373 nom->GetYaxis()->CenterTitle();
374 if(newaxis) nom->Draw(
"hist ][");
376 double yMax = nom->GetBinContent(nom->GetMaximumBin());
378 TGraphAsymmErrors*
g =
new TGraphAsymmErrors;
380 for(
int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
381 const double y = nom->GetBinContent(binIdx);
382 g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx),
y);
384 const double w = nom->GetXaxis()->GetBinWidth(binIdx);
389 for(
unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
390 double hi = ups[systIdx]->GetBinContent(binIdx)-
y;
391 double lo = dns[systIdx]->GetBinContent(binIdx)-
y;
403 g->SetPointError(binIdx, w/2, w/2,
sqrt(errDn),
sqrt(errUp));
406 g->SetFillColorAlpha(errCol, alpha);
408 g->GetYaxis()->SetRangeUser(0.00001, headroom*yMax);
409 nom->GetYaxis()->SetRangeUser(0.0001, headroom*yMax);
411 nom->Draw(
"hist same");
420 const std::vector<const ISyst*>&
systs,
423 int col,
int errCol,
float headroom,
429 std::vector<Spectrum> ups, dns;
431 for(
const ISyst* syst: systs){
446 std::vector<Spectrum> upShifts,
447 std::vector<Spectrum> downShifts,
448 double pot,
int col,
int errCol,
449 float headroom,
bool newaxis,
452 double nomIntegral = nominal.
Integral(pot);
456 for (
auto & collection : std::vector<std::reference_wrapper<std::vector<Spectrum>>>{ upShifts, downShifts })
458 for (
Spectrum & shiftSpec : collection.get())
461 shiftSpec.OverridePOT(shiftSpec.POT() * shiftSpec.Integral(pot) / nomIntegral);
471 const std::vector<Spectrum>& upShifts,
472 const std::vector<Spectrum>& downShifts,
474 const std::vector<Spectrum>& upShiftsAlt,
475 const std::vector<Spectrum>& downShiftsAlt,
476 double pot,
int col1,
int col2,
477 float headroom,
bool newaxis,
480 TH1* nom = nominal.
ToTH1(pot, kBlack, kSolid,
kPOT, bintype);
481 TH1* alt = alternative.
ToTH1(pot, kBlack, kSolid,
kPOT, bintype);
482 std::vector<TH1*> ups, dns;
483 std::vector<TH1*> upsA, dnsA;
484 for(
const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot, kBlack, kSolid,
kPOT, bintype));
485 for(
const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot, kBlack, kSolid,
kPOT, bintype));
486 for(
const auto& upShiftAlt:upShiftsAlt) upsA.push_back(upShiftAlt.ToTH1(pot, kBlack, kSolid,
kPOT, bintype));
487 for(
const auto& downShiftAlt:downShiftsAlt) dnsA.push_back(downShiftAlt.ToTH1(pot, kBlack, kSolid,
kPOT, bintype));
488 nom->SetLineColor(col1);
489 nom->GetXaxis()->CenterTitle();
490 nom->GetYaxis()->CenterTitle();
491 if(newaxis) nom->Draw(
"hist");
492 alt->SetLineColor(col2);
493 double yMax = nom->GetBinContent(nom->GetMaximumBin());
494 TGraphAsymmErrors*
g =
new TGraphAsymmErrors;
495 TGraphAsymmErrors* g2 =
new TGraphAsymmErrors;
496 for(
int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
497 const double y = nom->GetBinContent(binIdx);
498 const double y2 = alt->GetBinContent(binIdx);
499 g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx),
y);
500 g2->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx),
y2);
501 const double w = nom->GetXaxis()->GetBinWidth(binIdx);
502 double errUp = 0, errUp2 = 0;
503 double errDn = 0, errDn2 = 0;
504 for(
unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
505 double hi = ups[systIdx]->GetBinContent(binIdx)-
y;
506 double hi2 = upsA[systIdx]->GetBinContent(binIdx)-
y2;
507 double lo = y-dns[systIdx]->GetBinContent(binIdx);
508 double lo2 = y2-dnsA[systIdx]->GetBinContent(binIdx);
509 if(lo <= 0 && hi <= 0)
std::swap(lo, hi);
510 if(lo2 <= 0 && hi2 <= 0)
std::swap(lo2, hi2);
517 g->SetPointError(binIdx, w/2, w/2,
sqrt(errDn),
sqrt(errUp));
518 g2->SetPointError(binIdx, w/2, w/2,
sqrt(errDn2),
sqrt(errUp2));
520 g2->SetFillColorAlpha(col2, 0.3);
523 g2->GetYaxis()->SetRangeUser(0.0, 19.0);
524 g->SetFillColorAlpha(col1, 0.3);
528 g->GetYaxis()->SetRangeUser(0.0, yMax*headroom);
529 nom->GetYaxis()->SetRangeUser(0.0, yMax*headroom);
530 alt->Draw(
"hist same");
531 nom->Draw(
"hist same");
532 for(TH1* up: ups)
delete up;
533 for(TH1* dn: dns)
delete dn;
534 for(TH1* up2: upsA)
delete up2;
535 for(TH1* dn2: dnsA)
delete dn2;
545 const std::vector<const ISyst*>&
systs,
555 std::vector<Spectrum> ups, dns;
557 for(
const ISyst* syst: systs){
573 const std::vector<Spectrum>& upShifts,
574 const std::vector<Spectrum>& downShifts,
575 double pot,
int col,
int errCol,
576 float maxy,
bool newaxis)
579 TH1* nom = nominal.
ToTH1(pot);
580 nom->GetYaxis()->SetTitle(
"Events/0.1 GeV");
581 std::vector<TH1*> ups, dns;
583 for(
const auto& upShift:upShifts) ups.push_back(upShift.ToTH1(pot));
584 for(
const auto& downShift:downShifts) dns.push_back(downShift.ToTH1(pot));
588 for(TH1* up: ups)
delete up;
589 for(TH1* dn: dns)
delete dn;
597 std::vector<TH1*>& ups,
598 std::vector<TH1*>& dns,
600 float maxy,
bool newaxis)
606 else if(errCol == -1) errCol = col-7;
608 nom->GetYaxis()->SetTitle(
"Events/0.1 GeV");
609 nom->SetLineColor(col);
610 nom->GetXaxis()->CenterTitle();
611 nom->GetYaxis()->CenterTitle();
612 if(newaxis) nom->Draw(
"hist ][");
614 TGraphAsymmErrors*
g =
new TGraphAsymmErrors;
616 for(
int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
617 const double y = nom->GetBinContent(binIdx);
618 g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx),
y);
620 const double w = nom->GetXaxis()->GetBinWidth(binIdx);
625 for(
unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
626 double hi = ups[systIdx]->GetBinContent(binIdx)-
y;
627 double lo = dns[systIdx]->GetBinContent(binIdx)-
y;
639 g->SetPointError(binIdx, w/2, w/2,
sqrt(errDn),
sqrt(errUp));
644 if(quant==1 || quant==2){
645 nom->GetXaxis()->SetLabelSize(0.);
647 if(quant==2 || quant==4){
648 nom->GetYaxis()->SetLabelSize(0.);
650 float minx = 0.99 * nom->GetXaxis()->GetXmin();
651 float maxx = 0.99 * nom->GetXaxis()->GetXmax();
652 nom->GetXaxis()->SetLabelOffset(0.01);
653 nom->GetYaxis()->SetLabelOffset(0.01);
654 nom->GetYaxis()->SetTitleOffset(0.50);
655 nom->GetXaxis()->SetTitleSize(0.0);
656 nom->GetYaxis()->SetTitleSize(0.0);
657 nom->GetXaxis()->SetTitleSize(0.0);
658 nom->GetYaxis()->SetTitleSize(0.0);
659 nom->GetYaxis()->SetRangeUser(0.00001, maxy);
660 nom->GetXaxis()->SetRangeUser(minx, maxx);
661 g->GetYaxis() ->SetRangeUser(0.00001, maxy);
662 g->GetXaxis() ->SetRangeUser(minx, maxx);
663 g->SetFillColor(errCol);
665 nom->Draw(
"hist ][ same");
676 THStack*
ToTHStack(
const std::vector<std::pair<const Spectrum&, Color_t>>&
s,
679 THStack*
ret =
new THStack;
681 TH1*
h =
it.first.ToTH1(pot);
682 h->SetFillColor(
it.second);
691 double x0,
double y0,
double x1,
double y1)
694 if(x > x0 && x < x1 && y > y0 && y < y1)
return 0;
703 if(x > x0 && x < x1){
708 if(y > y0 && y < y1){
722 dx *= (gPad->GetX2()-gPad->GetX1());
723 dy *= (gPad->GetY2()-gPad->GetY1());
726 const double x0 = gPad->GetUxmin();
727 const double x1 = gPad->GetUxmax();
728 const double y0 = gPad->GetUymin();
729 const double y1 = gPad->GetUymax();
731 const double X = x1-x0;
732 const double Y = y1-y0;
738 if(yPin >= 0) besty = yPin * (gPad->GetY2()-gPad->GetY1());
740 for(
bool fallback: {
false,
true}){
741 for(
double x = x0+dx/2;
x < x1-dx/2;
x += X/50){
742 for(
double y = y0+dy/2;
y < y1-dy/2;
y += Y/50){
750 if(d < bestd)
continue;
752 TIter
next(gPad->GetListOfPrimitives());
753 while(TObject* obj =
next()){
755 if(!obj->InheritsFrom(TH1::Class()))
continue;
756 if( obj->InheritsFrom(TH2::Class()))
continue;
760 for(
int n = 1;
n <= h->GetNbinsX(); ++
n){
761 const double px = h->GetBinCenter(
n);
762 const double py = h->GetBinContent(
n);
769 (
x-dx/2)/X, (
y-dy/2)/Y,
770 (
x+dx/2)/X, (
y+dy/2)/Y));
781 if (yPin < 0) besty =
y;
786 if(bestd != 0)
break;
790 const double nx = (bestx-gPad->GetX1())/(gPad->GetX2()-gPad->GetX1());
791 const double ny = (besty-gPad->GetY1())/(gPad->GetY2()-gPad->GetY1());
793 const double ndx = dx/(gPad->GetX2()-gPad->GetX1());
794 const double ndy = dy/(gPad->GetY2()-gPad->GetY1());
796 return new TLegend(nx-ndx/2, ny-ndy/2, nx+ndx/2, ny+ndy/2);
801 const std::pair<double, double> & statErr,
804 bool doStat = statErr.first != 0 && statErr.second != 0;
805 int nbins = systErrors.size() + 1;
809 int firstSystBin = 1;
817 for(
auto it_systPairs: systErrors)
819 systTot.first +=
util::sqr(it_systPairs.second.first);
820 systTot.second +=
util::sqr(it_systPairs.second.second);
822 systTot.first =
sqrt(systTot.first);
823 systTot.second =
sqrt(systTot.second);
827 while(xrange <
std::max(
std::max(systTot.first, systTot.second),
std::max(statErr.first, statErr.second))+2) xrange += 5;
829 TH2*
axes =
new TH2F(
"", (
";" + label).c_str(),
830 10, -xrange, +xrange, nbins, 0, nbins);
831 axes->GetXaxis()->CenterTitle();
834 std::vector<std::pair<std::pair<double, double>,
std::string>> systList;
836 std::sort(systList.begin(),
838 [](
const auto & errNamePair1,
const auto & errNamePair2)
840 return (
std::max(errNamePair1.first.first, errNamePair1.first.second)
841 <
std::max(errNamePair2.first.first, errNamePair2.first.second) );
846 TAxis* yax = axes->GetYaxis();
847 for(
unsigned int i = 0;
i < systList.size(); ++
i)
848 yax->SetBinLabel(
i+firstSystBin+1, systList[
i].second.c_str());
850 yax->SetBinLabel(firstSystBin,
"Total syst. error");
851 if(doStat) yax->SetBinLabel(1,
"Statistical error");
854 yax->SetLabelSize(.06);
859 for(
unsigned int i = 0;
i < systList.size(); ++
i){
860 TBox* box =
new TBox(-systList[
i].first.first, i+firstSystBin+.2,
861 +systList[i].first.second, i+firstSystBin+.8);
862 std::cout<<
"ERROR = (" << systList[
i].first.first <<
", " << systList[
i].first.second <<
")" <<
std::endl;
863 box->SetFillColor(kAzure-8);
868 TBox* syst =
new TBox(-systTot.first, firstSystBin-.8,
869 +systTot.second, firstSystBin-.2);
870 std::cout<<
"TOTAL SYSTEMATIC ERROR = (" << systTot.first <<
", " << systTot.second <<
")" <<
std::endl;
871 syst->SetFillColor(kAzure-8);
876 TBox*
stat =
new TBox(-statErr.first, .2, +statErr.second, .8);
877 stat->SetFillColor(
kRed-7);
882 TLine* div =
new TLine(-xrange, firstSystBin, +xrange, firstSystBin);
883 div->SetLineWidth(2);
887 TLine*
zero =
new TLine(0, 0, 0, nbins);
888 zero->SetLineStyle(7);
889 zero->SetLineWidth(2);
893 gPad->SetLeftMargin(.25);
900 std::map<std::string, double>
systs = systbars;
902 std::map<std::string, std::pair<double, double>> systVals;
903 for (
const auto & systPair : systbars)
904 systVals[systPair.first] =
std::make_pair(systPair.second, systPair.second);
914 TGraphAsymmErrors* gr =
new TGraphAsymmErrors(h);
916 TFeldmanCousins
fc(0.6827);
918 for(
int i = 0;
i < h->GetNbinsX(); ++
i){
920 gr->GetPoint(
i, x, y);
922 if ( drawEmptyBins || y!=0 ){
923 if ( y < 50 ) gr->SetPointEYlow(
i, y-fc.CalculateLowerLimit(y,0));
924 else gr->SetPointEYlow(
i,
sqrt(y));
925 if ( y < 30 ) gr->SetPointEYhigh(
i,fc.CalculateUpperLimit(y,0)-
y);
926 else gr->SetPointEYhigh(
i,
sqrt(y));
930 gr->SetPointEXlow(
i, 0);
931 gr->SetPointEXhigh(
i, 0);
934 gr->SetMarkerStyle(kFullCircle);
935 gr->SetMarkerColor(1);
936 gr->SetMarkerSize(1);
946 TGraphAsymmErrors* gr =
new TGraphAsymmErrors(h);
947 TGraphAsymmErrors*
gr2 =
new TGraphAsymmErrors();
949 TFeldmanCousins
fc(0.6827);
952 for(
int i = 0;
i < h->GetNbinsX(); ++
i){
954 gr->GetPoint(
i, x, y);
958 if ( (drawEmptyBins && h2->GetBinContent(
i+1) > 0.05) || y!=0 ){
959 gr2->SetPoint(j, x, y);
960 gr2->SetPointEXlow(j, gr->GetErrorXlow(
i));
961 gr2->SetPointEXhigh(j, gr->GetErrorXhigh(
i));
962 if ( y < 50 ) gr2->SetPointEYlow(j, y-fc.CalculateLowerLimit(y,0));
963 else gr2->SetPointEYlow(j,
sqrt(y));
964 if ( y < 30 ) gr2->SetPointEYhigh(j,fc.CalculateUpperLimit(y,0)-
y);
965 else gr2->SetPointEYhigh(j,
sqrt(y));
970 gr2->SetPointEXlow(
i, 0);
971 gr2->SetPointEXhigh(
i, 0);
982 for(
int bin = 1;
bin <= hist->GetNbinsX();
bin++) {
983 if(hist->GetBinContent(
bin)==0){
984 hist->SetBinContent(
bin,0.001);
985 histScaled->SetBinContent(
bin,0.001);
989 TFeldmanCousins
fc(0.6827);
990 TGraphAsymmErrors *datapoisson =
new TGraphAsymmErrors(histScaled);
992 for(
int bin = 1;
bin <= hist->GetNbinsX();
bin++) {
994 double binWidth = hist->GetBinWidth(
bin);
995 double scaleFactor = overallScale / binWidth;
996 double y = hist->GetBinContent(
bin);
1003 double errup = 0, errdn = 0;
1004 if(y < 30) errup = (fc.CalculateUpperLimit(y,0) -
y);
1005 else errup =
sqrt(y);
1006 if(y < 50) errdn = (y - fc.CalculateLowerLimit(y,0));
1007 else errdn =
sqrt(y);
1009 datapoisson->SetPointEYhigh(
bin-1, scaleFactor * errup);
1010 datapoisson->SetPointEYlow (
bin-1, scaleFactor * errdn);
1012 datapoisson->SetMarkerStyle(kFullCircle);
1013 datapoisson->SetMarkerColor(1);
1014 datapoisson->SetMarkerSize(1);
1015 datapoisson->SetLineWidth(2);
1024 for(
int bin = 1;
bin <= data->GetNbinsX();
bin++) {
1025 if(data->GetBinContent(
bin)==0){
1026 data->SetBinContent(
bin,0.001);
1030 TFeldmanCousins
fc(0.6827);
1031 TGraphAsymmErrors *gr =
new TGraphAsymmErrors(data);
1033 for(
int bin = 1;
bin <= data->GetNbinsX();
bin++) {
1035 double width = data->GetBinWidth(
bin);
1036 double scale = width / overallScale;
1039 gr->GetPoint(
bin-1, x, y);
1040 double bkg = bkgd->GetBinContent(
bin);
1049 double errup = 0, errdn = 0;
1050 if(y < 30) errup = (fc.CalculateUpperLimit(y,bkg) -
y);
1051 else errup =
sqrt(y);
1052 if(y < 50) errdn = (y - fc.CalculateLowerLimit(y,bkg));
1053 else errdn =
sqrt(y);
1055 gr->SetPointEYhigh(
bin-1, 1/scale * errup);
1056 gr->SetPointEYlow (
bin-1, 1/scale * errdn);
1059 gr->SetMarkerStyle(kFullCircle);
1060 gr->SetMarkerColor(1);
1061 gr->SetMarkerSize(1);
1062 gr->SetLineWidth(2);
1072 TFeldmanCousins
fc(0.6827);
1076 for (
int bin = 1;
bin <= data->GetNbinsX(); ++
bin) {
1078 gRat->GetPoint(
bin-1, x, y);
1079 double den = pred->GetBinContent(
bin);
1086 if (y <= 0) y = 0.01;
1088 double errup = 0, errdn = 0;
1089 if (y < 30) errup = (fc.CalculateUpperLimit(y, 0) -
y) / den;
1091 if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, 0)) / den;
1094 gRat->SetPoint(
bin-1, x, y / den);
1095 gRat->SetPointEYhigh(
bin-1, errup);
1096 gRat->SetPointEYlow(
bin-1, errdn);
1098 gRat->SetMarkerStyle(kFullCircle);
1099 gRat->SetMarkerColor(1);
1100 gRat->SetMarkerSize(1);
1101 gRat->SetLineWidth(2);
1110 TFeldmanCousins
fc(0.6827);
1114 for (
int bin = 1;
bin <= data->GetNbinsX(); ++
bin) {
1115 double width = data->GetBinWidth(
bin);
1116 double scale = width / overallScale;
1119 gRat->GetPoint(
bin-1, x, y);
1120 double den = pred->GetBinContent(
bin);
1129 if (y <= 0) y = 0.01;
1131 double errup = 0, errdn = 0;
1132 if (y < 30) errup = (fc.CalculateUpperLimit(y, 0) -
y) / den;
1134 if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, 0)) / den;
1137 gRat->SetPoint(
bin-1, x, y / den);
1138 gRat->SetPointEYhigh(
bin-1, errup);
1139 gRat->SetPointEYlow (
bin-1, errdn);
1141 gRat->SetMarkerStyle(kFullCircle);
1142 gRat->SetMarkerColor(1);
1143 gRat->SetMarkerSize(1);
1144 gRat->SetLineWidth(2);
1153 TFeldmanCousins
fc(0.6827);
1157 for (
int bin = 1;
bin <= data->GetNbinsX(); ++
bin) {
1159 gRat->GetPoint(
bin-1, x, y);
1160 double den = pred->GetBinContent(
bin);
1161 double bkg = bkgd->GetBinContent(
bin);
1168 if (y <= 0) y = 0.01;
1170 double errup = 0, errdn = 0;
1171 if (y < 30) errup = (fc.CalculateUpperLimit(y, bkg) -
y) / den;
1173 if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, bkg)) / den;
1176 gRat->SetPoint(
bin-1, x, y / den);
1177 gRat->SetPointEYhigh(
bin-1, errup);
1178 gRat->SetPointEYlow (
bin-1, errdn);
1180 gRat->SetMarkerStyle(kFullCircle);
1181 gRat->SetMarkerColor(1);
1182 gRat->SetMarkerSize(1);
1183 gRat->SetLineWidth(2);
1193 TFeldmanCousins
fc(0.6827);
1197 for (
int bin = 1;
bin <= data->GetNbinsX(); ++
bin) {
1198 double width = data->GetBinWidth(
bin);
1199 double scale = width / overallScale;
1202 gRat->GetPoint(
bin-1, x, y);
1203 double den = pred->GetBinContent(
bin);
1204 double bkg = bkgd->GetBinContent(
bin);
1215 if (y <= 0) y = 0.01;
1217 double errup = 0, errdn = 0;
1218 if (y < 30) errup = (fc.CalculateUpperLimit(y, bkg) -
y) / den;
1220 if (y < 50) errdn = (y - fc.CalculateLowerLimit(y, bkg)) / den;
1223 gRat->SetPoint(
bin-1, x, y / den);
1224 gRat->SetPointEYhigh(
bin-1, errup);
1225 gRat->SetPointEYlow (
bin-1, errdn);
1227 gRat->SetMarkerStyle(kFullCircle);
1228 gRat->SetMarkerColor(1);
1229 gRat->SetMarkerSize(1);
1230 gRat->SetLineWidth(2);
1238 int n = hmin->GetNbinsX();
1239 TGraph* gr =
new TGraph(4*n);
1241 for(
int i=1;
i<=
n;
i++)
1243 double xdown = hmax->GetBinLowEdge(
i);
1244 double xup = hmax->GetBinLowEdge(
i+1);
1245 double ydown = hmin->GetBinContent(
i);
1246 double yup = hmax->GetBinContent(
i);
1248 gr->SetPoint(2*(
i-1), xdown, ydown);
1249 gr->SetPoint(2*(
i-1)+1, xup, ydown);
1250 gr->SetPoint(4*n-2*
i, xup, yup);
1251 gr->SetPoint(4*n-2*
i+1, xdown, yup);
1254 gr->SetFillColor(hmin->GetLineColor() - 9);
1255 gr->SetLineColor(hmin->GetLineColor() - 9);
1263 const std::pair<double, double> & quantileDivisions)
1265 if (hist->GetDimension() != 2)
1266 throw std::runtime_error(Form(
"Can't profile a histogram with other than 2 dimensions. Yours is a %d-D histogram...", hist->GetDimension()));
1268 TAxis
const*
axis =
nullptr;
1269 std::function<TH1D*(const char *, Int_t, Int_t, Option_t*)> projectionMethod;
1270 if (axisName ==
"x" || axisName ==
"X")
1272 axis = hist->GetXaxis();
1273 projectionMethod = std::bind(&
TH2::ProjectionY, hist, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
1275 else if (axisName ==
"y" || axisName ==
"Y")
1277 axis = hist->GetYaxis();
1278 projectionMethod = std::bind(&
TH2::ProjectionX, hist, std::placeholders::_1, std::placeholders::_2, std::placeholders::_3, std::placeholders::_4);
1281 throw std::runtime_error( Form(
"Can't profile onto unknown axis: '%s'", axisName.c_str()) );
1284 std::vector<double> points_x, points_y, errors_x_low, errors_x_high, errors_y_low, errors_y_high;
1287 double _quantileDivisions[2] = {quantileDivisions.first, quantileDivisions.second};
1288 for (
int binNum = 1; binNum <= axis->GetNbins(); binNum++)
1290 std::unique_ptr<TH1D>
proj ( projectionMethod(Form(
"%s_tmp_%d", hist->GetName(),
int(
std::time(
nullptr))), binNum, binNum,
"") );
1292 double y =
proj->GetMean();
1293 points_x.push_back(axis->GetBinCenter(binNum));
1294 points_y.push_back(y);
1297 unsigned int n_qs = 0;
1298 if (
proj->Integral() == 0)
1301 quantiles[0] = quantiles[1] = 0;
1304 n_qs =
proj->GetQuantiles(2, quantiles, _quantileDivisions);
1306 throw std::runtime_error( Form(
"GetQuantiles() didn't compute all the quantiles in HistoTools.ProfileQuantile(). I requested 2 quantiles, but got %d of them...", n_qs) );
1308 double binWidth = axis->GetBinWidth(binNum);
1309 errors_x_low.push_back(binWidth/2.);
1310 errors_x_high.push_back(binWidth/2.);
1311 errors_y_low.push_back(y-quantiles[0]);
1312 errors_y_high.push_back(quantiles[1]-y);
1315 TGraphAsymmErrors * outGraph =
new TGraphAsymmErrors(
1324 std::string name = (graphName.size()) ? graphName : Form(
"%s_quantile_errors", hist->GetName());
1325 outGraph->SetName(name.c_str());
1336 TMarker* manMarker =
new TMarker(bfSin, bfDm, marker);
1337 manMarker->SetMarkerSize(
size);
1338 manMarker->SetMarkerColor(color);
1347 double mirror = bfSin;
1348 TMarker* manMarker =
new TMarker(bfSin, bfDm, marker);
1349 manMarker->SetMarkerSize(
size);
1350 manMarker->SetMarkerColor(color);
1353 if(bfSin > 0.514) mirror = 0.514 - (bfSin - 0.514);
1354 if(bfSin < 0.514) mirror = 0.514 + (0.514 - bfSin );
1355 TMarker* mirrorMarker =
new TMarker(mirror, bfDm, marker);
1356 mirrorMarker->SetMarkerSize(
size);
1357 mirrorMarker->SetMarkerColor(color);
1358 mirrorMarker->Draw();
1369 hist->GetXaxis()->SetLabelOffset(0.01);
1370 hist->GetYaxis()->SetLabelOffset(0.01);
1371 hist->GetYaxis()->SetTitleOffset(0.50);
1372 hist->GetXaxis()->SetTitleSize(0.0);
1373 hist->GetYaxis()->SetTitleSize(0.0);
1374 hist->GetYaxis()->SetRangeUser(0.0, ymax);
1376 if(quant==1 || quant==2) hist->GetXaxis()->SetLabelSize(0.);
1377 if(quant==2 || quant==4) hist->GetYaxis()->SetLabelSize(0.);
1385 double optmin = 0., optmax = 0., optwi = 0.;
1388 double xmin = hist->GetXaxis()->GetXmin();
1389 double xmax = hist->GetXaxis()->GetXmax();
1390 int xdiv = hist->GetXaxis()->GetNdivisions() % 100;
1391 THLimitsFinder::Optimize(xmin, xmax, xdiv, optmin, optmax, optstp, optwi);
1393 if ( (xmax - optmax) < optwi/2 )
1394 hist->GetXaxis()->ChangeLabel( -1, -1, 0);
1396 optmin = 0.; optmax = 0.; optwi = 0.; optstp = 0;
1398 int ydiv = hist->GetYaxis()->GetNdivisions() % 100;
1399 THLimitsFinder::Optimize(ymin, ymax, ydiv, optmin, optmax, optstp, optwi);
1401 hist->GetYaxis()->SetRangeUser(ymin, ymax);
1402 if ( (ymax - optmax) < optwi/2 )
1403 hist->GetYaxis()->ChangeLabel( -1, -1, 0);
1408 void PimpHist(TH1*
hist, Style_t linestyle, Color_t linecolor,
int linewidth, Style_t markerstyle, Color_t markercolor,
double markersize){
1409 hist->SetLineColor(linecolor);
1410 hist->SetLineStyle(linestyle);
1411 hist->SetLineWidth(linewidth);
1412 hist->SetMarkerColor(markercolor);
1413 hist->SetMarkerStyle(markerstyle);
1414 hist->SetMarkerSize(markersize);
1421 if(!gPad)
new TCanvas;
1422 p1 =
new TPad(
"",
"", 0, 0, 1, 1);
1423 p2 =
new TPad(
"",
"", 0, 0, 1, 1);
1425 p1->SetBottomMargin(ysplit);
1426 p2->SetTopMargin(1-ysplit);
1430 for(TPad*
p: {
p2, p1}){
1441 canvas =
new TCanvas(
UniqueName().c_str(),
"",960,800);
1443 pad3 =
new TPad(
UniqueName().c_str(),
"pad3",0,0,1,1);
1445 pad3->SetTopMargin(0.5);
1446 pad3->SetBottomMargin(0.1);
1447 pad3->SetLeftMargin(0.1);
1448 pad3->SetRightMargin(0.49);
1449 pad3->SetFillStyle(0);
1452 pad4 =
new TPad(
UniqueName().c_str(),
"pad4",0,0,1,1);
1454 pad4->SetTopMargin(0.5);
1455 pad4->SetBottomMargin(0.1);
1456 pad4->SetLeftMargin(0.51);
1457 pad4->SetRightMargin(0.08);
1458 pad4->SetFillStyle(0);
1461 pad1 =
new TPad(
UniqueName().c_str(),
"pad1",0,0,1,1);
1463 pad1->SetTopMargin(0.1);
1464 pad1->SetBottomMargin(0.5);
1465 pad1->SetLeftMargin(0.1);
1466 pad1->SetRightMargin(0.49);
1467 pad1->SetFillStyle(0);
1470 pad2 =
new TPad(
UniqueName().c_str(),
"pad2",0,0,1,1);
1472 pad2->SetTopMargin(0.1);
1473 pad2->SetBottomMargin(0.5);
1474 pad2->SetLeftMargin(0.51);
1475 pad2->SetRightMargin(0.08);
1476 pad2->SetFillStyle(0);
1485 histo->GetXaxis()->CenterTitle();
1486 histo->GetYaxis()->CenterTitle();
1487 histo->GetZaxis()->CenterTitle();
1505 auto h =
new TH1D (
"sh",
";; Pull (N#sigma)",
1506 nsysts, -0.5, nsysts + 0.5);
1508 for (
int systIdx = 0; systIdx <
nsysts; ++systIdx){
1510 h->SetBinContent( systIdx + 1, shiftThis);
1513 if( tempName.Contains(
"Neutron") ) h->GetXaxis()->SetBinLabel( systIdx + 1,
"NeutronSyst2018");
1514 else h->GetXaxis()->SetBinLabel( systIdx + 1,
systs[systIdx]->ShortName().c_str());
1517 h->SetMarkerStyle(kFullCircle);
1519 h->SetTitleOffset(3.0);
1520 h->GetXaxis()->SetNdivisions(nsysts,kFALSE);
1521 h->GetXaxis()->SetLabelSize(0.065);
1522 h->GetXaxis()->LabelsOption(
"v");
1523 h->GetYaxis()->SetRangeUser(-1,1);
1524 auto c1 =
new TCanvas (
"c1",
"c1",700,350);
1525 c1->SetBottomMargin(0.5);
1534 for(
int i = 0;
i < b->GetN(); ++
i){
1536 b->GetPoint(
i, x, y);
1537 ret->SetPoint(ret->GetN(),
x,
y);
1540 ret->SetFillColor(fillcol);
1548 g->SetPoint(g->GetN(),
xmax,
y);
1549 g->SetPoint(g->GetN(),
xmin,
y);
1551 g->SetFillColor(col);
1560 TPaveText *pText =
new TPaveText (0.1, 0.90, 0.16, 0.94,
"NDC");
1561 pText->SetFillStyle(0);
1562 pText->SetLineColor(0);
1563 pText->SetLineWidth(0);
1564 pText->SetBorderSize(1);
1565 pText->SetTextColor(kGray+2);
1566 if (isFHC) pText->AddText(
"#nu-beam");
1567 else pText->AddText(
"#bar{#nu}-beam");
1568 pText->SetTextSize(2/40.);
1569 pText->SetTextAlign(11);
1583 else if (quant == 1) {
1584 ret =
new TPaveText(0.37, 0.75, 0.47, 0.90,
"NDC");
1585 ret->SetBorderSize(0);
1586 ret->SetFillStyle(0);
1587 TText *
text = ret->AddText(
TString::Format(
"#splitline{#bf{Quartile %d}}{best resolution}", quant));
1588 text ->SetTextAlign(22);
1589 text ->SetTextSize(0.025);
1592 else if (quant == 2) {
1593 ret =
new TPaveText(0.75, 0.75, 0.85, 0.90,
"NDC");
1594 ret->SetBorderSize(0);
1595 ret->SetFillStyle(0);
1597 text ->SetTextAlign(22);
1598 text ->SetTextSize(0.025);
1601 else if (quant == 3) {
1602 ret =
new TPaveText(0.37, 0.35, 0.47, 0.50,
"NDC");
1603 ret->SetBorderSize(0);
1604 ret->SetFillStyle(0);
1606 text ->SetTextAlign(22);
1607 text ->SetTextSize(0.025);
1610 else if (quant == 4) {
1611 ret =
new TPaveText(0.75, 0.35, 0.85, 0.50,
"NDC");
1612 ret->SetBorderSize(0);
1613 ret->SetFillStyle(0);
1614 TText *
text = ret->AddText(
TString::Format(
"#splitline{#bf{Quartile %d}}{worst resolution}", quant));
1615 text ->SetTextAlign(22);
1616 text ->SetTextSize(0.025);
1620 std::cout <<
"\n***********************************************\n" 1621 <<
"\n No valid quantile chosen for DrawQuantLabel \n" 1622 <<
"\n***********************************************\n" T max(const caf::Proxy< T > &a, T b)
void PlotWithAreaSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype)
Plot prediction with +/-1sigma shape-only error band.
std::vector< double > quantiles(TH1D *h)
TSpline3 lo("lo", xlo, ylo, 12,"0")
void DataMCRatio(const Spectrum &data, const IPrediction *mc, osc::IOscCalc *calc, double miny, double maxy)
Plot data/MC ratio for the given spectrum. Normalize MC to Data by POT.
Cuts and Vars for the 2020 FD DiF Study.
std::map< std::string, double > xmax
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
TH1 * PlotSystShifts(const SystShifts &shifts, bool sortName)
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
void MakeHistCanvasReady_Quant(const int quant, TH1 *hist, double ymax)
Eigen::ArrayXd ProjectionX(const Eigen::MatrixXd &mat)
Helper for Unweighted.
Float_t y1[n_points_granero]
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.
virtual const std::string & ShortName() const final
The name printed out to the screen.
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Simple record of shifts applied to systematic parameters.
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN'T A BETTER WAY!
Float_t x1[n_points_granero]
TGraphAsymmErrors * PlotWithSystErrorBand(IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float headroom, bool newaxis, EBinType bintype, double alpha)
Plot prediction with +/-1sigma error band.
General interface to oscillation calculators.
virtual Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &syst, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
TPaveText * DrawQuantLabel(int quant)
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
void CenterTitles(TH1 *histo)
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
void SplitCanvas(double ysplit, TPad *&p1, TPad *&p2)
Split the current pad into two vertically stacked pieces, suitable for ratio plots and so on...
TCanvas * canvas(const char *nm, const char *ti, const char *a)
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
const Color_t kTotalMCErrorBandColor
static SystShifts Nominal()
fvar< T > round(const fvar< T > &x)
::xsd::cxx::tree::time< char, simple_type > time
T sqr(T x)
More efficient square function than pow(x,2)
Encapsulate code to systematically shift a caf::SRProxy.
const Color_t kNumuBackgroundColor
void GetSystBands(osc::IOscCalc *calc, const IPrediction *pred, std::vector< const ISyst * > systs, std::vector< TH1 * > &hUps, std::vector< TH1 * > &hDns, double pot, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign)
TGraph * RatioAsymmError(TH1 *data, TH1 *pred)
TGraphAsymmErrors * ProfileQuantile(const TH2 *hist, const std::string &axisName, const std::string &graphName, const std::pair< double, double > &quantileDivisions)
Calculate profile with error bars corresponding to specified quantiles of a 2D distribution (by defau...
Representation of a spectrum in any variable, with associated POT.
void PimpHist(TH1 *hist, Style_t linestyle, Color_t linecolor, int linewidth, Style_t markerstyle, Color_t markercolor, double markersize)
Pimp histogram once and for all.
const XML_Char const XML_Char * data
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
void drawBFSingle(double bfSin, double bfDm, Color_t color, Style_t marker, double size=1.5)
Draw a single BF point.
Charged-current interactions.
TSpline3 hi("hi", xhi, yhi, 18,"0")
void PlotWithSystErrorBandTwoPreds(const Spectrum &nominal, const std::vector< Spectrum > &upShifts, const std::vector< Spectrum > &downShifts, const Spectrum &alternative, const std::vector< Spectrum > &upShiftsAlt, const std::vector< Spectrum > &downShiftsAlt, double pot, int col1, int col2, float headroom, bool newaxis, EBinType bintype)
Plot two different predictions with +/-1sigma shape-only error bands.
TH1 * DataMCComparisonAreaNormalized(const Spectrum &data, const Spectrum &mc)
TGraph * graphAsymmErrorScaled(TH1 *histScaled, TH1 *hist, double overallScale)
TPaveText * DrawBeamLabel(bool isFHC)
Put the standardized beam label in the left corner of the active canvas.
const Color_t kBeamNueBackgroundColor
TGraph * ShadeBetweenHistograms(TH1 *hmin, TH1 *hmax)
TGraphAsymmErrors * GraphWithPoissonErrors2(const TH1 *h, const TH1 *h2, bool noErrorsXaxis, bool drawEmptyBins)
Same as above but use a reference histogram to determine which empty bins to draw.
T GetShift(const ISyst *syst) const
virtual Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &syst) const
TGraph * ExtendGraphToTop(TGraph *g, int col, double xmin, double xmax, double y)
Represent the ratio between two spectra.
TGraph * JoinGraphs(TGraph *a, TGraph *b, int fillcol)
Join graphs and set a fill color. Useful for contours.
static float min(const float a, const float b, const float c)
double PointDistanceToBox(double x, double y, double x0, double y0, double x1, double y1)
Helper for AutoPlaceLegend.
TGraph * RatioAsymmErrorWithBkg(TH1 *data, TH1 *pred, TH1 *bkgd)
void CountingExperimentErrorBarChart(const std::map< std::string, double > &systbars, double statErr, bool bkgdOrSig, bool shortchart)
Make a simple plot of relative size of different errors.
void DataMCAreaNormalizedRatio(const Spectrum &data, const IPrediction *mc, osc::IOscCalc *calc, double miny, double maxy)
Plot data/MC ratio for the given spectrum. Normalize MC to Data by area.
TGraph * graphAsymmErrorWithBkgScaled(TH1 *data, TH1 *bkgd, double overallScale)
TH1 * DataMCComparisonComponents(const Spectrum &data, const IPrediction *mc, osc::IOscCalc *calc)
Plot MC broken down into flavour components, overlayed with data.
Eigen::ArrayXd ProjectionY(const Eigen::MatrixXd &mat)
Helper for WeightingVariable.
TGraphAsymmErrors * PlotWithSystErrorBand_Quant(const int quant, IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float maxy, bool newaxis)
void GetBFSystBands(osc::IOscCalc *calc, const IPrediction *pred, std::vector< const ISyst * > systs, const SystShifts &bfshifts, std::vector< TH1 * > &hUps, std::vector< TH1 * > &hDns, double pot, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign)
TGraph * RatioAsymmErrorWithBkgScaled(TH1 *data, TH1 *pred, TH1 *bkgd, double overallScale)
Neutral-current interactions.
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.
THStack * ToTHStack(const std::vector< std::pair< const Spectrum &, Color_t >> &s, double pot)
Can call like ToTHStack({{h1, kRed}, {h2, kBlue}}, pot)
bool SortSystsName(const ISyst *systA, const ISyst *systB)
Both neutrinos and antineutrinos.
Standard interface to all prediction techniques.
std::vector< const ISyst * > ActiveSysts() const
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
const Color_t kTotalMCColor
void SplitCanvasQuant(TCanvas *&canvas, TPad *&pad1, TPad *&pad2, TPad *&pad3, TPad *&pad4)
TGraph * RatioAsymmErrorScaled(TH1 *data, TH1 *pred, double overallScale)
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
const Color_t kNCBackgroundColor
TCanvas * RatioPlot(std::vector< TH1 * > topHistos, std::vector< TString > topOption, std::vector< TH1 * > bottomHistos, std::vector< TString > bottomOption, TString pidtag, bool pidaxis=false)
void SetShift(const ISyst *syst, double shift, bool force=false)
std::string UniqueName()
Return a different string each time, for creating histograms.
void drawBFMirror(double bfSin, double bfDm, Color_t color, Style_t marker, double size=1.5)
Draw best fit at both octants. Truth for neutrinos, not sure about antineutrinos. ...