16 #include "TGraphAsymmErrors.h" 35 TLatex*
prelim =
new TLatex(.9, .95,
"NOvA Preliminary");
36 prelim->SetTextColor(
kBlue);
38 prelim->SetTextSize(2/36.);
39 prelim->SetTextAlign(32);
45 histo->GetXaxis()->CenterTitle();
46 histo->GetYaxis()->CenterTitle();
47 histo->GetZaxis()->CenterTitle();
52 TLatex*
tag =
new TLatex(.83, .92, pid);
53 tag->SetTextColor(kAzure+3);
55 tag->SetTextSize(0.038);
57 tag->SetTextAlign(12);
63 TLatex*
tag =
new TLatex(.18, .825, pid);
64 if ( pid ==
"IH") tag->SetTextColor(
kOrange+2);
65 if ( pid ==
"NH") tag->SetTextColor(kAzure+3);
67 tag->SetTextSize(0.043);
69 tag->SetTextAlign(12);
78 ret->GetXaxis()->SetTitle(h->GetXaxis()->GetTitle());
79 ret->GetYaxis()->SetTitle(h->GetYaxis()->GetTitle());
81 std::vector<std::pair<int, int>> binMap = {{3, 2},
94 for(
auto it: binMap) ret->SetBinContent(
it.second, h->GetBinContent(
it.first));
102 std::vector<TString>
in = {
"0",
"1",
"2",
"3",
"4",
"5",
103 "6",
"7",
"8",
"9",
"_",
" ",
"(",
")"};
104 std::vector<TString>
out = {
"zero",
"one",
"two",
"three",
"four",
"five",
105 "six",
"seven",
"eight",
"nine",
"",
"",
"",
""};
106 for(
unsigned int i=0;
i<in.size();
i++)
107 mystring.ReplaceAll(in[
i],out[i]);
115 TGraph*
l1 =
new TGraph;
116 l1->SetPoint(0, 10, 0);
117 l1->SetPoint(1, 10, 50000);
118 TGraph* l2 =
new TGraph;
119 l2->SetPoint(0, 20, 0);
120 l2->SetPoint(1, 20, 50000);
121 for(TGraph*
l: {
l1, l2}){
123 l->SetLineColor(kGray+3);
131 const double y = gPad->GetUymax() * 0.955;
139 for(
int i = 0;
i < 3; ++
i){
145 ltx->SetTextSize(0.032);
146 ltx->SetTextColor(kGray+3);
147 ltx->SetTextAlign(22);
156 double xmax = axes->GetXaxis()->GetXmax();
157 double xmin = axes->GetXaxis()->GetXmin();
158 double xrange = xmax-
xmin;
160 double ymin = axes->GetMinimum();
161 TGaxis* gax1 =
new TGaxis(0, ymin, 8, ymin, 0, 4, 4,
"");
162 TGaxis* gax2 =
new TGaxis(10, ymin, 18, ymin, 0, 4, 4,
"");
163 TGaxis* gax3 =
new TGaxis(20, ymin, 30, ymin, 0, 5, 5,
"");
165 float size = axes->GetXaxis()->GetLabelSize();
166 float offs = axes->GetXaxis()->GetLabelOffset();
167 int font = axes->GetXaxis()->GetLabelFont();
169 axes->GetXaxis()->SetLabelSize(0);
171 axes->GetXaxis()->SetNdivisions(3000);
172 axes->GetXaxis()->SetTitle(
"Reconstructed Neutrino Energy (GeV)");
174 axes->GetYaxis()->SetTitleOffset(0.8);
179 gax1->SetLabelSize(size);
180 gax1->SetLabelOffset(offs);
181 gax1->SetLabelFont(font);
185 for(TGaxis*
ax: {gax1, gax2, gax3}){
187 ax->SetLabelSize(size-0.001);
188 ax->SetLabelOffset(offs);
189 ax->SetLabelFont(font);
193 TLine* xline =
new TLine();
194 xline->SetLineWidth(2);
195 xline->DrawLine(xmin, 0, xmax, 0);
203 TGraph*
l1 =
new TGraph;
204 l1->SetPoint(0, 6, 0);
205 l1->SetPoint(1, 6, 50000);
206 TGraph* l2 =
new TGraph;
207 l2->SetPoint(0, 12, 0);
208 l2->SetPoint(1, 12, 50000);
209 for(TGraph*
l: {
l1, l2}){
212 l->SetLineColor(kGray+3);
217 TGaxis* gax1 =
new TGaxis(6, 0, 6, 18, 0, 18, 504,
"U+-");
218 TGaxis* gax2 =
new TGaxis(12, 0, 12, 18, 0, 18, 504,
"U+-");
228 const double y = gPad->GetUymax() * .94;
236 for(
int i = 0;
i < 3; ++
i){
241 ltx->SetTextSize(.04);
242 ltx->SetTextAlign(22);
248 double ymin = hist->GetMinimum();
250 TGaxis* gax1 =
new TGaxis(1, ymin, 5, ymin, 1, 3, 3,
"");
251 TGaxis* gax2 =
new TGaxis(7, ymin, 11, ymin, 1, 3, 3,
"");
252 TGaxis* gax3 =
new TGaxis(13, ymin, 17, ymin, 1, 3, 3,
"");
254 float tick = hist->GetXaxis()->GetTickLength();
255 float size = hist->GetXaxis()->GetLabelSize();
256 float offs = hist->GetXaxis()->GetLabelOffset();
257 int font = hist->GetXaxis()->GetLabelFont();
259 hist->GetXaxis()->SetLabelSize(0);
261 hist->GetXaxis()->SetNdivisions(1800);
262 hist->GetXaxis()->SetTitle(
"Reconstructed neutrino energy (GeV)");
263 hist->GetYaxis()->SetTitle(
"Events / 0.5 GeV Bin");
264 hist->GetYaxis()->SetTitleOffset(0.8);
265 hist->GetXaxis()->CenterTitle();
266 hist->GetYaxis()->CenterTitle();
270 gax1->SetTickSize(tick);
271 gax1->SetLabelSize(size);
272 gax1->SetLabelOffset(offs);
273 gax1->SetLabelFont(font);
277 for(TGaxis*
ax: {gax1, gax2, gax3}){
279 ax->SetTickSize(tick);
280 ax->SetLabelSize(size);
281 ax->SetLabelOffset(offs);
282 ax->SetLabelFont(font);
288 double FindLimitY (std::vector<TH1*> histos,
bool doMax =
true)
290 double max = histos[0]->GetMaximum();
291 double min = histos[0]->GetMinimum();
293 for(
unsigned int i=0;
i<histos.size();
i++){
295 if ( histos[
i]->GetMinimum() <
min ) min = histos[
i]->GetMinimum();
298 if ( doMax )
return max;
307 ,TString ltxcommand =
""){
308 unsigned int nentries = mcnom.size();
309 if(mcshift.size()!=nentries || labels.size()!=
nentries){
std::cerr <<
"incompatible sizes\n";
return;}
312 std::cout <<
"{\\begin{tabular}{ l | r r r } \n" 313 << std::setw(20) <<
" " <<
" & " 314 << std::setw(10) <<
"Nominal" <<
" & " 315 << std::setw(10) <<
"Shift" <<
" & " 316 << std::setw(12) <<
"\\% Diff." 317 <<
"\t \\\\ \\hline \n";
319 std::cout.unsetf(std::ios_base::floatfield);
321 << std::setw(20) << labels[
i] <<
" & " 322 << std::setw(10) << mcnom[
i] ->Integral() <<
" & " 323 << std::setw(10) << mcshift[
i]->Integral() <<
" & " 324 << std::setw(10) << std::setprecision(1) << std::fixed
325 << ( mcnom[
i]->Integral()>0 ? (((mcshift[
i]->Integral()/mcnom[
i]->Integral()) -1.)*100.):0.)
332 void ComparisonTable(std::vector<TH1*> mcnom, std::vector<TH1*> mcplus, std::vector<TH1*> mcminus,
333 std::vector<TString>
labels,TString ltxcommand =
""){
334 unsigned int nentries = mcnom.size();
335 if(mcplus.size()!=nentries || mcminus.size()!=nentries || labels.size()!=
nentries){
std::cerr <<
"incompatible sizes\n";
return;}
337 std::cout <<
"{ \\begin{tabular}{ l | r r r | r r } \n" 338 << std::setw(20) <<
" " <<
" & " 339 << std::setw(10) <<
"Nominal" <<
" & " 340 << std::setw(10) <<
"Shift (+)" <<
" & " 341 << std::setw(10) <<
"Shift (-)" <<
" & " 342 << std::setw(12) <<
"\\% Diff. (+)" <<
" & " 343 << std::setw(12) <<
"\\% Diff. (-)" 344 <<
"\t \\\\ \\hline \n";
346 std::cout.unsetf(std::ios_base::floatfield);
348 << std::setw(20) << labels[
i] <<
" & " 349 << std::setw(10) << mcnom[
i] ->Integral() <<
" & " 350 << std::setw(10) << mcplus[
i] ->Integral() <<
" & " 351 << std::setw(10) << mcminus[
i]->Integral() <<
" & " 352 << std::setw(10) << std::setprecision(1) << std::fixed
353 << ( mcnom[
i]->Integral()>0 ? (((mcplus [
i]->Integral()/mcnom[
i]->Integral()) -1.)*100.):0.)
356 << ( mcnom[
i]->Integral()>0 ? (((mcminus[
i]->Integral()/mcnom[
i]->Integral()) -1.)*100.):0.)
370 leg->SetFillStyle(0);
372 leg->SetX1(kLegendDefaultX1); leg->SetX2(kLegendDefaultX2);
373 leg->SetY1(kLegendDefaultY1); leg->SetY2(kLegendDefaultY2);
376 leg->SetX1(0.405); leg->SetX2(0.565);
377 leg->SetY1(0.59); leg->SetY2(0.82);
378 leg->SetBorderSize(0);
379 leg->SetTextSize(0.037);
381 if(
opt==
"ExPIDPads"){
382 leg->SetX1(0.55); leg->SetX2(0.89);
383 leg->SetY1(0.75); leg->SetY2(0.89);
384 leg->SetTextSize(0.045);
389 double x1=kLegendDefaultX1,
double y1=kLegendDefaultY1,
390 double x2=kLegendDefaultX2,
double y2=kLegendDefaultY2)
393 leg->SetFillStyle(0);
394 TH1* hdata =
new TH1F(
"",
"", 1, 0, 1);
395 TH1*
hsig =
new TH1F(
"",
"", 1, 0, 1);
396 TH1* hNC =
new TH1F(
"",
"", 1, 0, 1);
397 TH1* hCC =
new TH1F(
"",
"", 1, 0, 1);
398 TH1* hBeam =
new TH1F(
"",
"", 1, 0, 1);
399 TH1* htot =
new TH1F(
"",
"", 1, 0, 1);
400 hdata->SetMarkerStyle(kFullCircle);
406 if(sig) leg->AddEntry(hsig,
"#nu_{e} Signal",
"l");
407 if(
data) leg->AddEntry(hdata,
"ND Data",
"epl");
408 if(
tot)leg->AddEntry(htot,
"Total Bkg",
"l");
409 leg->AddEntry(hBeam,
"Beam #nu_{e}",
"l");
410 leg->AddEntry(hNC,
"NC",
"l");
411 leg->AddEntry(hCC,
"#nu_{#mu} CC",
"l");
417 std::vector<TString>
title,
418 double x1=kLegendDefaultX1,
double y1=kLegendDefaultY1,
419 double x2=kLegendDefaultX2,
double y2=kLegendDefaultY2)
422 leg->SetFillStyle(0);
423 for (
unsigned int v = 0;
v<h.size(); ++
v){
424 if (isdata[
v]) TLegendEntry*
le = leg->AddEntry(h[v],title[v],
"epl");
426 if (h[v]->GetLineColor() ==
kRed) TLegendEntry* le = leg->AddEntry(h[v],title[v],
"l");
427 else TLegendEntry* le = leg->AddEntry(h[v],title[v],
"f");
435 std::vector<TString> topOption)
437 TCanvas *cpid =
new TCanvas(
"cpid",
"cpid", 410, 850);
442 p[0] =
new TPad(
"",
"", 0, 0, 1, 1);
443 p[1] =
new TPad(
"",
"", 0, 0, 1, 1);
444 p[2] =
new TPad(
"",
"", 0, 0, 1, 1);
445 p[0]->SetTopMargin(0.1);
446 p[0]->SetBottomMargin(.634);
447 p[1]->SetTopMargin(.369);
448 p[1]->SetBottomMargin(.369);
449 p[2]->SetTopMargin(.634);
450 p[2]->SetBottomMargin(0.1);
451 p[0]->SetFillStyle(0);
452 p[1]->SetFillStyle(0);
453 p[2]->SetFillStyle(0);
460 for(
int k=pidBins-1; k>=0; --k){
461 he[k] = (TH1*) topHistos[0]->
Clone();
463 he[k]->GetXaxis()->SetRangeUser(Ebins*k,Ebins*(k+1));
464 if ( k>0 ) he[k]->GetXaxis()->SetLabelSize(0);
465 p[pidBins-1-k]->cd();
466 he[k]->Draw(topOption[0]);
467 he[k]->GetYaxis()->SetTitleSize(26);
468 he[k]->GetYaxis()->SetTitleFont(43);
469 he[k]->GetYaxis()->SetTitleOffset(1.25);
470 he[k]->GetYaxis()->SetLabelFont(43);
471 he[k]->GetYaxis()->SetLabelSize(14);
472 he[k]->GetXaxis()->SetLabelOffset(-0.008);
473 he[k]->GetXaxis()->SetTitleOffset(0.5);
474 p[pidBins-1-k]->SetGridy();
476 for(
unsigned int ii=0;ii<topHistos.size();ii++){
477 for(
int j=pidBins-1;
j>=0; --
j){
478 h[ii][
j] = (TH1*) topHistos[ii]->
Clone();
479 h[ii][
j]->GetXaxis()->SetRangeUser(Ebins*
j,Ebins*(j+1));
480 p[pidBins-1-
j]->cd();
481 h[ii][
j]->Draw(
"same "+topOption[ii]);
484 for(
int h=pidBins-1; h>=0; --
h){
486 p[pidBins-1-
h]->cd();
487 he[
h]->Draw(
"same "+topOption[0]);
496 std::vector<TString> topOption,
497 std::vector<TH1*> bottomHistos,
498 std::vector<TString> bottomOption,
502 TCanvas *
c =
new TCanvas(
"c",
"canvas", 600+(150*pidaxis), 700-(50*pidaxis));
503 auto h1 = (TH1*) topHistos[0]->
Clone();
504 auto h3 = (TH1*) bottomHistos[0]->
Clone();
509 TPad*
p1 =
new TPad(
"",
"", 0, 0, 1, 1);
510 TPad*
p2 =
new TPad(
"",
"", 0, 0, 1, 1);
511 p1->SetBottomMargin(.3);
512 p2->SetTopMargin(.7);
520 h1->Draw(topOption[0]);
521 for(
unsigned int ii=0;ii<topHistos.size();ii++){
522 topHistos[ii]->Draw(
"same "+topOption[ii]);}
523 h1->Draw(
"same "+topOption[0]);
526 h1->GetYaxis()->SetTitleSize(26);
527 h1->GetYaxis()->SetTitleFont(43);
528 h1->GetYaxis()->SetTitleOffset(1.2);
529 h1->GetYaxis()->SetLabelFont(43);
530 h1->GetYaxis()->SetLabelSize(18);
531 h1->GetXaxis()->SetLabelSize(0);
532 h1->GetXaxis()->SetTitleSize(0);
541 h3->GetYaxis()->SetTitleSize(26);
542 h3->GetYaxis()->SetTitleFont(43);
543 h3->GetYaxis()->SetTitleOffset(1.2);
544 h3->GetYaxis()->SetLabelFont(43);
545 h3->GetYaxis()->SetLabelSize(15);
546 h3->GetYaxis()->SetDecimals();
548 h3->GetXaxis()->SetTitle(
h1->GetXaxis()->GetTitle());
549 h3->GetXaxis()->SetTitleSize(26);
550 h3->GetXaxis()->SetTitleFont(43);
551 h3->GetXaxis()->SetLabelFont(43);
552 h3->GetXaxis()->SetLabelSize(18);
553 h3->GetXaxis()->SetLabelOffset(0.005);
554 h3->Draw(bottomOption[0]);
555 for(
unsigned int ii=0;ii<bottomHistos.size();ii++){
556 bottomHistos[ii]->Draw(
"same "+bottomOption[ii]);}
557 h3->Draw(
"same "+bottomOption[0]);
559 auto lone =
new TLine();
560 lone->SetLineStyle(3);
562 lone->DrawLine(p2->GetUxmin(),1.0,p2->GetUxmax(),1.0);
563 lone->DrawLine(p2->GetUxmin(),1.2,p2->GetUxmax(),1.2);
564 lone->DrawLine(p2->GetUxmin(),0.8,p2->GetUxmax(),0.8);
565 float p2Y = p2->GetUymax()-p2->GetUymin();
566 h3->GetYaxis()->SetRangeUser(p2->GetUymin()-(0.01*p2Y),p2->GetUymax()+(0.01*p2Y));
575 void PlotDataMC(std::vector<TH1*>vnom, std::vector<bool>isdata,
576 TLegend *
leg, TString pidtag=
"",TString htag=
"",
577 TString
out_name=
"plot_FD",
bool ratioplot=
false,
578 bool ratioerrors=
false,
bool pidaxis = 0)
580 gStyle->SetPadLeftMargin(0.15);
581 gStyle->SetPadTopMargin(0.08);
584 for (
unsigned int bin = 1;
bin<=vnom[0]->GetNbinsX(); ++
bin){
585 bool filledbin=
false;
586 for (
unsigned int v = 0;
v<vnom.size(); ++
v){
588 if (vnom[v]->GetBinContent(
bin) > 0 ) filledbin=
true;
590 for (
unsigned int v = 0;
v<vnom.size(); ++
v){
591 if (vnom[
v]->GetBinContent(
bin) == 0 ){
592 if ( !filledbin ) vnom[
v]->SetBinContent(
bin,-0.01);
593 else vnom[
v]->SetBinContent(
bin,0.001);
595 vnom[
v]->GetYaxis()->SetTitleOffset(1.05);
601 for (
unsigned int v = 0;
v<vnom.size(); ++
v)
602 if (vnom[
v]->GetLineColor() == kBlack ) datahist=
v;
606 auto c2 =
new TCanvas (
"c2",
"c2",715,500);
607 for (
unsigned int v = 0;
v<vnom.size(); ++
v){
610 if (v == datahist ) gr->Draw(
"ep same");
611 else vnom[
v]->Draw(
"p same");
613 else vnom[
v]->Draw(
"hist same");
627 TLatex* ltx =
new TLatex(.175, 0.74,
630 ltx->SetTextAlign(11);
631 ltx->SetTextSize(0.0365);
637 gPad->SetFillStyle(0);
650 std::vector<TH1*> vratio;
651 for(
unsigned int ii=1;ii<vnom.size();ii++){
653 htemp->Divide(vnom[datahist],vnom[ii]);
654 htemp->SetMarkerColor(vnom[ii]->GetLineColor());
655 htemp->SetMarkerStyle(34);
656 htemp->SetMarkerSize(1.2);
657 htemp->GetYaxis()->SetTitle(
"Data / MC");
658 htemp->GetYaxis()->SetRangeUser(0.5,1.5);
659 vratio.push_back(htemp);
661 std::vector<TString> topOption = {
" ",
"hist"};
662 std::vector<TString> bottomOption (vratio.size(),ratioerrors ?
"p":
"hist p");
663 auto cratio2=
RatioPlot(vnom, topOption,vratio,bottomOption,pidtag,pidaxis);
667 cratio2->SaveAs(
"plots/"+
out_name+
"_datamc_ratio.pdf");
668 cratio2->SaveAs(
"plots/png_"+
out_name+
"_datamc_ratio.png");
669 cratio2->SaveAs(
"plots/rootfiles/"+
out_name+
"_datamc_ratio.root");
675 TLegend *
leg,TString pidtag=
"",
677 bool ratioplot=
false,
bool pidaxis=
false)
680 std::vector<TString> h_opt;
681 for(
auto h:htots) h_opt.push_back(
"hist");
682 h_opt.push_back(
"p");
685 auto c1 =
new TCanvas (
"c1",
"c1");
686 htots[0]->SetMaximum(1.3*htots[0]->
GetMaximum());
687 for(
auto h:htots)
h->Draw(
"hist same");
688 htots[0]->Draw(
"hist same");
689 hdata->Draw(
"p same");
698 c1->SaveAs(
"plots/png_"+
out_name+
"_datamc.png");
699 c1->SaveAs(
"plots/rootfiles/"+
out_name+
"_datamc.root");
702 std::vector<TH1*> ratios;
703 std::vector<TString> r_opt;
705 auto ratio = (TH1*)
h->Clone();
707 ratio->GetYaxis()->SetRangeUser(0.5,1.5);
708 ratios.push_back(
ratio);
709 r_opt.push_back(
"e1");
711 ratios[0]->GetYaxis()->SetTitle(
"Data / MC");
712 ratios[0]->GetYaxis()->SetRangeUser(0.5,1.5);
713 htots.push_back(hdata);
714 auto cratio1 =
RatioPlot(htots,h_opt,ratios,r_opt,pidtag,pidaxis);
718 cratio1->SaveAs(
"plots/"+
out_name+
"_datamc_ratio.pdf");
719 cratio1->SaveAs(
"plots/png_"+
out_name+
"_datamc_ratio.png");
720 cratio1->SaveAs(
"plots/rootfiles/"+
out_name+
"_datamc_ratio.root");
723 htots.push_back(hdata);
737 TLegend *
leg, TString pidtag=
"",
738 TString
out_name=
"plot_nd",
bool ratioplot=
false,
739 bool ratioerrors=
false,
bool pidaxis = 0)
744 auto c2 =
new TCanvas (
"c2",
"c2");
745 for(
auto v:vnom)
v->Draw(
"hist same");
746 for(
auto v:vshi)
v->Draw(
"hist same");
755 c2->SaveAs(
"plots/png_"+
out_name+
"_mccomp.png");
756 c2->SaveAs(
"plots/rootfiles/"+
out_name+
"_mccomp.root");
760 std::vector<TH1*> vratio;
761 for(
unsigned int ii=0;ii<vshi.size();ii++){
763 htemp->Divide(vshi[ii],vnom[ii]);
764 htemp->SetMarkerColor(vshi[ii]->GetLineColor());
765 htemp->SetMarkerStyle(34);
766 htemp->SetMarkerSize(1.2);
767 htemp->GetYaxis()->SetTitle(
"Shift / Nom");
768 htemp->GetYaxis()->SetRangeUser(0.5,1.5);
769 vratio.push_back(htemp);
771 vnom.insert(vnom.end(),vshi.begin(),vshi.end());
772 std::vector<TString> topOption (vnom.size(),
"hist");
773 std::vector<TString> bottomOption (vratio.size(),ratioerrors ?
"p":
"hist p");
774 auto cratio2=
RatioPlot(vnom, topOption,vratio,bottomOption,pidtag,pidaxis);
778 cratio2->SaveAs(
"plots/"+
out_name+
"_mccomp_ratio.pdf");
779 cratio2->SaveAs(
"plots/png_"+
out_name+
"_mccomp_ratio.png");
780 cratio2->SaveAs(
"plots/rootfiles/"+
out_name+
"_mccomp_ratio.root");
783 vnom.insert(vnom.end(),vshi.begin(),vshi.end());
784 std::vector<TString> topOption (vnom.size(),
"hist");
797 void CompareNDDataMC (TDirectory * d_no, TDirectory* d_sh,TString plottitle,TString
out_name, TString pidtag,
bool printtable=
true,
bool pidaxis = 0){
799 auto kNDPOT =
data->POT();
800 auto hdata =
data->ToTH1(kNDPOT);
801 hdata ->SetMarkerStyle(kFullCircle);
808 auto htot_no = (TH1*)hnue_no->Clone();
809 htot_no->Add(hnumu_no);htot_no->Add(hnc_no);
816 auto htot_sh = (TH1*)hnue_sh->Clone();
817 htot_sh->Add(hnumu_sh);htot_sh->Add(hnc_sh);
822 htot_no->SetMaximum(
FindLimitY({htot_no,hdata}));
823 htot_no->SetTitle(plottitle);
824 hnumu_no->SetMaximum(
FindLimitY({hnumu_no,hnc_no,hnue_no,hnumu_sh,hnc_sh,hnue_sh}));
825 hnumu_no->SetTitle(plottitle);
829 TGaxis::SetMaxDigits(3);
830 htot_no->GetXaxis()->SetDecimals();
832 auto leg1=
new TLegend(0.63,0.67,0.85,0.87);
833 leg1->SetFillStyle(0);
834 leg1->AddEntry(hdata,
"ND Data",
"epl");
835 leg1->AddEntry(htot_no,
"Nominal",
"l");
836 leg1->AddEntry(htot_sh,
"Shift",
"l");
841 std::vector<TH1*>vnom={hnumu_no,hnc_no,hnue_no};
842 std::vector<TH1*>vshi={hnumu_sh,hnc_sh,hnue_sh};
845 TH1* htemp =
new TH1F(
"",
"", 1, 0, 1);
846 htemp->SetLineColor(kGray+1);
847 htemp->SetLineStyle(2);
848 leg2->AddEntry(htemp,
"Shifted",
"l");
854 std:: cout <<
"% - " << (plottitle.Remove(plottitle.First(
";"))).Data() <<
" - ND - " << pidtag <<
" - " <<
std::endl;
856 {htot_sh,hnue_sh,hnc_sh,hnumu_sh},
857 {
"Total MC",
"Beam $\\nu_{e}$",
"NC",
"$\\nu_{\\mu}$ CC"},
858 (plottitle.Remove(plottitle.First(
";")))+pidtag
902 hNueBeam_sh->SetLineStyle(2);
903 hNumu_sh->SetLineStyle(2);
904 hNC_sh->SetLineStyle(2);
905 hNueSig_sh->SetLineStyle(2);
907 hNueSig_no->SetMaximum(1.3*
max(hNueSig_no->GetMaximum(),hNueSig_sh->GetMaximum()));
908 hNueSig_no->SetMaximum(
FindLimitY({hNueSig_no,hNueBeam_no,hNC_no,hNumu_no,
909 hNueSig_sh,hNueBeam_sh,hNC_sh,hNumu_sh}));
910 hNueSig_no->SetTitle(plottitle);
912 TGaxis::SetMaxDigits(3);
914 std::vector<TH1*>vnom={hNueSig_no,hNueBeam_no,hNC_no,hNumu_no};
915 std::vector<TH1*>vshi={hNueSig_sh,hNueBeam_sh,hNC_sh,hNumu_sh};
917 TH1* htemp =
new TH1F(
"",
"", 1, 0, 1);
918 htemp->SetLineColor(kGray+1);
919 htemp->SetLineStyle(2);
920 leg->AddEntry(htemp,
"Shifted",
"l");
925 auto hTotBkg_no = (TH1*) hNueBeam_no->Clone();
926 hTotBkg_no->Add(hNC_no); hTotBkg_no->Add(hNumu_no); hTotBkg_no->Add(hTau_no);
928 auto hTotBkg_sh = (TH1*) hNueBeam_sh->Clone();
929 hTotBkg_sh->Add(hNC_sh); hTotBkg_sh->Add(hNumu_sh); hTotBkg_sh->Add(hTau_sh);
932 std::cout <<
"% - " << (plottitle.Remove(plottitle.First(
";"))).Data() <<
" - " << pidtag <<
" - " <<
std::endl;
933 ComparisonTable ({hNueSig_no, hTotBkg_no, hNueBeam_no,hNC_no,hNumu_no,hTau_no},
934 {hNueSig_sh, hTotBkg_sh, hNueBeam_sh,hNC_sh,hNumu_sh,hTau_sh},
935 {
"$\\nu_e$ signal",
"Total Beam Bkg",
"Beam $\\nu_{e}$CC",
"NC",
"$\\nu_{\\mu}$CC",
"$\\nu_{\\tau}$CC"},
936 (plottitle.Remove(plottitle.First(
";")))+pidtag
946 hplus ->SetLineColor(mccolor);
947 hminus->SetLineColor(mccolor);
948 hplus ->SetFillColor(mccolor);
949 hminus->SetFillColor(10);
950 hplus ->SetMarkerColor(
kRed-7); hminus->SetMarkerColor(
kBlue-7);
953 for(
int i=1;
i<=hplus->GetNbinsX();
i++){
954 double tplus = hplus->GetBinContent(
i);
955 double tminus = hminus->GetBinContent(
i);
956 hplus->SetBinContent(
i,
max(tplus, tminus));
957 hminus->SetBinContent(
i,
min(tplus,tminus));
964 TLegend *
leg, TString pidtag=
"",
965 TString
out_name=
"plot_nd",
bool ratioplot=
false,
966 bool ratioerrors=
false,
bool pidaxis = 0)
968 auto hplu_plot = (TH1*) vshi[0]->
Clone();
969 auto hmin_plot = (TH1*) vshi[1]->
Clone();
971 std::vector<TH1*>vshi_plot={hplu_plot,hmin_plot};
975 auto c2 =
new TCanvas (
"c2",
"c2");
976 vnom[0]->SetMaximum(1.2*vnom[0]->
GetMaximum());
977 vnom[0]->Draw(
"hist");
978 for(
auto v:vshi_plot)
v->Draw(
"hist same");
979 for(
auto v:vnom)
v->Draw(
"hist same");
988 c2->SaveAs(
"plots/png_"+
out_name+
"_mccomp.png");
989 c2->SaveAs(
"plots/rootfiles/"+
out_name+
"_mccomp.root");
993 std::vector<TH1*> vratio;
994 for(
unsigned int ii=0;ii<vshi.size();ii++){
996 htemp->Divide(vshi[ii],vnom[0]);
997 htemp->SetMarkerStyle(34);
998 htemp->SetMarkerSize(1.2);
999 htemp->GetYaxis()->SetTitle(
"Shift / Nom");
1000 vratio.push_back(htemp);
1002 vratio[0]->GetYaxis()->SetRangeUser(0.5,1.5);
1003 vnom.insert(vnom.begin()+1,vshi_plot.begin(),vshi_plot.end());
1004 std::vector<TString> topOption (vnom.size(),
"hist");
1005 std::vector<TString> bottomOption (vratio.size(),ratioerrors ?
"p":
"hist p");
1006 auto cratio2=
RatioPlot(vnom, topOption,vratio,bottomOption,pidtag,pidaxis);
1010 cratio2->SaveAs(
"plots/"+
out_name+
"_mccomp_ratio.pdf");
1011 cratio2->SaveAs(
"plots/png_"+
out_name+
"_mccomp_ratio.png");
1012 cratio2->SaveAs(
"plots/rootfiles/"+
out_name+
"_mccomp_ratio.root");
1029 ret.
bkg = (TH1*)ret.
beam->Clone();
1039 void CompareNDDataMC (TDirectory * d_no, TDirectory* d_pl,TDirectory*d_mi,TString plottitle,TString
out_name, TString pidtag,
bool printtable=
true,
bool pidaxis = 0){
1041 auto kNDPOT =
data->POT();
1042 auto hdata =
data->ToTH1(kNDPOT);
1043 hdata ->SetMarkerStyle(kFullCircle);
1051 hnom.bkg->SetMaximum(
FindLimitY({hnom.bkg,hdata,hplu.bkg,hmin.bkg}));
1052 hnom.bkg->SetTitle(plottitle);
1055 TGaxis::SetMaxDigits(3);
1056 hnom.bkg->GetXaxis()->SetDecimals();
1058 auto leg1=
new TLegend(kLegendDefaultX1,kLegendDefaultY1,kLegendDefaultX2,kLegendDefaultY2);
1059 leg1->SetFillStyle(0);
1060 leg1->AddEntry(hdata,
"ND Data",
"epl");
1061 leg1->AddEntry(hnom.bkg,
"Nominal",
"l");
1062 leg1->AddEntry(hplu.bkg,
"Syst.",
"f");
1069 leg2->AddEntry(hplu.bkg,
"Syst.",
"f");
1076 std::cout <<
"% - " << (plottitle.Remove(plottitle.First(
";"))).Data() <<
" - ND - " << pidtag <<
" - " <<
std::endl;
1078 {hplu.bkg,hplu.beam,hplu.nc,hplu.numu},
1079 {hmin.bkg,hmin.beam,hmin.nc,hmin.numu},
1080 {
"Total MC",
"Beam $\\nu_{e}$",
"NC",
"$\\nu_{\\mu}$ CC"},
1081 (plottitle.Remove(plottitle.First(
";")))+pidtag
1101 ret.
bkg = (TH1*) ret.
beam->Clone();
1117 for (
int i=0;
i<s.Length();
i++){
1133 TGaxis::SetMaxDigits(3);
1137 hnom.bkg->SetMaximum(
FindLimitY({hnom.bkg,hplu.bkg,hmin.bkg}));
1138 hnom.bkg->SetTitle(plottitle);
1141 leg->AddEntry(hplu.bkg,
"Syst.",
"f");
1148 hnom.sig->SetMaximum(
FindLimitY({hnom.sig,hplu.sig,hmin.sig}));
1149 hnom.sig->SetTitle(plottitle);
1151 auto leg2 =
new TLegend(kLegendDefaultX1,kLegendDefaultY1,kLegendDefaultX2,kLegendDefaultY2);
1152 leg2->AddEntry(hnom.sig,
"#nu_{e} Signal",
"l");
1153 leg2->AddEntry(hplu.sig,
"Syst.",
"f");
1160 std::cout <<
"% - " << (plottitle.Remove(plottitle.First(
";"))).Data() <<
" - " << pidtag <<
" - " <<
std::endl;
1161 ComparisonTable ({hnom.sig, hnom.bkg, hnom.beam, hnom.nc, hnom.numu, hnom.tau},
1162 {hplu.sig, hplu.bkg, hplu.beam, hplu.nc, hplu.numu, hplu.tau},
1163 {hmin.sig, hmin.bkg, hmin.beam, hmin.nc, hmin.numu, hmin.tau},
1164 {
"$\\nu_e$ signal",
"Total Beam Bkg",
"Beam $\\nu_{e}$CC",
"NC",
"$\\nu_{\\mu}$CC",
"$\\nu_{\\tau}$CC"},
1176 std::vector<TString>
channels = {
"numutonumucc" ,
"numutonuecc",
1177 "numutonutaucc",
"nuetonumucc",
1178 "nuetonuecc" ,
"nuetonutaucc" ,
1180 std::vector<TH1*> vshi = {hshi.numutonumu, hshi.sig,
1181 hshi.numutonutau, hshi.nuetonumu,
1182 hshi.beam, hshi.nuetonutau,
1184 std::vector<TH1*> vnom = {hnom.numutonumu, hnom.sig,
1185 hnom.numutonutau, hnom.nuetonumu,
1186 hnom.beam, hnom.nuetonutau,
1188 if(hshi.sig->GetNbinsX()!=30) {
1190 <<
"MakeNueSystematicsFile doesn't know what to do with this binning\n" 1191 <<
"Nothing will be saved\n";
1202 for(
unsigned int i=0;
i<vshi.size();
i++){
1203 vshi[
i]->Divide(vnom[
i]);
1204 TString hname =
"fdratio_"+pid+
"_"+channels[
i]+
"_"+
sigma;
1205 TString htitle =
"FD "+channels[
i]+
";Reconstructed Energy (GeV);"+
pid;
1207 TH2D* h_epid =
new TH2D(hname,htitle,3,pidEdges,10,0,5);
1208 for(
int xx=1;
xx<=3;
xx++){
1209 for(
int yy=1;yy<=10;yy++){
1210 double ratio = vshi[
i]->GetBinContent(yy+10*(
xx-1));
1212 if(ratio>0)h_epid->SetBinContent(
xx,yy,ratio-1);
1215 nue_syst_file->cd();
1223 std::vector<int> ndcounts;
1224 std::vector<int> fdcounts;
1225 std::vector<TString>
labels;
1227 labels.push_back(
"Data - $\\nu_\\mu $");
1228 ndcounts.push_back(((TH1*)dpred->Get(
"extrap/MEextrap/Decomp/data_comp/hist"))->
GetEntries());
1229 fdcounts.push_back(0);
1230 labels.push_back(
"Data - $\\nu_e $");
1231 ndcounts.push_back(((TH1*)dpred->Get(
"extrap/EEextrap/Decomp/data_comp/hist"))->
GetEntries());
1232 fdcounts.push_back(0);
1235 labels.push_back(
"Signal - $\\nu_\\mu \\to \\nu_e $");
1236 labels.push_back(
"Signal - $\\bar{\\nu}_\\mu \\to \\bar{\\nu}_e $");
1237 for(TString xp:{
"ME",
"MEAnti"}){
1238 ndcounts.push_back(((TH1*)dpred->Get(
"extrap/"+xp+
"extrap/RecoToTrueND/hist"))->
GetEntries());
1239 fdcounts.push_back(((TH1*)dpred->Get(
"extrap/"+xp+
"extrap/TrueToRecoFD/hist"))->
GetEntries());
1243 labels.push_back(
"Bkg. - $\\nu_e \\to \\nu_e $");
1244 labels.push_back(
"Bkg. - $\\nu_\\mu \\to \\nu_\\mu $");
1245 labels.push_back(
"Bkg. - $NC \\to NC $");
1246 for(TString xp:{
"EE",
"MM",
"NC"}){
1247 ndcounts.push_back(((TH1*)dpred->Get(
"extrap/"+xp+
"extrap/RecoND/hist"))->
GetEntries());
1248 fdcounts.push_back(((TH1*)dpred->Get(
"extrap/"+xp+
"extrap/TrueToRecoFD/hist"))->
GetEntries());
1252 labels.push_back(
"Bkg. - $\\nu_e \\to \\nu_\\tau $");
1253 labels.push_back(
"Bkg. - $\\nu_\\mu \\to \\nu_\\tau $");
1254 for(TString xp: {
"ET",
"MT"}){
1255 ndcounts.push_back(0);
1256 fdcounts.push_back(((TH1*)dpred->Get(
"extrap/"+xp+
"extrap/RecoFD/hist"))->
GetEntries());
1258 labels.push_back(
"Bkg. - Other");
1259 double fdcount_minor=0;
1260 for(TString xp: {
"EM",
"EEAnti",
"EMAnti",
"ETAnti",
"MMAnti",
"MTAnti"}){
1261 fdcount_minor+=(((TH1*)dpred->Get(
"extrap/"+xp+
"extrap/RecoFD/hist"))->
GetEntries());
1263 ndcounts.push_back(0);
1264 fdcounts.push_back(fdcount_minor);
1269 std::cout <<
"\\begin{tabular}{ l | r | r }\n" 1270 << std::setw(45) <<
" " <<
" & " 1271 << std::setw(10) <<
"ND" <<
" & " 1272 << std::setw(10) <<
"FD" 1273 <<
" \\\\ \\hline \n";
1274 for(
unsigned int i=0;
i<labels.size();
i++){
1275 std::cout << std::setw(45) << labels[
i] <<
" & " 1276 << std::setw(10)<< ndcounts[
i] <<
" & " 1277 << std::setw(10)<< fdcounts[
i] <<
" \\\\ ";
struct DataMCComponets GetFDMCComponents(osc::IOscCalc *calc, IPrediction *pred_no, TString output_name="nue", int linestyle=1, bool bkgdetails=false)
struct NueMCComponents GetNDMCComponents(TDirectory *d_no, double kNDPOT, int linestyle=1)
double FindLimitY(std::vector< TH1 * > histos, bool doMax=true)
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
TString Latexify(TString s)
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
void PIDBinLabels(std::string pid)
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.
const Binning kCVNNLBinning
Float_t x1[n_points_granero]
void PlotMCComponentsComparison(std::vector< TH1 * >vnom, std::vector< TH1 * >vshi, TLegend *leg, TString pidtag="", TString out_name="plot_nd", bool ratioplot=false, bool ratioerrors=false, bool pidaxis=0)
static const double kFDPOT
void PIDBinLabels(bool shortaxis)
const double kLegendDefaultX2
void CenterTitles(TH1 *histo)
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
const Color_t kTotalMCErrorBandColor
const Color_t kNumuBackgroundColor
void PlotMCComponentsErrorBand(std::vector< TH1 * >vnom, std::vector< TH1 * >vshi, TLegend *leg, TString pidtag="", TString out_name="plot_nd", bool ratioplot=false, bool ratioerrors=false, bool pidaxis=0)
void PIDBinLabelsShortAxis()
const Color_t kNueSignalColor
void PlotNDDataTotalMCComparison(TH1 *hdata, std::vector< TH1 * > htots, TLegend *leg, TString pidtag="", TString out_name="plot_nd", bool ratioplot=false, bool pidaxis=false)
const XML_Char const XML_Char * data
void MakeNueSystematicsFile(IPrediction *pred_no, IPrediction *pred_sh, TString pid, TString sigma, TFile *nue_syst_file)
Charged-current interactions.
std::map< ToFCounter, std::vector< unsigned int > > channels
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
void ComparisonTable(const std::vector< TH1 * > &mcnom, const std::vector< TH1 * > &mcplus, const std::vector< TH1 * > &mcminus, std::vector< TString > labels, TString ltxcommand="")
void CompareNDDataMC(TDirectory *d_no, TDirectory *d_sh, TString plottitle, TString out_name, TString pidtag, bool printtable=true, bool pidaxis=0)
void FixLegend(TLegend *leg, TString opt="default")
const Color_t kBeamNueBackgroundColor
void PIDDivisionsShortAxis()
const double kLegendDefaultY1
const double kLegendDefaultX1
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it's s...
TH1 * RebinToShortAxis(TH1 *h)
void PID2DShortAxis(TH1 *hist)
void FormatErrorBand(TH1 *hplus, TH1 *hminus, bool signal=false, bool fixbins=false)
const std::vector< double > & Edges() const
TString MakeLatexCommandName(TString mystring)
const bool printratioplots
static float min(const float a, const float b, const float c)
const double kLegendDefaultY2
TString MakeLatexCommandName(TString mystring)
TLegend * CustomLegend(std::vector< TH1 * >h, std::vector< bool >isdata, std::vector< TString >title, double x1=kLegendDefaultX1, double y1=kLegendDefaultY1, double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
const Binning kLIDNLBinning
const Binning kLEMNLBinning
Neutral-current interactions.
const double kSecondAnaRunMatchedMCPOTScale
assert(nhit_max >=nhit_nbins)
void PlotDataMC(std::vector< TH1 * >vnom, std::vector< bool >isdata, TLegend *leg, TString pidtag="", TString htag="", TString out_name="plot_FD", bool ratioplot=false, bool ratioerrors=false, bool pidaxis=0)
Both neutrinos and antineutrinos.
Standard interface to all prediction techniques.
void CenterTitles(TH1 *histo)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
const Color_t kTotalMCColor
void PID2DAxis(TH1 *axes, bool third=false)
All neutrinos, any flavor.
TLegend * DefaultNueLegend(bool sig=true, bool tot=true, bool data=false, double x1=kLegendDefaultX1, double y1=kLegendDefaultY1, double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
void ComparePredictions(IPrediction *pred_no, IPrediction *pred_sh, TString plottitle, TString out_name, TString pidtag, bool printtable=true, bool pidaxis=false)
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)
std::string UniqueName()
Return a different string each time, for creating histograms.
void PID2DAxis(TH1 *hist, bool shortaxis)
void PrintRawEventCounts(TDirectory *dpred, TString title)
TCanvas * ExPIDPlot(std::vector< TH1 * > topHistos, std::vector< TString > topOption)
void PIDDivisions(bool shortaxis)