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