draw_plots_util.h
Go to the documentation of this file.
6 #include "CAFAna/Core/Binning.h"
7 #include "CAFAna/Core/HistAxis.h"
16 #include "CAFAna/Vars/FitVars.h"
17 #include "CAFAna/Decomp/IDecomp.h"
19 #include "CAFAna/Extrap/IExtrap.h"
22 
23 #include "OscLib/OscCalcDumb.h"
24 #include "OscLib/IOscCalc.h"
25 
26 #include "TStyle.h"
27 #include "TFile.h"
28 #include "TH1.h"
29 #include "TCanvas.h"
30 #include "TGraph.h"
31 #include "TGraphAsymmErrors.h"
32 #include "TGaxis.h"
33 #include "TLatex.h"
34 #include "TPaveText.h"
35 #include "TLegend.h"
36 #include "TLine.h"
37 #include "TSystem.h"
38 #include "TError.h"
39 
40 #include <iostream>
41 #include <fstream>
42 #include <iomanip>
43 #include <cxxabi.h>
44 #include <string.h>
45 
46 #include "tex_tables_util.h"
47 
48 static bool isFHC=true; // Global variable for FHC/RHC
49 
50 //static const double kFDPOT= ana::kAna2017POT;
51 const bool printratioplots=true;
52 
53 //enum AxisType {kNue1bin, kNue3bin, kNueNDbin, kNumuOpt, kNumuOther};
54 
55 bool usedumb = false;
56 bool mergePeripheral=false;
57 
58 osc::IOscCalc * GetCalculator (bool usingdumb){
59  double dCP_Pi = .0;
60  double ssth23 = 0.56;
61  double dmsq32 = 2.48e-3;
62 
64  if(usingdumb){
65  osc::OscCalcDumb dumb;
66  calc = dumb.Copy();
67  }
68  else{
69  auto calc2 = DefaultOscCalc();
70  kFitDeltaInPiUnits.SetValue(calc2, dCP_Pi);
71  kFitSinSqTheta23.SetValue(calc2, ssth23);
72  kFitDmSq32.SetValue(calc2, dmsq32);
73  calc = calc2->Copy();
74  }
75  return calc;
76 }
77 
78 
79 //////////////////////////////////////////////////////////////////////
80 // Histogram
81 //////////////////////////////////////////////////////////////////////
82 
83 TH1* PrettyRatio(TH1* shi, TH1* nom, int color,
84  double alpha=1., TString titley="1#sigma Shift / Nom")
85 {
86  auto htemp = reinterpret_cast<TH1*>(nom->Clone(ana::UniqueName().c_str()));
87  htemp->Divide(shi,nom);
88  htemp->SetMarkerColorAlpha(color,alpha);
89  htemp->SetLineColorAlpha(color,alpha);
90  htemp->SetMarkerStyle(34);
91  htemp->SetMarkerSize(1.4);
92  htemp->GetYaxis()->SetTitle(titley);
93  return htemp;
94 }
95 
96 //////////////////////////////////////////////////////////////////////
97 
98 double FindLimitY (std::vector<TH1*> histos, bool doMax = true)
99 {
100  double max = histos[0]->GetMaximum();
101  double min = histos[0]->GetMinimum();
102 
103  for(unsigned int i=0; i<histos.size(); i++){
104  if ( histos[i]->GetMaximum() > max ) max = histos[i]->GetMaximum();
105  if ( histos[i]->GetMinimum() < min ) min = histos[i]->GetMinimum();
106  }
107 
108  if ( doMax ) return max;
109  else return min;
110 }
111 
112 //////////////////////////////////////////////////////////////////////
113 void FormatErrorBand(TH1* hnom, TH1* hplus, TH1* hminus,
114  TString signal="bkg", bool fixbins=false)
115 {
116  int mccolor;
117  // why oh why the switch does not work for strings
118  if(signal == "sig") mccolor = kViolet -9; // signal color
119  else if(signal == "bkg") mccolor = kTotalMCErrorBandColor; // total bkg color
120  else if(signal == "beam") mccolor = kBeamNueBackgroundColor -10; // beam bkg color
121  else if(signal == "numu") mccolor = kNumuBackgroundColor -10; // numu bkg color
122  else if(signal == "nc") mccolor = kNCBackgroundColor -9; // nc bkg color
123  else mccolor = kTotalMCErrorBandColor;
124 
125  std::cout << "hplus: " << hplus->Integral() << std::endl;
126  std::cout << "hminus: " << hminus->Integral() << std::endl;
127  hplus ->SetLineColor(mccolor);
128  hminus->SetLineColor(mccolor);
129  hplus ->SetFillColor(mccolor);
130  hminus->SetFillColor(10);
131  hplus ->SetMarkerColor(kRed-7);
132  hminus->SetMarkerColor(kBlue-7);
133 
134  if (fixbins)
135  {
136  for(int i=1; i <= hplus->GetNbinsX(); i++)
137  {
138  double tnom = hnom->GetBinContent(i);
139  double tplus = hplus->GetBinContent(i);
140  double tminus = hminus->GetBinContent(i);
141  hplus->SetBinContent(i,max(tplus, max(tnom, tminus)));
142  hminus->SetBinContent(i,min(tplus, min(tnom, tminus)));
143  //hplus->SetBinContent(i,max(tplus, tminus));
144  //hminus->SetBinContent(i,min(tplus, tminus));
145  }
146  }
147 }
148 
149 //////////////////////////////////////////////////////////////////////
150 TH1* HistNumuOpt (TH1* orig){
151  auto ret = (TH1*) orig->Clone(UniqueName().c_str());
152  ret->Scale(0.1,"width");
153  ret->GetYaxis()->SetTitle("Events / 0.1 GeV");
154  ret->SetBit(kCanDelete);
155  ret->SetDirectory(gROOT);
156  return ret;
157 }
158 
159 //////////////////////////////////////////////////////////////////////
160 // Adornments
161 //////////////////////////////////////////////////////////////////////
162 
163 void PrettyTag(TString pid, int color, double xndc, double yndc){
164  TLatex* tag = new TLatex(xndc, yndc, pid);
165  tag->SetTextColor(color);
166  tag->SetNDC();
167  tag->SetTextSize(0.04);
168  tag->SetTextFont(42);
169  tag->SetTextAlign(22);
170  tag->Draw();
171  tag->SetBit(kCanDelete);
172 }
173 
174 void PIDTag(TString pid) { PrettyTag (pid, kAzure+3, .83, .92); }
175 void HTag (TString hie) { PrettyTag (hie, (hie=="NH")?kAzure+3:kOrange+2, .18, .825);}
176 
177 //////////////////////////////////////////////////////////////////////
178 // Legends
179 //////////////////////////////////////////////////////////////////////
180 
181 const double kLegendDefaultX1 = 0.63;
182 const double kLegendDefaultX2 = 0.85;
183 const double kLegendDefaultY1 = 0.68;
184 const double kLegendDefaultY2 = 0.88;
185 
186 TLegend * CustomLegend(
187  std::vector<TH1*>h, std::vector<TString> titles, std::vector<TString> option, int n=1,
188  double x1=kLegendDefaultX1, double y1=kLegendDefaultY1,
189  double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
190 {
191  TLegend* leg = new TLegend(x1, y1, x2, y2);
192  leg->SetFillStyle(0);
193  leg->SetNColumns(n);
194  for (unsigned int v = 0; v < h.size(); ++v)
195  leg->AddEntry(h[v], titles[v], option[v]);
196  leg->SetBit(kCanDelete);
197  return leg;
198 }
199 
200 //////////////////////////////////////////////////////////////////////
201 
202 void FixLegend(TLegend * leg, TString opt="default")
203 {
204  leg->SetFillStyle(0);
205  if(opt=="default")
206  {
207  leg->SetX1(kLegendDefaultX1); leg->SetX2(kLegendDefaultX2);
208  leg->SetY1(kLegendDefaultY1); leg->SetY2(kLegendDefaultY2);
209  }
210  if(opt=="ExPID")
211  {
212  leg->SetX1(0.355); leg->SetX2(0.515);
213  leg->SetY1(0.65); leg->SetY2(0.82);
214  leg->SetBorderSize(0);
215  leg->SetTextSize(0.03);
216  }
217  if(opt=="ExPIDNearDet")
218  {
219  leg->SetX1(0.695); leg->SetX2(0.855);
220  leg->SetY1(0.65); leg->SetY2(0.82);
221  leg->SetBorderSize(0);
222  leg->SetTextSize(0.03);
223  }
224  if(opt=="ExPIDPads")
225  {
226  leg->SetX1(0.55); leg->SetX2(0.89);
227  leg->SetY1(0.75); leg->SetY2(0.89);
228  leg->SetTextSize(0.045);
229  }
230  if(opt=="ExRatio")
231  {
232  leg->SetTextSize(0.035);
233  leg->SetX1(0.1); leg->SetX2(0.3);
234  leg->SetY1(0.05); leg->SetY2(0.35);
235  }
236 }
237 //////////////////////////////////////////////////////////////////////
238 
239 TGraph * DummyGraph(int color=kBlack, int style =1, int linewidth=2, int markerstyle = kFullCircle)
240 {
241  auto gr = new TGraph();
242  gr->SetLineColor(color);
243  gr->SetLineStyle(style);
244  gr->SetLineWidth(linewidth);
245  gr->SetMarkerStyle(markerstyle);
246  gr->SetMarkerColor(color);
247  return gr;
248 }
249 //////////////////////////////////////////////////////////////////////
250 
252  bool sig=true, bool tot=true, bool data=false,
253  double x1=kLegendDefaultX1, double y1=kLegendDefaultY1,
254  double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
255 {
256  TLegend* leg = new TLegend(x1,y1,x2,y2);
257  leg->SetFillStyle(0);
258  if(sig) leg->AddEntry( DummyGraph(kNueSignalColor), "#nu_{e} Signal","l");
259  if(data) leg->AddEntry( DummyGraph(), "ND Data","epl");
260  if(tot) leg->AddEntry( DummyGraph(kTotalMCColor), "Total Bkg","l");
261  leg->AddEntry( DummyGraph(kBeamNueBackgroundColor), "Beam #nu_{e}","l");
262  leg->AddEntry( DummyGraph(kNCBackgroundColor), "NC","l");
263  leg->AddEntry( DummyGraph(kNumuBackgroundColor), "#nu_{#mu} CC","l");
264  leg->SetBit(kCanDelete);
265  return leg;
266 
267 }
268 //////////////////////////////////////////////////////////////////////
269 
271  bool sig=true, bool tot=true, bool data=false,
272  double x1=kLegendDefaultX1, double y1=kLegendDefaultY1,
273  double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
274 {
275  TLegend* leg = new TLegend(x1,y1,x2,y2);
276  leg->SetFillStyle(0);
277  if(sig)
278  leg->AddEntry(DummyGraph(kNueSignalColor), "#nu_{#mu} Signal", "l");
279  if(data)
280  leg->AddEntry(DummyGraph(), "ND Data", "epl");
281  if(tot)
282  leg->AddEntry(DummyGraph(kTotalMCColor), "Total Bkg", "l");
283  if(isNumuAna)
284  leg->AddEntry(DummyGraph(kNumuBackgroundColor), "#nu_{#mu} app", "l");
285  else
286  leg->AddEntry(DummyGraph(kNumuBackgroundColor), "#nu_{#mu} CC", "l");
287  leg->AddEntry(DummyGraph(kNCBackgroundColor), "NC","l");
288  leg->AddEntry(DummyGraph(kBeamNueBackgroundColor), "#nu_{e} CC","l");
289  leg->AddEntry(DummyGraph(kNuTauBackgroundColor), "#nu_{#tau} CC","l");
290  leg->SetBit(kCanDelete);
291  return leg;
292 }
293 //////////////////////////////////////////////////////////////////////
294 
295 void PrintCanvasAll (TCanvas * c, const TString outname, const TString outdir = "plots/",
296  const std::vector <TString> exts = {"pdf","png","root"} )
297 {
298  for(auto const & ext:exts){
299  gSystem -> MakeDirectory( outdir );
300  gSystem -> MakeDirectory( outdir+"rootfiles/" );
301  if(ext=="png") c->SaveAs(outdir+"png_"+outname + "." + ext);
302  else if(ext=="root") c->SaveAs(outdir+"rootfiles/"+outname + "." + ext);
303  else c->SaveAs(outdir+outname + "." + ext);
304  }
305  delete c;
306 }
307 
308 //////////////////////////////////////////////////////////////////////
309 // Nue 3 bins functions
310 //////////////////////////////////////////////////////////////////////
311 
312 // NB, loses all the style info, so call early
313 /*TH1* RebinToShortAxis(TH1* h)
314 {
315  TH1* ret = new TH1D(UniqueName().c_str(), h->GetTitle(),
316  18, 0, 18);
317  ret->GetXaxis()->SetTitle(h->GetXaxis()->GetTitle());
318  ret->GetYaxis()->SetTitle(h->GetYaxis()->GetTitle());
319 
320  std::vector<std::pair<int, int>> binMap = {{3, 2},
321  {4, 3},
322  {5, 4},
323  {6, 5}, // end of 1st spectrum
324  {13, 8},
325  {14, 9},
326  {15, 10},
327  {16, 11}, // end of 2nd
328  {23, 14},
329  {24, 15},
330  {25, 16},
331  {26, 17}};
332 
333  for(auto it: binMap) ret->SetBinContent(it.second, h->GetBinContent(it.first));
334 
335  return ret;
336 }
337 */
338 //////////////////////////////////////////////////////////////////////
340  TH1* h, bool coreOnlyND=false, bool binLabels=true,
341  const double textSize=0.04, const int textColor=kGray+3,
342  const int lineColor=kBlack, const int lineStyle=3, const double yNDC=0.96)
343 {
344  if (coreOnlyND)
345  Info("NuePID3bin", "specified core only/ND.");
346 
347  Nue2018ThreeBinDivisions(coreOnlyND, lineColor, lineStyle);
348  if (binLabels)
349  {
350  Nue2018ThreeBinLabels(yNDC, textSize, textColor, coreOnlyND);
351  }
352  else
353  {
354  if (textSize != 0.04)
355  Warning("NuePID3bin",
356  "binLabels == false --> specified text size has no effect");
357  if (textColor != kGray+3)
358  Warning("NuePID3bin",
359  "binLabels == false --> specified text color has no effect");
360  if (yNDC != 0.96)
361  Warning("NuePID3bin",
362  "binLabels == false --> specified yNDC has no effect");
363  }
364 }
365 //////////////////////////////////////////////////////////////////////
366 // TCanvas * ExPIDPlot (std::vector<TH1*> topHistos,
367 // std::vector<TString> topOption)
368 // {
369 // TCanvas *cpid = new TCanvas("cpid", "cpid", 410, 850);
370 // float Ebins = 10;
371 // int pidBins = 3;
372 //
373 // TPad* p[pidBins];
374 // p[0] = new TPad("", "", 0, 0, 1, 1);
375 // p[1] = new TPad("", "", 0, 0, 1, 1);
376 // p[2] = new TPad("", "", 0, 0, 1, 1);
377 // p[0]->SetTopMargin(0.1);
378 // p[0]->SetBottomMargin(.634);
379 // p[1]->SetTopMargin(.369);
380 // p[1]->SetBottomMargin(.369);
381 // p[2]->SetTopMargin(.634);
382 // p[2]->SetBottomMargin(0.1);
383 // p[0]->SetFillStyle(0);
384 // p[1]->SetFillStyle(0);
385 // p[2]->SetFillStyle(0);
386 // p[0]->Draw();
387 // p[1]->Draw();
388 // p[2]->Draw();
389 //
390 // TH1* he[pidBins];
391 // TH1* h[topHistos.size()][pidBins];
392 // for(int k=pidBins-1; k>=0; --k){
393 // he[k] = (TH1*) topHistos[0]->Clone();
394 // he[k]->GetYaxis()->SetRangeUser(FindLimitY(topHistos,0) ,1.2*FindLimitY(topHistos,1));
395 // he[k]->GetXaxis()->SetRangeUser(Ebins*k,Ebins*(k+1));
396 // if ( k>0 ) he[k]->GetXaxis()->SetLabelSize(0);
397 // p[pidBins-1-k]->cd();
398 // he[k]->Draw(topOption[0]);
399 // he[k]->GetYaxis()->SetTitleSize(24);
400 // he[k]->GetYaxis()->SetTitleFont(43);
401 // he[k]->GetYaxis()->SetTitleOffset(1.25);
402 // he[k]->GetYaxis()->SetLabelFont(43);
403 // he[k]->GetYaxis()->SetLabelSize(14);
404 // he[k]->GetXaxis()->SetLabelOffset(-0.008);
405 // he[k]->GetXaxis()->SetTitleOffset(0.5);
406 // p[pidBins-1-k]->SetGridy();
407 // }
408 // for(unsigned int ii=0;ii<topHistos.size();ii++){
409 // for(int j=pidBins-1; j>=0; --j){
410 // h[ii][j] = (TH1*) topHistos[ii]->Clone();
411 // h[ii][j]->GetXaxis()->SetRangeUser(Ebins*j,Ebins*(j+1));
412 // p[pidBins-1-j]->cd();
413 // h[ii][j]->Draw("same "+topOption[ii]);
414 // }
415 // }
416 // for(int h=pidBins-1; h>=0; --h){
417 // he[h]->SetStats(0);
418 // p[pidBins-1-h]->cd();
419 // he[h]->Draw("same "+topOption[0]);
420 // }
421 // p[0]->RedrawAxis(); //needed for white histograms e.g. beam
422 // p[1]->RedrawAxis();
423 // p[2]->RedrawAxis();
424 // return cpid;
425 // }
426 
427 //////////////////////////////////////////////////////////////////////
428 // Pretty Canvas
429 //////////////////////////////////////////////////////////////////////
430 TCanvas * RatioPlot(
431  std::vector<TH1*> topHistos, std::vector<TString> topOption,
432  std::vector<TH1*> bottomHistos, std::vector<TString> bottomOption,
433  TString pidtag, AxisType pidaxis=kNue1bin, const bool nearDet=false)
434 {
435  gStyle->SetTitleStyle(0);
436  TCanvas *c = new TCanvas("c", "canvas", 600, 700);
437  // if(pidaxis == kNue3bin) c->SetCanvasSize(750,650);
438  //define top and bottom pads
439  TPad* p1 = new TPad("", "", 0, 0, 1, 1);
440  TPad* p2 = new TPad("", "", 0, 0, 1, 1);
441  p1->SetBottomMargin(.3);
442  p2->SetTopMargin(.7);
443 
444  for(auto p:{p1,p2}){
445  p->SetFillStyle(0);
446  p->Draw();
447  }
448 
449  // format histograms
450  auto h1 = (TH1*) topHistos[0]->Clone();
451  auto h3 = (TH1*) bottomHistos[0]->Clone();
452 
453  for(auto & h:{h1,h3}){
454  h->SetStats(0);
455 
456  h->GetYaxis()->SetTitleSize(26);
457  h->GetYaxis()->SetTitleFont(43);
458  h->GetYaxis()->SetTitleOffset(1.3);
459  h->GetYaxis()->SetLabelFont(43);
460 
461  h->GetXaxis()->SetTitleOffset(1.2);
462  h->GetXaxis()->SetTitleSize(28);
463  h->GetXaxis()->SetTitleFont(43);
464  h->GetXaxis()->SetLabelFont(43);
465  }
466 
467  if (pidaxis == kNue3bin)
468  h1->GetYaxis()->SetRangeUser(
469  FindLimitY(topHistos,0), 1.5*FindLimitY(topHistos, 1));
470  else
471  h1->GetYaxis()->SetRangeUser(
472  FindLimitY(topHistos,0), 1.2*FindLimitY(topHistos, 1));
473  h1->GetYaxis()->SetLabelSize(18);
474 
475  h1->GetXaxis()->SetLabelSize(0);
476  h1->GetXaxis()->SetTitleSize(0);
477 
478  h3->SetTitle("");
479  h3->GetYaxis()->SetRangeUser(max(0.,2-FindLimitY(bottomHistos,1)),
480  1.01*FindLimitY(bottomHistos,1));
481  h3->GetYaxis()->SetLabelSize(15);
482  h3->GetYaxis()->SetDecimals();
483 
484  h3->GetXaxis()->SetTitle(h1->GetXaxis()->GetTitle());
485  h3->GetXaxis()->SetLabelSize(18);
486  h3->GetXaxis()->SetLabelOffset(0.005);
487 
488  p1->cd();
489 
490  h1->Draw(topOption[0]);
491 
492  for(unsigned int ii = 1; ii < topHistos.size(); ++ii) {
493  topHistos[ii]->Draw("same "+topOption[ii]);}
494  h1->Draw("same "+topOption[0]);
495 
496  p1->RedrawAxis(); //white histograms cover axes
497 
498  if(pidaxis==kNue3bin)
499  NuePID3bin(h1, nearDet);
500 
501  p2->cd();
502  h3->Draw(bottomOption[0]);
503  for(unsigned int ii = 1; ii < bottomHistos.size(); ++ii){
504  bottomHistos[ii]->Draw("same " + bottomOption[ii]);}
505  h3->Draw("same " + bottomOption[0]);
506  if(bottomHistos.size() == 2 && strcmp(h3->GetYaxis()->GetTitle(),"Data / MC") != 0)
507  {
508  TLegend * leg_ratio = CustomLegend(bottomHistos,{" +1#sigma"," -1#sigma"},{bottomOption[0],bottomOption[0]});
509  FixLegend(leg_ratio, "ExRatio");
510  leg_ratio -> Draw();
511  }
512  auto lone = new TLine();
513  lone->SetLineStyle(3);
514  lone->SetLineColor(kGray+3);
515  p2->Update();
516  lone->DrawLine(p2->GetUxmin(),1.0,p2->GetUxmax(),1.0);
517  lone->DrawLine(p2->GetUxmin(),1.2,p2->GetUxmax(),1.2);
518  lone->DrawLine(p2->GetUxmin(),0.8,p2->GetUxmax(),0.8);
519  float p2Y = p2->GetUymax()-p2->GetUymin();
520  h3->GetYaxis()->SetRangeUser(p2->GetUymin()-(0.01*p2Y),p2->GetUymax()+(0.01*p2Y));
521  p2->RedrawAxis();
522  if(pidaxis==kNue3bin)
523  NuePID3bin(h3, nearDet, false); // hist, coreOnlyND, binLabels
524 
525  return c;
526 }
527 //////////////////////////////////////////////////////////////////////
528 // Pretty plot styles: Data MC, MC components, error bands
529 //////////////////////////////////////////////////////////////////////
530 
531 /*void PlotDataMC(std::vector<TH1*>vnom, std::vector<bool>isdata,
532  TLegend * leg, TString pidtag="",TString htag="",
533  TString out_name="plot_FD",bool ratioplot=false,
534  bool ratioerrors=false, AxisType pidaxis = kNue1bin)
535 {
536  gStyle->SetPadLeftMargin(0.15);
537  gStyle->SetPadTopMargin(0.08);
538 
539  // Fill empty bins with 0 if in range and <0 if out of range
540  for ( int bin = 1; bin<=vnom[0]->GetNbinsX(); ++bin){
541  bool filledbin=false;
542  for ( unsigned int v = 0; v<vnom.size(); ++v){
543  if ( !isdata[v] )
544  if (vnom[v]->GetBinContent(bin) > 0 ) filledbin=true;
545  }
546  for ( uint v = 0; v<vnom.size(); ++v){
547  if (vnom[v]->GetBinContent(bin) == 0 ){
548  if ( !filledbin ) vnom[v]->SetBinContent(bin,-0.01);
549  else vnom[v]->SetBinContent(bin,0.001);
550  }
551  vnom[v]->GetYaxis()->SetTitleOffset(1.05);
552  }
553  }
554 
555  // Draw marker w/error or line when appropriate
556  int datahist=-1;
557  for ( unsigned int v = 0; v<vnom.size(); ++v)
558  if (vnom[v]->GetLineColor() == kBlack ) datahist=v;
559 
560  // Option to make a ratio plot given
561  if(!ratioplot){
562  auto c2 = new TCanvas ("c2","c2",715,500);
563  for ( unsigned int v = 0; v<vnom.size(); ++v){
564  TGraphAsymmErrors* gr = GraphWithPoissonErrors(vnom[v], 0, 1);
565  if ( isdata[v] ){
566  if (v == datahist ) gr->Draw("ep same");
567  else vnom[v]->Draw("p same");
568  }
569  else vnom[v]->Draw("hist same");
570  }
571  gPad->RedrawAxis();
572 
573  // Adjust for ExPID axis
574  if(pidaxis==kNue3bin)
575  {
576  FixLegend(leg, "ExPID");
577  Nue2018ThreeBinDivisions();
578  Nue2018ThreeBinLabels();
579  }
580 
581  // pot label
582  PrettyTag(TString::Format("%g#times10^{20} POT equiv.", 8.85).Data(),
583  kBlack, .175, .74);
584  HTag(htag);
585  leg->Draw("same");
586 
587  gPad->SetFillStyle(0);
588  c2->SetFillStyle(0);
589  PrintCanvasAll (c2,out_name+"_datamc","plots/",{"pdf","eps","png","C","root"});
590  PrintLatexFigure(out_name+"_datamc");
591  //removed: prelim tag
592  }
593  else{
594  std::vector<TH1*> vratio;
595  for(unsigned int ii=1;ii<vnom.size();ii++){
596  auto htemp = PrettyRatio (vnom[datahist], vnom[ii], vnom[ii]->GetLineColor(), 1.,
597  "Data / MC");
598  htemp->GetYaxis()->SetRangeUser(0.5,1.5);
599  vratio.push_back(htemp);
600  }
601  std::vector<TString> topOption = {" ","hist"};
602  std::vector<TString> bottomOption (vratio.size(),ratioerrors ?"p":"hist p");
603  auto cratio2= RatioPlot(vnom, topOption,vratio,bottomOption,pidtag,pidaxis);
604  cratio2->cd();
605  leg->Draw("same");
606  HTag(pidtag);
607  PrintCanvasAll(cratio2,out_name+"_datamc_ratio");
608  PrintLatexFigure("/png_"+out_name+"_datamc_ratio");
609  }
610 }
611 */
612 //////////////////////////////////////////////////////////////////////
613 
615  TH1* hdata, std::vector<TH1*> htots, TLegend * leg,
616  TString pidtag="", TString out_name="plot_nd",
617  bool ratioplot=false, AxisType pidaxis = kNue1bin)
618 {
619  Info("PlotNDDataTotalMCComparison",
620  TString::Format("pidaxis = %d", pidaxis));
621 
622  if(pidaxis==kNue3bin)
623  FixLegend(leg, "ExPIDNearDet");
624  std::vector<TString> h_opt;
625  for(auto const & h:htots) h_opt.push_back("hist");
626  h_opt.push_back("p");
627 
628  if(!ratioplot)
629  {
630  auto c1 = new TCanvas ("c1","c1");
631  htots[0]->SetMaximum(1.3*htots[0]->GetMaximum());
632  for(auto const & h:htots) h->Draw("hist same");
633  htots[0]->Draw("hist same"); //need to redraw for error bands
634  hdata->Draw("p same");
635  PIDTag(pidtag);
636  if(pidaxis==kNue3bin)
637  NuePID3bin(htots[0], true); // hist, coreOnlyND, binLabels
638  leg->Draw("same");
639  PrintCanvasAll(c1, out_name+"_datamc");
640  PrintLatexFigure("/png_"+out_name+"_datamc");
641  }
642  else
643  {
644  std::vector<TH1*> ratios;
645  std::vector<TString> r_opt;
646  auto idx = 0;
647  for(auto const & h:htots)
648  {
649  auto ratio = PrettyRatio(hdata, h, h->GetLineColor(),
650  idx>0 ? 0.4:1., "Data / MC");
651  ratio->GetYaxis()->SetRangeUser(0.5,1.5);
652  ratios.push_back(ratio);
653  r_opt.push_back("e");
654  idx++;
655  }
656  htots.push_back(hdata);
657  auto cratio1 = RatioPlot(htots, h_opt, ratios,
658  r_opt, pidtag, pidaxis, true);
659  cratio1->cd();
660  leg->Draw("same");
661  PIDTag(pidtag);
662 
663  PrintCanvasAll(cratio1,out_name+"_datamc_ratio");
664  PrintLatexFigure("/png_"+out_name+"_datamc_ratio");
665  }
666 }
667 //////////////////////////////////////////////////////////////////////
669  std::vector<TH1*>vnom, std::vector<TH1*>vshi,
670  TLegend * leg, TString pidtag="",
671  TString out_name="plot_nd"/* Why ND? */, bool ratioplot=false,
672  bool ratioerrors=false, AxisType pidaxis=kNue1bin)
673 {
674  // Is it NearDet?
675  Info("PlotMCComponentsErrorBand", "out_name = " + out_name);
676  const auto isNearDet = out_name.Contains("_ND_", TString::kIgnoreCase) && !out_name.Contains("_pred_", TString::kIgnoreCase);
677  if (isNearDet)
678  Info("PlotMCComponentsErrorBand", "isNearDet = true");
679  else
680  Info("PlotMCComponentsErrorBand", "isNearDet = false");
681 
682  if(pidaxis==kNue3bin)
683  {
684  if (isNearDet)
685  FixLegend(leg, "ExPIDNearDet");
686  else
687  FixLegend(leg, "ExPID");
688  }
689 
690 
691  if(!ratioplot)
692  {
693  auto c2 = new TCanvas(UniqueName().c_str(), "c2");
694  for(auto const & v:vnom)
695  v->Draw("hist same");
696  for(auto const & v:vshi)
697  v->Draw("hist same");
698  if(pidaxis==kNue3bin)
699  NuePID3bin(vnom[0], isNearDet);
700  PIDTag(pidtag);
701  leg->Draw("same");
702  PrintCanvasAll(c2, out_name+"_mccomp");
703  PrintLatexFigure("/png_"+out_name+"_mccomp");
704  c2->Close();
705  }
706  else
707  {
708  std::vector<TH1*> vratio;
709  for(unsigned int ii=0; ii<vshi.size();ii++)
710  {
711  auto htemp = PrettyRatio(vshi[ii], vnom[ii],
712  vnom[ii]->GetLineColor(), ii>0 ? 0.4:1.);
713  htemp->GetYaxis()->SetRangeUser(0.5, 1.5);
714  vratio.push_back (htemp);
715  }
716  vnom.insert(vnom.end(), vshi.begin(), vshi.end());
717  std::vector<TString>topOption(vnom.size(), "hist");
718  std::vector<TString>bottomOption(vratio.size(),
719  ratioerrors ? "p":"hist p");
720  auto cratio2 = RatioPlot(
721  vnom, topOption, vratio, bottomOption, pidtag, pidaxis, isNearDet);
722  cratio2->cd();
723  leg->Draw("same");
724  PIDTag(pidtag);
725  PrintCanvasAll(cratio2, out_name+"_mccomp_ratio");
726  PrintLatexFigure("/png_"+out_name+"_mccomp_ratio");
727  }
728 }
729 
730 //////////////////////////////////////////////////////////////////////
731 
733  std::vector<TH1*>vnom, std::vector<TH1*>vshi,
734  TLegend * leg, TString pidtag="", TString out_name="plot_nd",
735  bool ratioplot=false, bool ratioerrors=false, AxisType pidaxis=kNue1bin)
736 {
737  // Is it NearDet?
738  Info("PlotMCComponentsErrorBand", "out_name = " + out_name);
739  const auto isNearDet = out_name.Contains("_ND_", TString::kIgnoreCase) && !out_name.Contains("_pred_", TString::kIgnoreCase);
740  if (isNearDet)
741  Info("PlotMCComponentsErrorBand", "isNearDet = true");
742  else
743  Info("PlotMCComponentsErrorBand", "isNearDet = false");
744 
745  auto hplu_band = reinterpret_cast<TH1*>(vshi[0]->Clone());
746  auto hmin_band = reinterpret_cast<TH1*>(vshi[1]->Clone());
747  auto const hplu_hist = reinterpret_cast<TH1*>(vshi[0]->Clone());
748  auto const hmin_hist = reinterpret_cast<TH1*>(vshi[1]->Clone());
749 
750  TString out = "bkg";
751  if(out_name.Contains("_sig")) out = "sig";
752  else if(out_name.Contains("decomp_beam")) out = "beam";
753  else if(out_name.Contains("decomp_nc")) out = "nc";
754  else if(out_name.Contains("decomp_numu")) out = "numu";
755 
756  FormatErrorBand(vnom[0], hplu_band, hmin_band, out, true);
757 
758  std::vector<TH1*>vshi_band={hplu_band,hmin_band};
759  std::vector<TH1*>vshi_hist={hplu_hist,hmin_hist};
760 
761  if(pidaxis==kNue3bin)
762  {
763  if (isNearDet)
764  FixLegend(leg, "ExPIDNearDet");
765  else
766  FixLegend(leg, "ExPID");
767  }
768 
769  if(!ratioplot)
770  {
771  auto c2 = new TCanvas ("c2","c2");
772  vnom[0]->SetMaximum(1.2*vnom[0]->GetMaximum());
773  vnom[0]->Draw("hist");
774  for(auto v:vshi_band) v->Draw("hist same");
775  for(auto v:vshi_hist) v->Draw("e same");
776  for(auto v:vnom) v->Draw("hist same");
777  if(pidaxis==kNue3bin)
778  NuePID3bin(vnom[0], isNearDet);
779  PIDTag(pidtag);
780  leg->Draw("same");
781  PrintCanvasAll(c2, out_name+"_mccomp");
782  PrintLatexFigure("/png_"+out_name+"_mccomp");
783  c2->Close();
784  }
785  else
786  {
787  std::vector<TH1*> vratio;
788  for(unsigned int ii=0; ii<vshi.size(); ii++)
789  {
790  auto htemp = PrettyRatio(vshi[ii], vnom[0], vshi[ii]->GetMarkerColor());
791  vratio.push_back(htemp);
792  }
793  vratio[0]->GetYaxis()->SetRangeUser(0.5,1.5);
794  vnom.insert(vnom.begin()+1,vshi_band.begin(),vshi_band.end());
795  std::vector<TString> topOption (vnom.size(), "hist");
796  std::vector<TString> bottomOption (vratio.size(),ratioerrors ?"p":"hist p");
797  auto cratio2= RatioPlot(
798  vnom, topOption, vratio, bottomOption, pidtag, pidaxis, isNearDet);
799  cratio2->cd();
800  leg->Draw("same");
801  PIDTag(pidtag);
802  PrintCanvasAll(cratio2, out_name+"_mccomp_ratio");
803  PrintLatexFigure("/png_"+out_name+"_mccomp_ratio");
804  }
805 }
806 
807 //////////////////////////////////////////////////////////////////////
808 
809 
810 
811 //////////////////////////////////////////////////////////////////////
812 // Define container for histograms + format
813 //////////////////////////////////////////////////////////////////////
814 
815 namespace ana
816 {
817 
818  struct DataMCComponents
819  {
820  TH1* data;
821  TH1* mcTotal;
822  TH1* sig, *wrongSign, *bkg, *beam, *numu, *nc, *tau;
824  };
825 
826  void DefaultFormatNue (DataMCComponents & comp, int linestyle = 1){
827  std::vector <std::pair <TH1*, int>> pairs = {
828  {comp.sig, kNueSignalColor},
829  {comp.wrongSign, kNueSignalWSColor},
830  {comp.beam, kBeamNueBackgroundColor},
831  {comp.nc, kNCBackgroundColor},
832  {comp.numu, kNumuBackgroundColor},
833  {comp.bkg, kTotalMCColor},
834  {comp.mcTotal, kTotalMCColor},
835  {comp.nuetonumu, kNumuBackgroundColor},
836  {comp.numutonumu, kNumuBackgroundColor},
837  {comp.nuetonutau, kGray},
838  {comp.numutonutau, kGray}
839  };
840  for (auto const & mypair:pairs)
841  {
842  if(mypair.first == 0)
843  {
844  Info("DefaultFormatNue", "Skipping formatting for component...");
845  continue;
846  }
847  mypair.first->SetLineColor(mypair.second);
848  mypair.first->SetLineStyle(linestyle);
849  ana::CenterTitles(mypair.first);
850  }
851  if(comp.data) comp.data->SetMarkerStyle(kFullCircle);
852  return;
853  }
854 //////////////////////////////////////////////////////////////////////
855 
856  const DataMCComponents GetNDCompsFromDecomp(const IDecomp* decomp)
857  {
858  // Get correct POT for horn current
859  double ndPOT;
860 
861  if (isFHC)
862  ndPOT = kAna2020FHCPOT;
863  else
864  ndPOT = kAna2020RHCPOT;
865 
866  // Get components from decomp
867  Spectrum nueComp = decomp->MC_NueComponent();
868  Info("GetNDComponents", "Loaded NueComponent...");
869  Spectrum antinueComp = decomp->MC_AntiNueComponent();
870  Info("GetNDComponents", "Loaded AnitNueComponent...");
871  Spectrum numuComp = decomp->MC_NumuComponent();
872  Info("GetNDComponents", "Loaded NumuComponent...");
873  Spectrum antinumuComp = decomp->MC_AntiNumuComponent();
874  Info("GetNDComponents", "Loaded AntiNumuComponent...");
875  Spectrum ncTotComp = decomp->MC_NCTotalComponent();
876  Info("GetNDComponents", "Loaded NCTotalComponent...");
877  Spectrum dataComp = decomp->Data_Component();
878  Info("GetNDComponents", "Loaded Data_Component...");
879 
880  auto beam = nueComp + antinueComp;
881  auto numu = numuComp + antinumuComp;
882 
883  DataMCComponents ret;
884  ret.beam = beam.ToTH1(ndPOT);
885  std::cout << "beam component: " << ret.beam->Integral()
886  << ", in " << ret.beam->GetNbinsX() << " bins" << std::endl;
887  ret.mcTotal = reinterpret_cast<TH1*>(ret.beam->Clone());
888 
889  ret.numu = numu.ToTH1(ndPOT);
890  std::cout << "numu component: " << ret.numu->Integral()
891  << ", in " << ret.numu->GetNbinsX() << " bins" << std::endl;
892  ret.mcTotal->Add(ret.numu);
893 
894  ret.nc = ncTotComp.ToTH1(ndPOT);
895  std::cout << "nc component: " << ret.nc->Integral()
896  << ", in " << ret.nc->GetNbinsX() << " bins" << std::endl;
897  ret.mcTotal->Add(ret.nc);
898 
899  std::cout << "Total MC: " << ret.mcTotal->Integral()
900  << ", in " << ret.mcTotal->GetNbinsX() << " bins" << std::endl;
901 
902  ret.data = dataComp.ToTH1(ndPOT);
903  std::cout << "data component: " << ret.data->Integral()
904  << ", in " << ret.beam->GetNbinsX() << " bins" << std::endl;
905 
906  ret.bkg = reinterpret_cast<TH1*>(ret.mcTotal->Clone());
907  ret.sig = 0;
908  ret.wrongSign = 0;
909  ret.tau = 0;
910  ret.nuetonumu = 0;
911  ret.nuetonutau = 0;
912  ret.numutonumu = 0;
913  ret.numutonutau = 0;
914 
915  std::cout << "Loaded ND components" << std::endl;
916 
917  // Add formatting
918  // TODO: it's quite unintuitive that all formatting goes through
919  // DefaultFormatNue, even for Numu!
920  DefaultFormatNue(ret);
921 
922  return ret;
923  }
924 
925 //////////////////////////////////////////////////////////////////////
926 
928  {
929  // This is to get the MC components after the decomposition (Total MC = Data)
930  // Get correct POT for horn current
931  double ndPOT;
932 
933  if (isFHC)
934  ndPOT = kAna2020FHCPOT;
935  else
936  ndPOT = kAna2020RHCPOT;
937 
938  // Get components from decomp
939  Spectrum nueComp = decomp->NueComponent();
940  Info("GetNDComponents", "Loaded NueComponent...");
941  Spectrum antinueComp = decomp->AntiNueComponent();
942  Info("GetNDComponents", "Loaded AnitNueComponent...");
943  Spectrum numuComp = decomp->NumuComponent();
944  Info("GetNDComponents", "Loaded NumuComponent...");
945  Spectrum antinumuComp = decomp->AntiNumuComponent();
946  Info("GetNDComponents", "Loaded AntiNumuComponent...");
947  Spectrum ncTotComp = decomp->NCTotalComponent();
948  Info("GetNDComponents", "Loaded NCTotalComponent...");
949  Spectrum dataComp = decomp->Data_Component();
950  Info("GetNDComponents", "Loaded Data_Component...");
951 
952  auto beam = nueComp + antinueComp;
953  auto numu = numuComp + antinumuComp;
954 
956  ret.beam = beam.ToTH1(ndPOT);
957  std::cout << "beam component: " << ret.beam->Integral()
958  << ", in " << ret.beam->GetNbinsX() << " bins" << std::endl;
959  ret.mcTotal = reinterpret_cast<TH1*>(ret.beam->Clone());
960 
961  ret.numu = numu.ToTH1(ndPOT);
962  std::cout << "numu component: " << ret.numu->Integral()
963  << ", in " << ret.numu->GetNbinsX() << " bins" << std::endl;
964  ret.mcTotal->Add(ret.numu);
965 
966  ret.nc = ncTotComp.ToTH1(ndPOT);
967  std::cout << "nc component: " << ret.nc->Integral()
968  << ", in " << ret.nc->GetNbinsX() << " bins" << std::endl;
969  ret.mcTotal->Add(ret.nc);
970 
971  std::cout << "Total MC: " << ret.mcTotal->Integral()
972  << ", in " << ret.mcTotal->GetNbinsX() << " bins" << std::endl;
973 
974  ret.data = dataComp.ToTH1(ndPOT);
975  std::cout << "data component: " << ret.data->Integral()
976  << ", in " << ret.beam->GetNbinsX() << " bins" << std::endl;
977 
978  ret.bkg = reinterpret_cast<TH1*>(ret.mcTotal->Clone());
979  ret.sig = 0;
980  ret.wrongSign = 0;
981  ret.tau = 0;
982  ret.nuetonumu = 0;
983  ret.nuetonutau = 0;
984  ret.numutonumu = 0;
985  ret.numutonutau = 0;
986 
987  std::cout << "Loaded ND components" << std::endl;
988 
989  // Add formatting
990  // TODO: it's quite unintuitive that all formatting goes through
991  // DefaultFormatNue, even for Numu!
992  DefaultFormatNue(ret);
993 
994  return ret;
995  }
996 
997 //////////////////////////////////////////////////////////////////////
998 
999  struct DataMCComponents GetFDMCComponents(
1000  osc::IOscCalc * calc, IPrediction * pred_no, TString
1001  output_name="nue", int linestyle=1, bool bkgdetails=false)
1002  {
1004  // Set some defaults for the components.
1005  Flavors::Flavors_t MySig, MyBeam, MyNuMu, MyNC, MyTau;
1006  MySig = MyBeam = MyNuMu = MyNC = MyTau = Flavors::kAll;
1007  // Change the components depending on what I'm looking at.
1008  if (!isNumuAna)
1009  { // Nu e
1010  MySig = Flavors::kNuMuToNuE;
1011  MyBeam = Flavors::kNuEToNuE;
1012  MyNuMu = Flavors::kAllNuMu;
1013  MyNC = Flavors::kAll;
1014  MyTau = Flavors::kAllNuTau;
1015  }
1016  else
1017  { // Nu mu
1018  MySig = Flavors::kNuMuToNuMu;
1019  MyBeam = Flavors::kAllNuE;
1020  MyNuMu = Flavors::kNuEToNuMu;
1021  MyNC = Flavors::kAll;
1022  MyTau = Flavors::kAllNuTau;
1023  }
1024 
1025  // Specify sign depending on polarity
1027  sign = wrongSign = Sign::kBoth;
1028  if (!isNumuAna)
1029  {
1030  if (isFHC)
1031  {
1032  sign = Sign::kNu;
1033  wrongSign = Sign::kAntiNu;
1034  }
1035  else
1036  {
1037  sign = Sign::kAntiNu;
1038  wrongSign = Sign::kNu;
1039  }
1040  }
1041  const double pot = (isFHC) ? kAna2020FHCPOT:kAna2020RHCPOT;
1042 
1043  // Predict components
1044  ret.sig = pred_no->PredictComponent(
1045  calc, MySig, Current::kCC, sign).ToTH1(pot);
1046  ret.wrongSign = pred_no->PredictComponent(
1047  calc, MySig, Current::kCC, wrongSign).ToTH1(pot);
1048  ret.beam = pred_no->PredictComponent(
1049  calc, MyBeam, Current::kCC, Sign::kBoth).ToTH1(pot);
1050  ret.numu = pred_no->PredictComponent(
1051  calc, MyNuMu, Current::kCC, Sign::kBoth).ToTH1(pot);
1052  ret.nc = pred_no->PredictComponent(
1053  calc, MyNC, Current::kNC, Sign::kBoth).ToTH1(pot);
1054  ret.tau = pred_no->PredictComponent(
1055  calc, MyTau, Current::kCC, Sign::kBoth).ToTH1(pot);
1056 
1057  ret.bkg = (pred_no->Predict(calc) - pred_no->PredictComponent(
1058  calc, MySig, Current::kCC, sign)).ToTH1(pot);
1059 
1060  ret.mcTotal = reinterpret_cast<TH1*>(ret.bkg->Clone());
1061  ret.mcTotal->Add(ret.sig);
1062 
1063  if(bkgdetails)
1064  {
1065  ret.nuetonumu = pred_no->PredictComponent(
1067  ret.numutonumu = pred_no->PredictComponent(
1069  ret.nuetonutau = pred_no->PredictComponent(
1071  ret.numutonutau = pred_no->PredictComponent(
1073  }
1074  else
1075  {
1076  ret.nuetonumu = 0; ret.numutonumu = 0;
1077  ret.nuetonutau = 0; ret.numutonutau = 0;
1078  }
1079 
1080  ret.data = 0;
1081  DefaultFormatNue(ret, linestyle);
1082 
1083  return ret;
1084  }
1085 
1086 //////////////////////////////////////////////////////////////////////
1087 // Compare ND Data MC
1088 //////////////////////////////////////////////////////////////////////
1089 
1091  TString plottitle, TString out_name, TString pidtag,
1092  AxisType pidaxis, bool printtable=true)
1093  {
1094  h_sh.bkg->SetLineColor(kTotalMCErrorBandColor);
1095  h_sh.bkg->SetLineStyle(1);
1096 
1097  h_no.bkg->SetMaximum(FindLimitY({h_no.mcTotal, h_no.data}));
1098  h_no.bkg->SetTitle(plottitle);
1099 
1100  TGaxis::SetMaxDigits(3);
1101  h_no.bkg->GetXaxis()->SetDecimals();
1102 
1103  auto leg1 = CustomLegend(
1104  {h_no.data, h_no.mcTotal, h_sh.mcTotal},
1105  {"ND Data","Nominal","Shift"},{"epl","l","l"});
1106 
1107  auto leg2 = DefaultNueLegend(false, true, false);
1108 
1109  TH1* mydata = h_no.data;
1110  std::vector <TH1*> comp_plot = {h_no.mcTotal, h_sh.mcTotal};
1111  std::vector <TH1*> vnom_plot = {h_no.mcTotal, h_no.beam, h_no.nc, h_no.numu};
1112  std::vector <TH1*> vshi_plot = {h_sh.mcTotal};
1113 
1114  std::vector <TH1*> vnom = {h_no.mcTotal, h_no.beam, h_no.nc, h_no.numu};
1115  std::vector <TH1*> vshi = {h_sh.mcTotal, h_sh.beam, h_sh.nc, h_sh.numu};
1116  std::vector<TString> labels = {"Total MC", "Beam $\\nu_{e}$", "NC", "$\\nu_{\\mu}$ CC"};
1117 
1118  // Formatting for shift --> do this before re-binning for Numu
1119  h_sh.mcTotal->SetLineStyle(2);
1120  h_sh.mcTotal->SetLineColor(kTotalMCColor+1);
1121 
1122  if(pidaxis == kNumuOpt)
1123  {
1124  leg2 = DefaultNumuLegend(false, true, false);
1125  mydata = HistNumuOpt(h_no.data);
1126  comp_plot = {HistNumuOpt(h_no.mcTotal), HistNumuOpt(h_sh.mcTotal)};
1127  vnom_plot = {HistNumuOpt(h_no.mcTotal), HistNumuOpt(h_no.numu),
1128  HistNumuOpt(h_no.nc), HistNumuOpt(h_no.beam)};
1129  vshi_plot = {HistNumuOpt(h_sh.mcTotal)};
1130 
1131  vnom = {HistNumuOpt(h_no.mcTotal), HistNumuOpt(h_no.numu),
1132  HistNumuOpt(h_no.nc), HistNumuOpt(h_no.beam)};
1133  vshi = {HistNumuOpt(h_sh.mcTotal), HistNumuOpt(h_sh.numu),
1134  HistNumuOpt(h_sh.nc), HistNumuOpt(h_sh.beam)};
1135  labels = {"Total MC", "$\\nu_{\\mu}$ CC", "NC", "$\\nu_{e}$ CC"};
1136  }
1137 
1139  mydata, comp_plot, leg1, pidtag, out_name, printratioplots, pidaxis);
1140 
1141  leg2->AddEntry(h_sh.mcTotal, "Shifted", "l");
1142 
1144  vnom_plot, vshi_plot, leg2, pidtag,
1145  out_name, printratioplots, false, pidaxis);
1146 
1147  if(printtable)
1148  {
1149  ComparisonTableOne(vnom, vshi, labels, out_name);
1150 
1151  int nbins = 1;
1152  if (pidaxis== kNue3bin)
1153  nbins = 2; // No peripheral bin for ND
1154  ComparisonTableOneNbins(vnom, vshi, labels, out_name, pidaxis);
1155  }
1156  }
1157 //////////////////////////////////////////////////////////////////////
1158 
1161  TString plottitle, TString out_name, TString pidtag,
1162  AxisType pidaxis, bool printtable=true, bool decomps = false)
1163  {
1164  Info("CompareNDDataTwoMC",
1165  TString::Format("pidaxis = %d", pidaxis));
1166  // std::cout << "DEBUG[CompareNDDataTwoMC]: hnom.bkg "
1167  // << hnom.bkg->Integral() << std::endl;
1168  // std::cout << "DEBUG[CompareNDDataTwoMC]: hplu.bkg "
1169  // << hplu.bkg->Integral() << std::endl;
1170  // std::cout << "DEBUG[CompareNDDataTwoMC]: hmin.bkg "
1171  // << hmin.bkg->Integral() << std::endl;
1172  //
1173  // std::cout << "DEBUG[CompareNDDataTwoMC]: hnom.data "
1174  // << hnom.data->Integral() << std::endl;
1175  // std::cout << "DEBUG[CompareNDDataTwoMC]: hplu.data "
1176  // << hplu.data->Integral() << std::endl;
1177  // std::cout << "DEBUG[CompareNDDataTwoMC]: hmin.data "
1178  // << hmin.data->Integral() << std::endl;
1179 
1180  // Call here to get correct formatting in legends
1181  // Need to call again with fixbins as true later
1182  FormatErrorBand(hnom.mcTotal, hplu.mcTotal, hmin.mcTotal, "bkg");
1183  if(decomps)
1184  {
1185  FormatErrorBand(hnom.beam, hplu.beam, hmin.beam, "beam");
1186  FormatErrorBand(hnom.numu, hplu.numu, hmin.numu, "numu");
1187  FormatErrorBand(hnom.nc, hplu.nc, hmin.nc, "nc");
1188  }
1189 
1190  hnom.mcTotal->SetMaximum(FindLimitY(
1191  {hnom.mcTotal, hnom.data, hplu.mcTotal, hmin.mcTotal}));
1192  hnom.mcTotal->SetTitle(plottitle);
1193 
1194  TGaxis::SetMaxDigits(3);
1195  hnom.mcTotal->GetXaxis()->SetDecimals();
1196 
1197  auto leg1 = CustomLegend(
1198  {hnom.data, hnom.mcTotal, hplu.mcTotal},
1199  {"ND Data","Nominal MC","Syst."}, {"epl","l","f"});
1200  auto leg2 = DefaultNueLegend(false, true, false);
1201 
1202  TH1 * mydata = hnom.data;
1203  std::vector <TH1*> comp_plot = {hnom.mcTotal, hplu.mcTotal, hmin.mcTotal};
1204  std::vector <TH1*> vnom_plot = {hnom.mcTotal, hnom.beam, hnom.nc, hnom.numu};
1205  std::vector <TH1*> vshi_plot = {hplu.mcTotal, hmin.mcTotal};
1206 
1207  std::vector <TH1*> vnom_plot_beam = {hnom.beam, hnom.beam, hnom.beam, hnom.beam};
1208  std::vector <TH1*> vnom_plot_nc = {hnom.nc, hnom.nc, hnom.nc, hnom.nc};
1209  std::vector <TH1*> vnom_plot_numu = {hnom.numu, hnom.numu, hnom.numu, hnom.numu};
1210 
1211  std::vector <TH1*> vshi_plot_beam = {hplu.beam, hmin.beam};
1212  std::vector <TH1*> vshi_plot_numu = {hplu.numu, hmin.numu};
1213  std::vector <TH1*> vshi_plot_nc = {hplu.nc, hmin.nc};
1214 
1215  std::vector <TH1*> vnom = {hnom.mcTotal, hnom.beam, hnom.nc, hnom.numu};
1216  std::vector <TH1*> vplu = {hplu.mcTotal, hplu.beam, hplu.nc, hplu.numu};
1217  std::vector <TH1*> vmin = {hmin.mcTotal, hmin.beam, hmin.nc, hmin.numu};
1218  std::vector<TString> labels = {
1219  "Total MC", "Beam $\\nu_{e}$", "NC", "$\\nu_{\\mu}$ CC"};
1220 
1221  if(pidaxis == kNumuOpt)
1222  {
1223  leg2 = DefaultNumuLegend(false, true, false);
1224  mydata = HistNumuOpt(hnom.data);
1225  comp_plot = {HistNumuOpt(hnom.mcTotal),
1226  HistNumuOpt(hplu.mcTotal), HistNumuOpt(hmin.mcTotal)};
1227  vnom_plot = {
1228  HistNumuOpt(hnom.mcTotal),
1229  HistNumuOpt(hnom.numu),
1230  HistNumuOpt(hnom.nc),
1231  HistNumuOpt(hnom.beam)
1232  };
1233  vshi_plot = {HistNumuOpt(hplu.mcTotal), HistNumuOpt(hmin.mcTotal)};
1234 
1235  // vnom = {hnom.mcTotal, hnom.numu, hnom.nc, hnom.beam};
1236  // vplu = {hplu.mcTotal, hplu.numu, hplu.nc, hplu.beam};
1237  // vmin = {hmin.mcTotal, hmin.numu, hmin.nc, hmin.beam};
1238  vnom = {HistNumuOpt(hnom.mcTotal), HistNumuOpt(hnom.numu),
1239  HistNumuOpt(hnom.nc), HistNumuOpt(hnom.beam)};
1240  vplu = {HistNumuOpt(hplu.mcTotal), HistNumuOpt(hplu.numu),
1241  HistNumuOpt(hplu.nc), HistNumuOpt(hplu.beam)};
1242  vmin = {HistNumuOpt(hmin.mcTotal), HistNumuOpt(hmin.numu),
1243  HistNumuOpt(hmin.nc), HistNumuOpt(hmin.beam)};
1244  labels = {"Total MC", "$\\nu_{\\mu}$ CC", "NC", "$\\nu_{e}$ CC"};
1245  }
1246 
1247  // Format error band on data/MC plot.
1248  // Do this here because function is shared with single shift case.
1249  FormatErrorBand(comp_plot[0], comp_plot[1], comp_plot[2],
1250  out_name.Contains("_sig"), "sig");
1251 
1253  mydata, comp_plot, leg1, pidtag,
1254  out_name, printratioplots, pidaxis);
1255 
1256  leg2->AddEntry(hplu.mcTotal, "Syst.", "f");
1257 
1259  vnom_plot, vshi_plot, leg2, pidtag,
1260  out_name, printratioplots, false, pidaxis);
1261 
1262  if(decomps)
1263  {
1265  vnom_plot_beam, vshi_plot_beam, leg2, pidtag,
1266  out_name+"_beam", printratioplots, false, pidaxis);
1268  vnom_plot_nc, vshi_plot_nc, leg2, pidtag,
1269  out_name+"_nc", printratioplots, false, pidaxis);
1271  vnom_plot_numu, vshi_plot_numu, leg2, pidtag,
1272  out_name+"_numu", printratioplots, false, pidaxis);
1273  }
1274 
1275  // TODO: seg fault when ratioplot is false
1276  // --> do we ever not want to plot the ratio?
1277 
1278  if(printtable)
1279  {
1280  ComparisonTable(vnom, vplu, vmin, labels, out_name);
1281  int nbins = 1;
1282  if (pidaxis== kNue3bin)
1283  nbins = 2; // no peripheral bin for ND
1284  ComparisonTableNbins(vnom, vplu, vmin, labels, out_name, kNueNDbin);
1285  }
1286  }
1287 
1288 //////////////////////////////////////////////////////////////////////
1289 
1291  {
1296  // NC
1298  // End NC
1301  };
1302 
1303  const IDecomp* GetDecomp(IPrediction* prediction,
1304  EModExtrapComps modExtrapComp)
1305  {
1306  // First cast IPrediction as a PredictionExtrap object and check it extist
1307  std::cout << "Checking if IPrediction exists" << std::endl;
1308  if (!prediction)
1309  abort();
1310 
1311  PredictionExtrap* predExtrap;
1312 
1313  // If numu then it's simmple;
1314  if (isNumuAna)
1315  {
1316  std::cout << abi::__cxa_demangle(typeid(*prediction).name(), 0, 0, 0) << std::endl;
1317  predExtrap = dynamic_cast<PredictionExtrap*>(prediction);
1318  }
1319  else
1320  {
1321  // If nue then we have to do some wrangling to get the correct object
1322  std::cout << abi::__cxa_demangle(typeid(*prediction).name(), 0, 0, 0) << std::endl;
1323  // This will be a PredictionExtendToPeripheral so cast as that
1324  PredictionExtendToPeripheral* predExtend = dynamic_cast<PredictionExtendToPeripheral*>(prediction);
1325  std::cout << abi::__cxa_demangle(typeid(*predExtend).name(), 0, 0, 0) << std::endl;
1326  std::cout << "Checking if PredictionExtendToPeripheral exists" << std::endl;
1327  if (!predExtend)
1328  abort();
1329 
1330  // Now we want to get the core part of this
1331  IPrediction* predCore = predExtend->GetCore();
1332  std::cout << "Checking if IPrediction (core) exists" << std::endl;
1333  if (!predCore)
1334  abort();
1335 
1336  // Now we want to cast this as a PredictionExtrap
1337  std::cout << abi::__cxa_demangle(typeid(*predCore).name(), 0, 0, 0) << std::endl;
1338  predExtrap = dynamic_cast<PredictionExtrap*>(predCore);
1339  std::cout << "Checking if PredictionExtrap exists" << std::endl;
1340  std::cout << abi::__cxa_demangle(typeid(*predExtrap).name(), 0, 0, 0) << std::endl;
1341  }
1342  if (!predExtrap)
1343  abort();
1344  // Get the Extrap object
1345  IExtrap* extrap = const_cast<IExtrap*>(predExtrap->GetExtrap());
1346 
1347  // Cast as modular extrap
1348  ModularExtrap* modExtrap = dynamic_cast<ModularExtrap*>(extrap);
1349 
1350  // Get mod extrap components
1351  std::vector<ModularExtrapComponent*> modExtrapComps =
1352  modExtrap->GetModExtrapComponents();
1353 
1354  // Get decomp for EEextrap
1355  // Use vector index to select appropriate set of extrap components
1356  std::cout << "*************************************" << std::endl;
1357  std::cout << "Getting modular extrap component: "
1358  << modExtrapComp << std::endl;
1359  std::cout << "*************************************" << std::endl;
1360  const IDecomp* decomp = modExtrapComps[modExtrapComp]->GetDecomp();
1361 
1362  return decomp;
1363  }
1364 
1365 //////////////////////////////////////////////////////////////////////
1366 
1368  PredictionSyst3Flavor2020* predictionSyst, const ISyst* syst,
1369  EModExtrapComps modExtrapComp, TString plottitle, TString out_name,
1370  TString tag, AxisType pidaxis, bool printtable=true)
1371  {
1372  Info("CompareNDDataMCFromVector",
1373  TString::Format("pidaxis = %d", pidaxis));
1374 
1375  // Get nominal from PredictionSyst3Flavor2020
1376  std::cout << std::endl;
1377  std::cout << "Getting Nominal prediction" << std::endl;
1378  auto pred_no = predictionSyst->GetNominalPrediction();
1379  std::cout << "In ComparePredictionsFromVector -- pred_no is " << pred_no << std::endl;
1380  // Get Nominal decomp
1381  // Option 1: through decomp
1382  const IDecomp* nomDecomp = GetDecomp(pred_no, modExtrapComp);
1383  DataMCComponents nomNDComps = GetNDCompsFromDecomp(nomDecomp);
1384  DataMCComponents nomNDDecomps = GetNDDecompsFromDecomp(nomDecomp);
1385 
1386  std::cout << "Got nomDecomp" << std::endl;
1387  std::cout << "Got nomNDComps" << std::endl;
1388  std::cout << "Got nomNDDecomps" << std::endl;
1389  // Get shifted predictions from PredictionSystJoint2020
1390  const auto& shiftedPred = predictionSyst->GetShiftedPrediction(syst);
1391 
1392  std::cout << std::endl;
1393  std::cout << "ShiftedPreds name: " << shiftedPred.systName << std::endl;
1394 
1395  int numSh = shiftedPred.preds.size();
1396  std::cout << "Number of shifts: " << numSh << std::endl;
1397  // for (auto shift:shiftedPred.preds) { std::cout << "Func shift: " << shift << std::endl; }
1398 
1399  std::vector<IPrediction*> pred_sh = shiftedPred.preds;
1400 
1401  // create vecor of components
1402  std::vector<DataMCComponents> comps;
1403  std::vector<DataMCComponents> decomps;
1404  for (auto pred:pred_sh)
1405  {
1406  const IDecomp* shiftedDecomp = GetDecomp(pred, modExtrapComp);
1407  std::cout << "Got shiftedDecomp" << std::endl;
1408  DataMCComponents shiftedComps = GetNDCompsFromDecomp(shiftedDecomp);
1409  std::cout << "Got shiftedComps" << std::endl;
1410  DataMCComponents shiftedDecomps = GetNDDecompsFromDecomp(shiftedDecomp);
1411  std::cout << "Got shiftedDecomps" << std::endl;
1412  std::cout << std::endl;
1413  comps.push_back(shiftedComps);
1414  decomps.push_back(shiftedDecomps);
1415  }
1416 
1417  std::cout << "comps has length: " << comps.size() << std::endl;
1418 
1419  if (numSh == 3)
1420  {
1421  // shift 0 = -1
1422  // shift 1 = Nominal
1423  // shift 2 = +1
1424  //std::cout << "Entering shift = 3" << std::endl;
1426  nomNDComps, comps[2], comps[0], plottitle,
1427  out_name+"_1sigma_undecomp", "#pm1#sigma", pidaxis, printtable);
1428 
1429  /*CompareNDDataTwoMC(
1430  nomNDDecomps, decomps[2], decomps[0], plottitle,
1431  out_name+"_1sigm_decomp", "#pm1#sigma", pidaxis, printtable, true);*/
1432  }
1433  else if (numSh == 5)
1434  {
1435  // shift 1 = -1
1436  // shift 2 = Nominal
1437  // shift 3 = +1
1438  //std::cout << "Entering shift = 5" << std::endl;
1440  nomNDComps, comps[3], comps[1], plottitle,
1441  out_name+"_1sigma", "#pm1#sigma", pidaxis, printtable);
1442 
1443  /*CompareNDDataTwoMC(
1444  nomNDDecomps, decomps[3], decomps[1], plottitle,
1445  out_name+"_1sigm_decomp", "#pm1#sigma", pidaxis, printtable, true);*/
1446  }
1447  else if (numSh == 4)
1448  {
1449  // fake negative shifts for one-sided systs
1450  // shift 0,1 = 0
1451  // shift 2 = Nominal
1452  // shift 3 = +1
1454  nomNDComps, comps[3], plottitle, out_name+"_1sigma",
1455  "#plus1#sigma", pidaxis, printtable);
1456  }
1457  else if (numSh == 2)
1458  {
1459  // shift 0 = Nominal
1460  // shift 1 = +1
1461  //std::cout << "Entering shift = 2" << std::endl;
1463  nomNDComps, comps[1], plottitle, out_name+"_1sigma",
1464  "#plus1#sigma", pidaxis, printtable);
1465  }
1466  }
1467 
1468 //////////////////////////////////////////////////////////////////////
1469 // Compare FD Predictions
1470 //////////////////////////////////////////////////////////////////////
1471 
1472  void CompareOneShiftPred(IPrediction * pred_no, IPrediction * pred_sh,
1473  TString plottitle, TString out_name, TString pidtag,
1474  bool printtable=true, AxisType pidaxis = kNue1bin, bool PrintParams=false){
1475 
1476  auto calc = GetCalculator (usedumb);
1477  auto h_no = GetFDMCComponents (calc, pred_no, out_name);
1478  auto h_sh = GetFDMCComponents (calc, pred_sh, out_name, 2);
1479 
1480  h_no.sig->SetMaximum(FindLimitY({h_no.sig,h_sh.sig}));
1481  h_no.bkg->SetMaximum(FindLimitY({h_no.bkg,h_sh.bkg}));
1482  h_no.sig->SetTitle(plottitle);
1483  h_no.bkg->SetTitle(plottitle);
1484  TGaxis::SetMaxDigits(3);
1485 
1486  auto leg = DefaultNueLegend(false,true,false);
1487  if (isNumuAna) leg = DefaultNumuLegend(false,true,false);
1488  leg->AddEntry(DummyGraph(kGray+1,2),"Shifted","l");
1489 
1490  std::vector<TString> legnames = {"#nu_{e} Signal","Shifted"};
1491  if (isNumuAna) legnames = {"#nu_{#mu} Signal","Shifted"};
1492  auto leg2 = CustomLegend({h_no.sig, h_sh.sig},legnames,{"l","l"});
1493 
1494  if(pidaxis == kNumuOpt)
1495  { //fix the plot but don't mess up the tables
1496 
1497  PlotMCComponentsComparison({ HistNumuOpt(h_no.bkg), HistNumuOpt(h_no.beam), HistNumuOpt(h_no.nc), HistNumuOpt(h_no.numu) },
1498  { HistNumuOpt(h_sh.bkg), HistNumuOpt(h_sh.beam), HistNumuOpt(h_sh.nc), HistNumuOpt(h_sh.numu) },
1499  leg,pidtag,out_name+"_bkg",printratioplots,false,pidaxis);
1501  { HistNumuOpt(h_sh.sig) },
1502  leg2,pidtag,out_name+"_sig",printratioplots,false,pidaxis);
1503  }
1504  else
1505  {
1506  PlotMCComponentsComparison({ h_no.bkg, h_no.beam, h_no.nc, h_no.numu },
1507  { h_sh.bkg, h_sh.beam, h_sh.nc, h_sh.numu },
1508  leg,pidtag, out_name+"_bkg", printratioplots, false, pidaxis);
1509  PlotMCComponentsComparison({ h_no.sig },
1510  { h_sh.sig },
1511  leg2,pidtag,out_name+"_sig",printratioplots,false,pidaxis);
1512  }
1513 
1514  if(PrintParams) PrintOscilationParameters (calc, usedumb);
1515 
1516  if(printtable){
1517  std::vector<TString> labels = {"$\\nu_e$ signal","Tot beam bkg","Beam $\\nu_{e}$ CC","NC", "$\\nu_{\\mu}$ CC", "$\\nu_{\\tau}$ CC"};
1518 
1519  if (isNumuAna) {
1520  labels = {"$\\nu_{\\mu}$ signal","Tot beam bkg","$\\nu_{e}$ CC","NC", "$\\nu_{\\mu} $App", "$\\nu_{\\tau}$ CC"};
1521  }
1522  ComparisonTableOne ({h_no.sig, h_no.bkg, h_no.beam, h_no.nc, h_no.numu, h_no.tau},
1523  {h_sh.sig, h_sh.bkg, h_sh.beam, h_sh.nc, h_sh.numu, h_sh.tau},
1524  labels,
1525  out_name
1526  );
1527  int nbins = 1;
1528  if (pidaxis== kNue3bin) nbins = 3;
1529  ComparisonTableOneNbins({h_no.sig, h_no.bkg, h_no.beam,h_no.nc,h_no.numu,h_no.tau},
1530  {h_sh.sig, h_sh.bkg, h_sh.beam,h_sh.nc,h_sh.numu,h_sh.tau},
1531  labels,
1532  out_name,
1533  pidaxis);
1534  }//end if table
1535  }
1536 //////////////////////////////////////////////////////////////////////
1537 
1538  void CompareTwoShiftPred(
1539  IPrediction * pred_no, IPrediction * pred_pl, IPrediction * pred_mi,
1540  TString plottitle, TString out_name, TString pidtag,
1541  bool printtable=true, AxisType pidaxis = kNue1bin, bool PrintParams=false)
1542  {
1543 
1544  auto calc = GetCalculator(usedumb);
1545 
1546  auto hnom = GetFDMCComponents(calc, pred_no, out_name);
1547  auto hplu = GetFDMCComponents(calc, pred_pl, out_name);
1548  auto hmin = GetFDMCComponents(calc, pred_mi, out_name);
1549 
1550  //FormatErrorBand(hnom.bkg, hplu.bkg, hmin.bkg, false);
1551  std::cout << "DEBUG! " << hnom.bkg->Integral() << " " << hplu.bkg->Integral() << " " << hmin.bkg->Integral() << std::endl;
1552  FormatErrorBand(hnom.bkg, hplu.bkg, hmin.bkg, "bkg");
1553  hnom.bkg->SetMaximum(FindLimitY({hnom.bkg, hplu.bkg, hmin.bkg}));
1554  hnom.bkg->SetTitle(plottitle);
1555  std::cout << "DEBUG! " << hnom.bkg->Integral() << " " << hplu.bkg->Integral() << " " << hmin.bkg->Integral() << std::endl;
1556 
1557  FormatErrorBand(hnom.sig, hplu.sig, hmin.sig, "sig");
1558  hnom.sig->SetMaximum(FindLimitY({hnom.sig, hplu.sig, hmin.sig}));
1559  hnom.sig->SetTitle(plottitle);
1560  std::cout << "DEBUG! " << hnom.bkg->Integral() << " " << hplu.bkg->Integral() << " " << hmin.bkg->Integral() << std::endl;
1561 
1562  TGaxis::SetMaxDigits(3);
1563 
1564  TLegend * leg, *leg2;
1565  if( pidaxis == kNumuOpt || pidaxis == kNumuOther)
1566  {
1567  leg = DefaultNumuLegend(false, true, false);
1568  leg2 = CustomLegend({hnom.sig, hplu.sig},
1569  {"#nu_{#mu} Signal","Shifted"},
1570  {"l","f"});
1571  }
1572  else
1573  {
1574  leg = DefaultNueLegend(false, true, false);
1575  leg2 = CustomLegend({hnom.sig, hplu.sig},
1576  {"#nu_{e} Signal","Shifted"},
1577  {"l","f"});
1578  }
1579  leg->AddEntry(hplu.bkg,"Syst.","f");
1580 
1581  std::vector <TH1*> nom_bkg_plot = {hnom.bkg, hnom.beam, hnom.nc, hnom.numu, hnom.tau};
1582  std::vector <TH1*> shi_bkg_plot = {hplu.bkg, hmin.bkg};
1583  std::vector <TH1*> nom_sig_plot = {hnom.sig};
1584  std::vector <TH1*> shi_sig_plot = {hplu.sig, hmin.sig};
1585 
1586  std::vector <TH1*> vnom = {hnom.sig, hnom.bkg, hnom.beam, hnom.nc, hnom.numu, hnom.tau};
1587  std::vector <TH1*> vplu = {hplu.sig, hplu.bkg, hplu.beam, hplu.nc, hplu.numu, hplu.tau};
1588  std::vector <TH1*> vmin = {hmin.sig, hmin.bkg, hmin.beam, hmin.nc, hmin.numu, hmin.tau};
1589 
1590  std::vector <TString> labels = {
1591  "$\\nu_e$ signal", "Tot beam bkg", "Beam $\\nu_{e}$ CC",
1592  "NC", "$\\nu_{\\mu}$ CC", "$\\nu_{\\tau}$ CC"};
1593  if(pidaxis == kNumuOpt || pidaxis == kNumuOther)
1594  {
1595  labels = {"$\\nu_{\\mu}$ signal", "Tot beam bkg", "NC",
1596  "$\\nu_{\\mu}$ App", "$\\nu_{e}$ CC", "$\\nu_{\\tau}$ CC"};
1597  vnom = {hnom.sig, hnom.bkg, hnom.nc, hnom.numu, hnom.beam, hnom.tau};
1598  vplu = {hplu.sig, hplu.bkg, hplu.nc, hplu.numu, hnom.beam, hplu.tau};
1599  vmin = {hmin.sig, hmin.bkg, hmin.nc, hmin.numu, hnom.beam, hmin.tau};
1600  }
1601  if(pidaxis == kNumuOpt)
1602  {
1603  nom_bkg_plot = {HistNumuOpt(hnom.bkg), HistNumuOpt(hnom.beam), HistNumuOpt(hnom.nc), HistNumuOpt(hnom.numu), HistNumuOpt(hnom.tau)};
1604  shi_bkg_plot = {HistNumuOpt(hplu.bkg), HistNumuOpt(hmin.bkg)};
1605  nom_sig_plot = {HistNumuOpt(hnom.sig)};
1606  shi_sig_plot = {HistNumuOpt(hplu.sig), HistNumuOpt(hmin.sig)};
1607  }
1608 
1609  PlotMCComponentsErrorBand(nom_bkg_plot, shi_bkg_plot,
1610  leg, pidtag, out_name + "_bkg", printratioplots, false, pidaxis);
1611  PlotMCComponentsErrorBand(nom_sig_plot, shi_sig_plot,
1612  leg2, pidtag, out_name + "_sig", printratioplots, false, pidaxis);
1613 
1614  if(PrintParams) PrintOscilationParameters (calc, usedumb);
1615 
1616  if(printtable){
1617 
1618  ComparisonTable (vnom, vplu, vmin, labels, out_name);
1619 
1620  int nbins = 1;
1621  if (pidaxis== kNue3bin) nbins = 3;
1622  ComparisonTableNbins (vnom, vplu, vmin, labels, out_name, pidaxis);
1623  }
1624  }
1625 
1626 //////////////////////////////////////////////////////////////////////
1627 
1629  PredictionSyst3Flavor2020* predictionSyst, const ISyst* syst,
1630  TString plottitle, TString out_name,
1631  bool printtable=true, AxisType pidaxis=kNue3bin)
1632  {
1633  // Get nominal from PredictionSyst3Flavor2020
1634  auto pred_no = predictionSyst->GetNominalPrediction();
1635 
1636  // Get shifted predictions from PredictionSyst3Flavor2020
1637  const auto& shiftedPred = predictionSyst->GetShiftedPrediction(syst);
1638 
1639  //std::cout << "ShiftedPreds name: " << shiftedPred.systName << std::endl;
1640 
1641  int numSh = shiftedPred.preds.size();
1642  //std::cout << "Number of shifts: " << numSh << std::endl;
1643  //for (auto shift:shiftedPred.preds) { std::cout << "Func shift: " << shift << std::endl; }
1644 
1645  std::vector<IPrediction*> pred_sh = shiftedPred.preds;
1646 
1647  if (numSh == 3)
1648  {
1649  // shift 0 = -1
1650  // shift 1 = Nominal
1651  // shift 2 = +1
1653  pred_no, pred_sh[2], pred_sh[0], plottitle,
1654  out_name+"_1sigma", "#pm1#sigma", printtable, pidaxis);
1655  }
1656  else if (numSh == 5)
1657  {
1658  // shift 1 = -1
1659  // shift 2 = Nominal
1660  // shift 3 = +1
1661  //std::cout << "Passing; " << pred_no << ", " << pred_sh[3] << ", " << pred_sh[1] << std::endl;
1663  pred_no, pred_sh[3], pred_sh[1], plottitle,
1664  out_name+"_1sigma", "#pm1#sigma", printtable, pidaxis);
1665  }
1666  else if (numSh == 4)
1667  {
1668  // fake negative shifts for one-sided systs
1669  // shift 0,1 = 0
1670  // shift 2 = Nominal
1671  // shift 3 = +1
1673  pred_no, pred_sh[3], plottitle, out_name+"_1sigma",
1674  "#plus1#sigma", printtable, pidaxis);
1675  }
1676  else if (numSh == 2)
1677  {
1678  // shift 0 = Nominal
1679  // shift 1 = +1
1680  //std::cout << "Passing; " << pred_no << ", " << pred_sh[1] << std::endl;
1682  pred_no, pred_sh[1], plottitle, out_name+"_1sigma",
1683  "#plus1#sigma", printtable, pidaxis);
1684  }
1685  }
1686 }
struct DataMCComponets GetFDMCComponents(osc::IOscCalc *calc, IPrediction *pred_no, TString output_name="nue", int linestyle=1, bool bkgdetails=false)
const double kLegendDefaultY2
static bool isFHC
tree Draw("slc.nhit")
enum BeamMode kOrange
const double kLegendDefaultY1
enum BeamMode kRed
string output_name
Pickle.
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ComparisonTableOneNbins(const std::vector< TH1 * > &mcnom, const std::vector< TH1 * > &mcshift, const std::vector< TString > &labels, TString ltxcommand, int N=3)
const IDecomp * GetDecomp(IPrediction *prediction, EModExtrapComps modExtrapComp)
virtual Spectrum AntiNueComponent() const =0
double ssth23
osc::IOscCalc * GetCalculator(bool usingdumb)
Antineutrinos-only.
Definition: IPrediction.h:50
Loads shifted spectra from files.
bool mergePeripheral
Float_t y1[n_points_granero]
Definition: compare.C:5
(&#39; appearance&#39;)
Definition: IPrediction.h:18
void CompareOneShiftPred(IPrediction *pred_no, IPrediction *pred_sh, TString plottitle, TString out_name, TString pidtag, bool printtable=true, AxisType pidaxis=kNue1bin, bool PrintParams=false)
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.
Definition: Spectrum.cxx:148
virtual Spectrum NumuComponent() const =0
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
const double kLegendDefaultX2
Float_t x1[n_points_granero]
Definition: compare.C:5
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)
TH1F * h3
Definition: berger.C:36
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:177
(&#39;beam &#39;)
Definition: IPrediction.h:15
const char * p
Definition: xmltok.h:285
const Color_t kNueSignalWSColor
Definition: Style.h:20
TH1 * ratio(TH1 *h1, TH1 *h2)
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
virtual Spectrum Data_Component() const
Definition: IDecomp.cxx:16
TLegend * CustomLegend(std::vector< TH1 * >h, std::vector< TString > title, std::vector< TString > option, double x1=kLegendDefaultX1, double y1=kLegendDefaultY1, double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
void ComparisonTableOne(const std::vector< TH1 * > &mcnom, const std::vector< TH1 * > &mcshift, const std::vector< TString > &labels, const TString ltxcommand="")
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
TLegend * leg1
Definition: plot_hist.C:105
void CompareTwoShiftPred(IPrediction *pred_no, IPrediction *pred_pl, IPrediction *pred_mi, TString plottitle, TString out_name, TString pidtag, bool printtable=true, AxisType pidaxis=kNue1bin, bool PrintParams=false)
TH1F * hplus
Definition: Xsec_final.C:75
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
void PrettyTag(TString pid, int color, double xndc, double yndc)
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const Color_t kNumuBackgroundColor
Definition: Style.h:30
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)
_IOscCalc< double > IOscCalc
Definition: IOscCalc.h:41
double FindLimitY(std::vector< TH1 * > histos, bool doMax=true)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Color_t kNueSignalColor
Definition: Style.h:19
void PlotNDDataTotalMCComparison(TH1 *hdata, std::vector< TH1 * > htots, TLegend *leg, TString pidtag="", TString out_name="plot_nd", bool ratioplot=false, bool pidaxis=false)
TLegend * DefaultNueLegend(bool sig=true, bool tot=true, bool data=false, double x1=kLegendDefaultX1, double y1=kLegendDefaultY1, double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
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, AxisType pidaxis=kNue1bin)
const XML_Char const XML_Char * data
Definition: expat.h:268
c2
Definition: demo5.py:33
const Color_t kNuTauBackgroundColor
Definition: Style.h:26
const int nbins
Definition: cellShifts.C:15
Charged-current interactions.
Definition: IPrediction.h:39
TLegend * DefaultNumuLegend(bool sig=true, bool tot=true, bool data=false, double x1=kLegendDefaultX1, double y1=kLegendDefaultY1, double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
void CompareNDDataMCFromVector(TDirectory *d_no, std::vector< TDirectory * > d_sh, TString plottitle, TString out_name, TString tag, bool printtable=true, AxisType pidaxis=kNue1bin)
const bool printratioplots
bool isNumuAna
void PIDTag(TString pid)
void Nue2018ThreeBinDivisions(bool coreOnly, const int color, const int style)
const DataMCComponents GetNDDecompsFromDecomp(const IDecomp *decomp)
void ComparisonTable(const std::vector< TH1 * > &mcnom, const std::vector< TH1 * > &mcplus, const std::vector< TH1 * > &mcminus, std::vector< TString > labels, TString ltxcommand="")
void FormatErrorBand(TH1 *hplus, TH1 *hminus, bool signal=false, bool fixbins=false)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
const DataMCComponents GetNDCompsFromDecomp(const IDecomp *decomp)
void NuePID3bin(TH1 *h, TString pidtag, bool binlabels=true)
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
#define pot
std::vector< ModularExtrapComponent * > GetModExtrapComponents() const
Definition: ModularExtrap.h:60
TGraph * DummyGraph(int color=kBlack, int style=1, int linewidth=2, int markerstyle=kFullCircle)
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
const double kLegendDefaultX1
std::vector< IPrediction * > preds
void HTag(TString hie)
std::vector< std::string > comps
const ShiftedPreds & GetShiftedPrediction(const ISyst *syst)
void FormatErrorBand(TH1 *hplus, TH1 *hminus, bool signal=false, bool fixbins=false)
const double kAna2020FHCPOT
Definition: Exposures.h:233
virtual Spectrum AntiNumuComponent() const =0
void ComparePredictionsFromVector(TDirectory *dpred_no, std::vector< TDirectory * > dpred_sh, TString plottitle, TString out_name, bool printtable=true, AxisType pidaxis=kNue1bin)
std::vector< float > Spectrum
Definition: Constants.h:610
void CompareNDDataOneMC(TDirectory *d_no, TDirectory *d_sh, TString plottitle, TString out_name, TString pidtag, bool printtable=true, AxisType pidaxis=kNue1bin)
const double kAna2020RHCPOT
Definition: Exposures.h:235
Neutrinos-only.
Definition: IPrediction.h:49
void PrintLatexFigure(TString name, TString plotdir="plots/")
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
virtual _IOscCalc * Copy() const =0
void CompareNDDataTwoMC(TDirectory *d_no, TDirectory *d_pl, TDirectory *d_mi, TString plottitle, TString out_name, TString pidtag, bool printtable=true, AxisType pidaxis=kNue1bin)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
TH1F * h1
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)
TH1 * PrettyRatio(TH1 *shi, TH1 *nom, int color, double alpha=1., TString titley="Shift / Nom")
Interface to extrapolation procedures.
Definition: IExtrap.h:8
AxisType
double dmsq32
const IExtrap * GetExtrap() const
TH1F * hminus
Definition: Xsec_final.C:76
void PrintOscilationParameters(osc::IOscCalc *calc, bool usingdumb)
TH1 * HistNumuOpt(TH1 *orig)
enum BeamMode kViolet
void ComparisonTableNbins(const std::vector< TH1 * > &mcnom, const std::vector< TH1 * > &mcplus, const std::vector< TH1 * > &mcminus, std::vector< TString > labels, TString ltxcommand="", int N=3)
virtual Spectrum NCTotalComponent() const
Definition: IDecomp.h:18
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
Standard interface to all decomposition techniques.
Definition: IDecomp.h:13
void Nue2018ThreeBinLabels(const double yNDC, const double textsize, const int color, const bool nd)
Neutral-current interactions.
Definition: IPrediction.h:40
void PlotNDDataTotalMCComparison(TH1 *hdata, std::vector< TH1 * > htots, TLegend *leg, TString pidtag="", TString out_name="plot_nd", bool ratioplot=false, AxisType pidaxis=kNue1bin)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
c1
Definition: demo5.py:24
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const Color_t kTotalMCColor
Definition: Style.h:16
void DefaultFormatNue(DataMCComponets &comp, int linestyle=1)
Take the output of an extrapolation and oscillate it as required.
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
TLegend * DefaultNueLegend(bool sig=true, bool tot=true, bool data=false, double x1=kLegendDefaultX1, double y1=kLegendDefaultY1, double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
bool usedumb
Extrapolate each component using a separate ModularExtrapComponent.
Definition: ModularExtrap.h:23
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
const Color_t kNCBackgroundColor
Definition: Style.h:22
enum BeamMode kBlue
virtual IOscCalc * Copy() const override
Definition: OscCalcDumb.h:19
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, AxisType pidaxis=kNue1bin)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
TCanvas * RatioPlot(std::vector< TH1 * > topHistos, std::vector< TString > topOption, std::vector< TH1 * > bottomHistos, std::vector< TString > bottomOption, TString pidtag, AxisType pidaxis=kNue1bin)
void FixLegend(TLegend *leg, TString opt="default")
void PrintCanvasAll(TCanvas *c, const TString outname, const TString outdir="plots/", const std::vector< TString > exts={"pdf","png","root"})
const std::string outdir
virtual Spectrum NueComponent() const =0
def sign(x)
Definition: canMan.py:197