temp_draw_plots_util.h
Go to the documentation of this file.
1 #pragma once
2 
10 
11 #include "TStyle.h"
12 #include "TFile.h"
13 #include "TH1.h"
14 #include "TCanvas.h"
15 #include "TGraph.h"
16 #include "TGraphAsymmErrors.h"
17 #include "TGaxis.h"
18 #include "TLatex.h"
19 #include "TLegend.h"
20 #include "TLine.h"
21 #include "TSystem.h"
22 
23 #include <iostream>
24 #include <fstream>
25 #include <iomanip>
26 #include <vector>
27 
30 const bool printratioplots=true;
31 
32 // Put a "NOvA Preliminary" tag in the corner
33 void Prelim()
34 {
35  TLatex* prelim = new TLatex(.9, .95, "NOvA Preliminary");
36  prelim->SetTextColor(kBlue);
37  prelim->SetNDC();
38  prelim->SetTextSize(2/36.);
39  prelim->SetTextAlign(32);
40  prelim->Draw();
41 }
42 
43 void CenterTitles(TH1* histo)
44 {
45  histo->GetXaxis()->CenterTitle();
46  histo->GetYaxis()->CenterTitle();
47  histo->GetZaxis()->CenterTitle();
48 }
49 
50 void PIDTag(TString pid)
51 {
52  TLatex* tag = new TLatex(.83, .92, pid);
53  tag->SetTextColor(kAzure+3);
54  tag->SetNDC();
55  tag->SetTextSize(0.038);
56  tag->SetTextFont(42);
57  tag->SetTextAlign(12);
58  tag->Draw();
59 }
60 
61 void HTag(TString pid)
62 {
63  TLatex* tag = new TLatex(.18, .825, pid);
64  if ( pid == "IH") tag->SetTextColor(kOrange+2);
65  if ( pid == "NH") tag->SetTextColor(kAzure+3);
66  tag->SetNDC();
67  tag->SetTextSize(0.043);
68  tag->SetTextFont(42);
69  tag->SetTextAlign(12);
70  tag->Draw();
71 }
72 
73 // NB, loses all the style info, so call early
74 TH1* RebinToShortAxis(TH1* h)
75 {
76  TH1* ret = new TH1D(ana::UniqueName().c_str(), h->GetTitle(),
77  18, 0, 18);
78  ret->GetXaxis()->SetTitle(h->GetXaxis()->GetTitle());
79  ret->GetYaxis()->SetTitle(h->GetYaxis()->GetTitle());
80 
81  std::vector<std::pair<int, int>> binMap = {{3, 2},
82  {4, 3},
83  {5, 4},
84  {6, 5}, // end of 1st spectrum
85  {13, 8},
86  {14, 9},
87  {15, 10},
88  {16, 11}, // end of 2nd
89  {23, 14},
90  {24, 15},
91  {25, 16},
92  {26, 17}};
93 
94  for(auto it: binMap) ret->SetBinContent(it.second, h->GetBinContent(it.first));
95 
96  return ret;
97 }
98 
99 
100  TString MakeLatexCommandName(TString mystring){
101  //There must be a better way
102  std::vector<TString> in = {"0","1","2","3","4","5",
103  "6","7","8","9","_"," ","(",")"};
104  std::vector<TString> out = {"zero","one","two","three","four","five",
105  "six","seven","eight","nine","","","",""};
106  for(unsigned int i=0;i<in.size();i++)
107  mystring.ReplaceAll(in[i],out[i]);
108  return mystring;
109 
110  }
111 
113  // Separators between the different PID ranges
114  // TGraph not TLine because TGraph is properly clipped by frame
115  TGraph* l1 = new TGraph;
116  l1->SetPoint(0, 10, 0);
117  l1->SetPoint(1, 10, 50000);
118  TGraph* l2 = new TGraph;
119  l2->SetPoint(0, 20, 0);
120  l2->SetPoint(1, 20, 50000);
121  for(TGraph* l: {l1, l2}){
122  l->SetLineWidth(1);
123  l->SetLineColor(kGray+3);
124  l->Draw("l same");
125  }
126 }
127 
129  // Put the labels near the top of the plot
130  gPad->Update();
131  const double y = gPad->GetUymax() * 0.955;
132 
133  const ana::Binning* bins = 0;
134  if(pid == "LID") bins = &ana::kLIDNLBinning;
135  if(pid == "LEM") bins = &ana::kLEMNLBinning;
136  if(pid == "CVN") bins = &ana::kCVNNLBinning;
137  assert(bins);
138 
139  for(int i = 0; i < 3; ++i){
140 // TLatex* ltx = new TLatex(3+6*i, y, TString::Format("%g < %s < %g",
141  TLatex* ltx = new TLatex(5+10*i, y, TString::Format("%g < %s < %g",
142  bins->Edges()[i],
143  pid.c_str(),
144  bins->Edges()[i+1]));
145  ltx->SetTextSize(0.032);
146  ltx->SetTextColor(kGray+3);
147  ltx->SetTextAlign(22);
148  ltx->Draw();
149  }
150 
151 }
152 
153 void PID2DAxis(TH1* axes, bool third = false)
154 {
155 
156  double xmax = axes->GetXaxis()->GetXmax();
157  double xmin = axes->GetXaxis()->GetXmin();
158  double xrange = xmax-xmin;
159 
160  double ymin = axes->GetMinimum();
161  TGaxis* gax1 = new TGaxis(0, ymin, 8, ymin, 0, 4, 4, "");
162  TGaxis* gax2 = new TGaxis(10, ymin, 18, ymin, 0, 4, 4, "");
163  TGaxis* gax3 = new TGaxis(20, ymin, 30, ymin, 0, 5, 5, "");
164 
165  float size = axes->GetXaxis()->GetLabelSize();
166  float offs = axes->GetXaxis()->GetLabelOffset();
167  int font = axes->GetXaxis()->GetLabelFont();
168  // And then hide its title
169  axes->GetXaxis()->SetLabelSize(0);
170  // and make all the tick marks the same size
171  axes->GetXaxis()->SetNdivisions(3000);
172  axes->GetXaxis()->SetTitle("Reconstructed Neutrino Energy (GeV)");
173 // axes->GetYaxis()->SetTitle("Events / 6.05 x 10^{20} POT-equiv");
174  axes->GetYaxis()->SetTitleOffset(0.8);
175  axes->SetTitle("");
176  CenterTitles(axes);
177 
178  if ( third ){
179  gax1->SetLabelSize(size);
180  gax1->SetLabelOffset(offs);
181  gax1->SetLabelFont(font);
182  gax1->Draw();
183  }
184  else
185  for(TGaxis* ax: {gax1, gax2, gax3}){
186  // Copy settings from pre-existing axis
187  ax->SetLabelSize(size-0.001);
188  ax->SetLabelOffset(offs);
189  ax->SetLabelFont(font);
190  ax->Draw();
191  }
192 
193  TLine* xline = new TLine();
194  xline->SetLineWidth(2);
195  xline->DrawLine(xmin, 0, xmax, 0);
196  gPad->RedrawAxis();
197 }
198 //----------------------------------------------------------------------
200 {
201  // Separators between the different PID ranges
202  // TGraph not TLine because TGraph is properly clipped by frame
203  TGraph* l1 = new TGraph;
204  l1->SetPoint(0, 6, 0);
205  l1->SetPoint(1, 6, 50000);
206  TGraph* l2 = new TGraph;
207  l2->SetPoint(0, 12, 0);
208  l2->SetPoint(1, 12, 50000);
209  for(TGraph* l: {l1, l2}){
210  // l->SetLineStyle(3);
211  l->SetLineWidth(1);
212  l->SetLineColor(kGray+3);
213  l->Draw("l same");
214  }
215 
216  // Hardcoding maxy is obviously unsustainable
217  TGaxis* gax1 = new TGaxis(6, 0, 6, 18, 0, 18, 504, "U+-");
218  TGaxis* gax2 = new TGaxis(12, 0, 12, 18, 0, 18, 504, "U+-");
219 }
220 //----------------------------------------------------------------------
222 {
223  const std::string fPID = "CVN"; // HACK
224 
225  // Put the labels near the top of the plot
226  std::string pid = (std::string)fPID;
227  gPad->Update();
228  const double y = gPad->GetUymax() * .94;
229 
230  const ana::Binning* bins = 0;
231  if(fPID == "LID") bins = &ana::kLIDNLBinning;
232  if(fPID == "LEM") bins = &ana::kLEMNLBinning;
233  if(fPID == "CVN") bins = &ana::kCVNNLBinning;
234  assert(bins);
235 
236  for(int i = 0; i < 3; ++i){
237  TLatex* ltx = new TLatex(3+6*i, y, TString::Format("%g < %s < %g",
238  bins->Edges()[i],
239  pid.c_str(),
240  bins->Edges()[i+1]));
241  ltx->SetTextSize(.04);//0.035);
242  ltx->SetTextAlign(22);
243  ltx->Draw();
244  }
245 }
247 {
248  double ymin = hist->GetMinimum();
249  // Paint our own x-axis fragments
250  TGaxis* gax1 = new TGaxis(1, ymin, 5, ymin, 1, 3, 3, "");
251  TGaxis* gax2 = new TGaxis(7, ymin, 11, ymin, 1, 3, 3, "");
252  TGaxis* gax3 = new TGaxis(13, ymin, 17, ymin, 1, 3, 3, "");
253 
254  float tick = hist->GetXaxis()->GetTickLength();
255  float size = hist->GetXaxis()->GetLabelSize();
256  float offs = hist->GetXaxis()->GetLabelOffset();
257  int font = hist->GetXaxis()->GetLabelFont();
258  // And then hide its title
259  hist->GetXaxis()->SetLabelSize(0);
260  // and make all the tick marks the same size
261  hist->GetXaxis()->SetNdivisions(1800);
262  hist->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
263  hist->GetYaxis()->SetTitle("Events / 0.5 GeV Bin");
264  hist->GetYaxis()->SetTitleOffset(0.8);
265  hist->GetXaxis()->CenterTitle();
266  hist->GetYaxis()->CenterTitle();
267  hist->SetTitle("");
268 
269  if (false){// fThird ){
270  gax1->SetTickSize(tick);
271  gax1->SetLabelSize(size);
272  gax1->SetLabelOffset(offs);
273  gax1->SetLabelFont(font);
274  gax1->Draw();
275  }
276  else{
277  for(TGaxis* ax: {gax1, gax2, gax3}){
278  // Copy settings from pre-existing axis
279  ax->SetTickSize(tick);
280  ax->SetLabelSize(size);
281  ax->SetLabelOffset(offs);
282  ax->SetLabelFont(font);
283  ax->Draw();
284  }
285  }
286 }
287 
288 double FindLimitY (std::vector<TH1*> histos, bool doMax = true)
289 {
290  double max = histos[0]->GetMaximum();
291  double min = histos[0]->GetMinimum();
292 
293  for(unsigned int i=0; i<histos.size(); i++){
294  if ( histos[i]->GetMaximum() > max ) max = histos[i]->GetMaximum();
295  if ( histos[i]->GetMinimum() < min ) min = histos[i]->GetMinimum();
296  }
297 
298  if ( doMax ) return max;
299  else return min;
300 
301 }
302 
303 namespace ana
304 {
305 
306  void ComparisonTable(std::vector<TH1*> mcnom, std::vector<TH1*> mcshift, std::vector<TString> labels
307  ,TString ltxcommand = ""){
308  unsigned int nentries = mcnom.size();
309  if(mcshift.size()!=nentries || labels.size()!=nentries){std::cerr << "incompatible sizes\n"; return;}
310 
311  if(ltxcommand != "") std::cout << "\n\\newcommand{\\" << MakeLatexCommandName(ltxcommand) << "}\n";
312  std::cout << "{\\begin{tabular}{ l | r r r } \n"
313  << std::setw(20) <<" " << " & "
314  << std::setw(10) << "Nominal" << " & "
315  << std::setw(10) << "Shift" << " & "
316  << std::setw(12) << "\\% Diff."
317  <<"\t \\\\ \\hline \n";
318  for(unsigned int i=0;i<nentries;i++){
319  std::cout.unsetf(std::ios_base::floatfield);
320  std::cout << std::setprecision(3)
321  << std::setw(20) << labels[i] << " & "
322  << std::setw(10) << mcnom[i] ->Integral() << " & "
323  << std::setw(10) << mcshift[i]->Integral() << " & "
324  << std::setw(10) << std::setprecision(1) << std::fixed
325  << ( mcnom[i]->Integral()>0 ? (((mcshift[i]->Integral()/mcnom[i]->Integral()) -1.)*100.):0.)
326  << "\\% \t \\\\ \n";
327  }
328  std::cout << "\\end{tabular}\n}\n";
329  }
330 
331 
332  void ComparisonTable(std::vector<TH1*> mcnom, std::vector<TH1*> mcplus, std::vector<TH1*> mcminus,
333  std::vector<TString> labels,TString ltxcommand = ""){
334  unsigned int nentries = mcnom.size();
335  if(mcplus.size()!=nentries || mcminus.size()!=nentries || labels.size()!=nentries){std::cerr << "incompatible sizes\n"; return;}
336  if(ltxcommand != "") std::cout << "\\newcommand{\\" << MakeLatexCommandName(ltxcommand) << "}\n";
337  std::cout << "{ \\begin{tabular}{ l | r r r | r r } \n"
338  << std::setw(20) <<" " << " & "
339  << std::setw(10) << "Nominal" << " & "
340  << std::setw(10) << "Shift (+)" << " & "
341  << std::setw(10) << "Shift (-)" << " & "
342  << std::setw(12) << "\\% Diff. (+)" << " & "
343  << std::setw(12) << "\\% Diff. (-)"
344  <<"\t \\\\ \\hline \n";
345  for(unsigned int i=0;i<nentries;i++){
346  std::cout.unsetf(std::ios_base::floatfield);
347  std::cout << std::setprecision(3)
348  << std::setw(20) << labels[i] << " & "
349  << std::setw(10) << mcnom[i] ->Integral() << " & "
350  << std::setw(10) << mcplus[i] ->Integral() << " & "
351  << std::setw(10) << mcminus[i]->Integral() << " & "
352  << std::setw(10) << std::setprecision(1) << std::fixed
353  << ( mcnom[i]->Integral()>0 ? (((mcplus [i]->Integral()/mcnom[i]->Integral()) -1.)*100.):0.)
354  << "\\% & "
355  << std::setw(10)
356  << ( mcnom[i]->Integral()>0 ? (((mcminus[i]->Integral()/mcnom[i]->Integral()) -1.)*100.):0.)
357  << "\\% \t \\\\ \n";
358  }
359  std::cout << "\\end{tabular} } \n";
360  }
361 
362  ////////////////////////////////////////
363 
364  const double kLegendDefaultX1 = 0.63;
365  const double kLegendDefaultX2 = 0.85;
366  const double kLegendDefaultY1 = 0.68;
367  const double kLegendDefaultY2 = 0.88;
368 
369  void FixLegend(TLegend * leg, TString opt="default"){
370  leg->SetFillStyle(0);
371  if(opt=="default"){
372  leg->SetX1(kLegendDefaultX1); leg->SetX2(kLegendDefaultX2);
373  leg->SetY1(kLegendDefaultY1); leg->SetY2(kLegendDefaultY2);
374  }
375  if(opt=="ExPID"){
376  leg->SetX1(0.405); leg->SetX2(0.565);
377  leg->SetY1(0.59); leg->SetY2(0.82);
378  leg->SetBorderSize(0);
379  leg->SetTextSize(0.037);
380  }
381  if(opt=="ExPIDPads"){
382  leg->SetX1(0.55); leg->SetX2(0.89);
383  leg->SetY1(0.75); leg->SetY2(0.89);
384  leg->SetTextSize(0.045);
385  }
386  }
387  ////////////////////////////////////////
388  TLegend * DefaultNueLegend (bool sig=true, bool tot=true,bool data=false,
389  double x1=kLegendDefaultX1, double y1=kLegendDefaultY1,
390  double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
391  {
392  TLegend* leg = new TLegend(x1,y1,x2,y2);
393  leg->SetFillStyle(0);
394  TH1* hdata = new TH1F("", "", 1, 0, 1);
395  TH1* hsig = new TH1F("", "", 1, 0, 1);
396  TH1* hNC = new TH1F("", "", 1, 0, 1);
397  TH1* hCC = new TH1F("", "", 1, 0, 1);
398  TH1* hBeam = new TH1F("", "", 1, 0, 1);
399  TH1* htot = new TH1F("", "", 1, 0, 1);
400  hdata->SetMarkerStyle(kFullCircle);
401  hsig->SetLineColor(kNueSignalColor);
402  hNC->SetLineColor(kNCBackgroundColor);
403  hCC->SetLineColor(kNumuBackgroundColor);
404  hBeam->SetLineColor(kBeamNueBackgroundColor);
405  htot->SetLineColor(kTotalMCColor);
406  if(sig) leg->AddEntry(hsig,"#nu_{e} Signal","l");
407  if(data) leg->AddEntry(hdata,"ND Data","epl");
408  if(tot)leg->AddEntry(htot,"Total Bkg","l");
409  leg->AddEntry(hBeam,"Beam #nu_{e}","l");
410  leg->AddEntry(hNC,"NC","l");
411  leg->AddEntry(hCC,"#nu_{#mu} CC","l");
412  return leg;
413 
414  }
415  ////////////////////////////////////////
416  TLegend * CustomLegend (std::vector<TH1*>h, std::vector<bool>isdata,
417  std::vector<TString>title,
418  double x1=kLegendDefaultX1, double y1=kLegendDefaultY1,
419  double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
420  {
421  TLegend* leg = new TLegend(x1,y1,x2,y2);
422  leg->SetFillStyle(0);
423  for ( unsigned int v = 0; v<h.size(); ++v){
424  if (isdata[v]) TLegendEntry* le = leg->AddEntry(h[v],title[v],"epl");
425  else{
426  if (h[v]->GetLineColor() == kRed) TLegendEntry* le = leg->AddEntry(h[v],title[v],"l");
427  else TLegendEntry* le = leg->AddEntry(h[v],title[v],"f");
428  }
429  }
430  return leg;
431 
432  }
433  ////////////////////////////////////////
434  TCanvas * ExPIDPlot (std::vector<TH1*> topHistos,
435  std::vector<TString> topOption)
436  {
437  TCanvas *cpid = new TCanvas("cpid", "cpid", 410, 850);
438  float Ebins = 10;
439  int pidBins = 3;
440 
441  TPad* p[pidBins];
442  p[0] = new TPad("", "", 0, 0, 1, 1);
443  p[1] = new TPad("", "", 0, 0, 1, 1);
444  p[2] = new TPad("", "", 0, 0, 1, 1);
445  p[0]->SetTopMargin(0.1);
446  p[0]->SetBottomMargin(.634);
447  p[1]->SetTopMargin(.369);
448  p[1]->SetBottomMargin(.369);
449  p[2]->SetTopMargin(.634);
450  p[2]->SetBottomMargin(0.1);
451  p[0]->SetFillStyle(0);
452  p[1]->SetFillStyle(0);
453  p[2]->SetFillStyle(0);
454  p[0]->Draw();
455  p[1]->Draw();
456  p[2]->Draw();
457 
458  TH1* he[pidBins];
459  TH1* h[topHistos.size()][pidBins];
460  for(int k=pidBins-1; k>=0; --k){
461  he[k] = (TH1*) topHistos[0]->Clone();
462  he[k]->GetYaxis()->SetRangeUser(FindLimitY(topHistos,0) ,1.2*FindLimitY(topHistos,1));
463  he[k]->GetXaxis()->SetRangeUser(Ebins*k,Ebins*(k+1));
464  if ( k>0 ) he[k]->GetXaxis()->SetLabelSize(0);
465  p[pidBins-1-k]->cd();
466  he[k]->Draw(topOption[0]);
467  he[k]->GetYaxis()->SetTitleSize(26);
468  he[k]->GetYaxis()->SetTitleFont(43);
469  he[k]->GetYaxis()->SetTitleOffset(1.25);
470  he[k]->GetYaxis()->SetLabelFont(43);
471  he[k]->GetYaxis()->SetLabelSize(14);
472  he[k]->GetXaxis()->SetLabelOffset(-0.008);
473  he[k]->GetXaxis()->SetTitleOffset(0.5);
474  p[pidBins-1-k]->SetGridy();
475  }
476  for(unsigned int ii=0;ii<topHistos.size();ii++){
477  for(int j=pidBins-1; j>=0; --j){
478  h[ii][j] = (TH1*) topHistos[ii]->Clone();
479  h[ii][j]->GetXaxis()->SetRangeUser(Ebins*j,Ebins*(j+1));
480  p[pidBins-1-j]->cd();
481  h[ii][j]->Draw("same "+topOption[ii]);
482  }
483  }
484  for(int h=pidBins-1; h>=0; --h){
485  he[h]->SetStats(0);
486  p[pidBins-1-h]->cd();
487  he[h]->Draw("same "+topOption[0]);
488  }
489  p[0]->RedrawAxis(); //needed for white histograms e.g. beam
490  p[1]->RedrawAxis();
491  p[2]->RedrawAxis();
492  return cpid;
493  }
494  ////////////////////////////////////////
495  TCanvas * RatioPlot (std::vector<TH1*> topHistos,
496  std::vector<TString> topOption,
497  std::vector<TH1*> bottomHistos,
498  std::vector<TString> bottomOption,
499  TString pidtag,
500  bool pidaxis=false)
501  {
502  TCanvas *c = new TCanvas("c", "canvas", 600+(150*pidaxis), 700-(50*pidaxis));
503  auto h1 = (TH1*) topHistos[0]->Clone();
504  auto h3 = (TH1*) bottomHistos[0]->Clone();
505 
506  CenterTitles(h1);
507  CenterTitles(h3);
508 
509  TPad* p1 = new TPad("", "", 0, 0, 1, 1);
510  TPad* p2 = new TPad("", "", 0, 0, 1, 1);
511  p1->SetBottomMargin(.3);
512  p2->SetTopMargin(.7);
513  p1->SetFillStyle(0);
514  p2->SetFillStyle(0);
515  p1->Draw();
516  p2->Draw();
517 
518  p1->cd();
519  h1->SetStats(0);
520  h1->Draw(topOption[0]);
521  for(unsigned int ii=0;ii<topHistos.size();ii++){
522  topHistos[ii]->Draw("same "+topOption[ii]);}
523  h1->Draw("same "+topOption[0]);
524  p1->RedrawAxis(); //white histograms cover axes
525  h1->GetYaxis()->SetRangeUser(FindLimitY(topHistos,0) ,1.2*FindLimitY(topHistos,1));
526  h1->GetYaxis()->SetTitleSize(26);
527  h1->GetYaxis()->SetTitleFont(43);
528  h1->GetYaxis()->SetTitleOffset(1.2);
529  h1->GetYaxis()->SetLabelFont(43);
530  h1->GetYaxis()->SetLabelSize(18);
531  h1->GetXaxis()->SetLabelSize(0);
532  h1->GetXaxis()->SetTitleSize(0);
533  if (pidaxis){
534  PIDDivisions();
535  PIDBinLabels((std::string)pidtag);
536  }
537 
538  p2->cd();
539  h3->SetTitle("");
540  h3->GetYaxis()->SetRangeUser(max(0.,2-FindLimitY(bottomHistos,1)) ,1.01*FindLimitY(bottomHistos,1));
541  h3->GetYaxis()->SetTitleSize(26);
542  h3->GetYaxis()->SetTitleFont(43);
543  h3->GetYaxis()->SetTitleOffset(1.2);
544  h3->GetYaxis()->SetLabelFont(43);
545  h3->GetYaxis()->SetLabelSize(15);
546  h3->GetYaxis()->SetDecimals();
547 
548  h3->GetXaxis()->SetTitle(h1->GetXaxis()->GetTitle());
549  h3->GetXaxis()->SetTitleSize(26);
550  h3->GetXaxis()->SetTitleFont(43);
551  h3->GetXaxis()->SetLabelFont(43);
552  h3->GetXaxis()->SetLabelSize(18);
553  h3->GetXaxis()->SetLabelOffset(0.005);
554  h3->Draw(bottomOption[0]);
555  for(unsigned int ii=0;ii<bottomHistos.size();ii++){
556  bottomHistos[ii]->Draw("same "+bottomOption[ii]);}
557  h3->Draw("same "+bottomOption[0]);
558 
559  auto lone = new TLine();
560  lone->SetLineStyle(3);
561  p2->Update();
562  lone->DrawLine(p2->GetUxmin(),1.0,p2->GetUxmax(),1.0);
563  lone->DrawLine(p2->GetUxmin(),1.2,p2->GetUxmax(),1.2);
564  lone->DrawLine(p2->GetUxmin(),0.8,p2->GetUxmax(),0.8);
565  float p2Y = p2->GetUymax()-p2->GetUymin();
566  h3->GetYaxis()->SetRangeUser(p2->GetUymin()-(0.01*p2Y),p2->GetUymax()+(0.01*p2Y));
567 
568  if (pidaxis){
569  PIDDivisions();
570  PID2DAxis(h3);
571  }
572  return c;
573  }
574  ////////////////////////////////////////
575  void PlotDataMC(std::vector<TH1*>vnom, std::vector<bool>isdata,
576  TLegend * leg, TString pidtag="",TString htag="",
577  TString out_name="plot_FD",bool ratioplot=false,
578  bool ratioerrors=false, bool pidaxis = 0)
579  {
580  gStyle->SetPadLeftMargin(0.15);
581  gStyle->SetPadTopMargin(0.08);
582 
583  // Fill empty bins with 0 if in range and <0 if out of range
584  for ( unsigned int bin = 1; bin<=vnom[0]->GetNbinsX(); ++bin){
585  bool filledbin=false;
586  for ( unsigned int v = 0; v<vnom.size(); ++v){
587  if ( !isdata[v] )
588  if (vnom[v]->GetBinContent(bin) > 0 ) filledbin=true;
589  }
590  for ( unsigned int v = 0; v<vnom.size(); ++v){
591  if (vnom[v]->GetBinContent(bin) == 0 ){
592  if ( !filledbin ) vnom[v]->SetBinContent(bin,-0.01);
593  else vnom[v]->SetBinContent(bin,0.001);
594  }
595  vnom[v]->GetYaxis()->SetTitleOffset(1.05);
596  }
597  }
598 
599  // Draw marker w/error or line when appropriate
600  int datahist=-1;
601  for ( unsigned int v = 0; v<vnom.size(); ++v)
602  if (vnom[v]->GetLineColor() == kBlack ) datahist=v;
603 
604  // Option to make a ratio plot given
605  if(!ratioplot){
606  auto c2 = new TCanvas ("c2","c2",715,500);
607  for ( unsigned int v = 0; v<vnom.size(); ++v){
608  TGraphAsymmErrors* gr = GraphWithPoissonErrors(vnom[v], 0, 1);
609  if ( isdata[v] ){
610  if (v == datahist ) gr->Draw("ep same");
611  else vnom[v]->Draw("p same");
612  }
613  else vnom[v]->Draw("hist same");
614  }
615  gPad->RedrawAxis();
616 
617  // Adjust for ExPID axis
618  if (pidaxis){
619  FixLegend(leg,"ExPID");
620  NuePlotStyle style(pidtag,0);
621  style.PIDDivisions(true);
622  style.PID2DAxis(vnom[0],true);
623  style.PIDBinLabels(true);
624  }
625 
626  // pot label
627  TLatex* ltx = new TLatex(.175, 0.74,
628  TString::Format("%g#times10^{20} POT equiv.", 6.05).Data());
629  ltx->SetNDC();
630  ltx->SetTextAlign(11);
631  ltx->SetTextSize(0.0365);
632 // ltx->Draw();
633  // HTag(htag);
634  leg->Draw("same");
635 
636  // Save with & without prelim tag
637  gPad->SetFillStyle(0);
638  c2->SetFillStyle(0);
639  c2->SaveAs("plots/"+out_name+"_datamc.pdf");
640  c2->SaveAs("plots/"+out_name+"_datamc.eps");
641  c2->SaveAs("plots/"+out_name+"_datamc.C");
642  c2->SaveAs("plots/"+out_name+"_datamc.root");
643 // Prelim();
644  // c2->SaveAs("plots/"+out_name+"_pre_datamc.pdf");
645  //c2->SaveAs("plots/"+out_name+"_pre_datamc.eps");
646  //c2->SaveAs("plots/"+out_name+"_pre_datamc.C");
647  //c2->SaveAs("plots/"+out_name+"_pre_datamc.root");
648  }
649  else{
650  std::vector<TH1*> vratio;
651  for(unsigned int ii=1;ii<vnom.size();ii++){
652  auto htemp = (TH1*) vnom[ii]->Clone(ana::UniqueName().c_str());
653  htemp->Divide(vnom[datahist],vnom[ii]);
654  htemp->SetMarkerColor(vnom[ii]->GetLineColor());
655  htemp->SetMarkerStyle(34);
656  htemp->SetMarkerSize(1.2);
657  htemp->GetYaxis()->SetTitle("Data / MC");
658  htemp->GetYaxis()->SetRangeUser(0.5,1.5);
659  vratio.push_back(htemp);
660  }
661  std::vector<TString> topOption = {" ","hist"};
662  std::vector<TString> bottomOption (vratio.size(),ratioerrors ?"p":"hist p");
663  auto cratio2= RatioPlot(vnom, topOption,vratio,bottomOption,pidtag,pidaxis);
664  cratio2->cd();
665  leg->Draw("same");
666  //HTag(pidtag);
667  cratio2->SaveAs("plots/"+out_name+"_datamc_ratio.pdf");
668  cratio2->SaveAs("plots/png_"+out_name+"_datamc_ratio.png");
669  cratio2->SaveAs("plots/rootfiles/"+out_name+"_datamc_ratio.root");
670  }
671  }
672 
673  ////////////////////////////////////////
674  void PlotNDDataTotalMCComparison(TH1* hdata,std::vector<TH1*> htots,
675  TLegend * leg,TString pidtag="",
676  TString out_name="plot_nd",
677  bool ratioplot=false, bool pidaxis=false)
678  {
679  if(pidaxis) FixLegend(leg,"ExPID");
680  std::vector<TString> h_opt;
681  for(auto h:htots) h_opt.push_back("hist");
682  h_opt.push_back("p");
683 
684  if(!ratioplot){
685  auto c1 = new TCanvas ("c1","c1");
686  htots[0]->SetMaximum(1.3*htots[0]->GetMaximum());
687  for(auto h:htots) h->Draw("hist same");
688  htots[0]->Draw("hist same"); //need to redraw for error bands
689  hdata->Draw("p same");
690  PIDTag(pidtag);
691  if (pidaxis){
692  PIDDivisions();
693  PIDBinLabels((std::string)pidtag);
694  PID2DAxis(htots[0]);
695  }
696  leg->Draw("same");
697  c1->SaveAs("plots/"+out_name+"_datamc.pdf");
698  c1->SaveAs("plots/png_"+out_name+"_datamc.png");
699  c1->SaveAs("plots/rootfiles/"+out_name+"_datamc.root");
700  }
701  else{
702  std::vector<TH1*> ratios;
703  std::vector<TString> r_opt;
704  for(auto h:htots){
705  auto ratio = (TH1*) h->Clone();
706  ratio->Divide(hdata,h);
707  ratio->GetYaxis()->SetRangeUser(0.5,1.5);
708  ratios.push_back(ratio);
709  r_opt.push_back("e1");
710  }
711  ratios[0]->GetYaxis()->SetTitle("Data / MC");
712  ratios[0]->GetYaxis()->SetRangeUser(0.5,1.5);
713  htots.push_back(hdata);
714  auto cratio1 =RatioPlot(htots,h_opt,ratios,r_opt,pidtag,pidaxis);
715  cratio1->cd();
716  leg->Draw("same");
717  PIDTag(pidtag);
718  cratio1->SaveAs("plots/"+out_name+"_datamc_ratio.pdf");
719  cratio1->SaveAs("plots/png_"+out_name+"_datamc_ratio.png");
720  cratio1->SaveAs("plots/rootfiles/"+out_name+"_datamc_ratio.root");
721  }
722  if(pidaxis){
723  htots.push_back(hdata);
724  auto cpads =ExPIDPlot(htots,h_opt);
725  cpads->cd();
726  FixLegend(leg,"ExPIDPads");
727  leg->Draw("same");
728  PIDTag(pidtag);
729 // cpads->SaveAs("plots/"+out_name+"_datamc_pads.pdf");
730  // cpads->SaveAs("plots/png_"+out_name+"_datamc_pads.png");
731  //cpads->SaveAs("plots/rootfiles/"+out_name+"_datamc_pads.root");
732  }
733  }
734 
735  ////////////////////////////////////////
736  void PlotMCComponentsComparison(std::vector<TH1*>vnom,std::vector<TH1*>vshi,
737  TLegend * leg, TString pidtag="",
738  TString out_name="plot_nd",bool ratioplot=false,
739  bool ratioerrors=false, bool pidaxis = 0)
740  {
741 
742  if(pidaxis) FixLegend(leg,"ExPID");
743  if(!ratioplot){
744  auto c2 = new TCanvas ("c2","c2");
745  for(auto v:vnom) v->Draw("hist same");
746  for(auto v:vshi) v->Draw("hist same");
747  if (pidaxis){
748  PIDDivisions();
749  PIDBinLabels((std::string)pidtag);
750  PID2DAxis(vnom[0]);
751  }
752  PIDTag(pidtag);
753  leg->Draw("same");
754  c2->SaveAs("plots/"+out_name+"_mccomp.pdf");
755  c2->SaveAs("plots/png_"+out_name+"_mccomp.png");
756  c2->SaveAs("plots/rootfiles/"+out_name+"_mccomp.root");
757  c2->Close();
758  }
759  else{
760  std::vector<TH1*> vratio;
761  for(unsigned int ii=0;ii<vshi.size();ii++){
762  auto htemp = (TH1*) vnom[ii]->Clone(ana::UniqueName().c_str());
763  htemp->Divide(vshi[ii],vnom[ii]);
764  htemp->SetMarkerColor(vshi[ii]->GetLineColor());
765  htemp->SetMarkerStyle(34);
766  htemp->SetMarkerSize(1.2);
767  htemp->GetYaxis()->SetTitle("Shift / Nom");
768  htemp->GetYaxis()->SetRangeUser(0.5,1.5);
769  vratio.push_back(htemp);
770  }
771  vnom.insert(vnom.end(),vshi.begin(),vshi.end());
772  std::vector<TString> topOption (vnom.size(),"hist");
773  std::vector<TString> bottomOption (vratio.size(),ratioerrors ?"p":"hist p");
774  auto cratio2= RatioPlot(vnom, topOption,vratio,bottomOption,pidtag,pidaxis);
775  cratio2->cd();
776  leg->Draw("same");
777  PIDTag(pidtag);
778  cratio2->SaveAs("plots/"+out_name+"_mccomp_ratio.pdf");
779  cratio2->SaveAs("plots/png_"+out_name+"_mccomp_ratio.png");
780  cratio2->SaveAs("plots/rootfiles/"+out_name+"_mccomp_ratio.root");
781  }
782  if(pidaxis){
783  vnom.insert(vnom.end(),vshi.begin(),vshi.end());
784  std::vector<TString> topOption (vnom.size(),"hist");
785  auto cpads =ExPIDPlot(vnom,topOption);
786  cpads->cd();
787  FixLegend(leg,"ExPIDPads");
788  leg->Draw("same");
789  PIDTag(pidtag);
790 // cpads->SaveAs("plots/"+out_name+"_mccomp_pads.pdf");
791  // cpads->SaveAs("plots/png_"+out_name+"_mccomp_pads.png");
792  //cpads->SaveAs("plots/rootfiles/"+out_name+"_mccomp_pads.root");
793 
794  }
795  }
796 
797  void CompareNDDataMC (TDirectory * d_no, TDirectory* d_sh,TString plottitle,TString out_name, TString pidtag, bool printtable=true, bool pidaxis = 0){
798  auto data = Spectrum::LoadFrom(d_no->GetDirectory("data_comp"));
799  auto kNDPOT = data->POT();
800  auto hdata = data->ToTH1(kNDPOT);
801  hdata ->SetMarkerStyle(kFullCircle);
802 
803  auto hnue_no = (*Spectrum::LoadFrom(d_no->GetDirectory("nue_comp")) +
804  *Spectrum::LoadFrom(d_no->GetDirectory("antinue_comp"))).ToTH1(kNDPOT, kBeamNueBackgroundColor, 1);
805  auto hnumu_no = (*Spectrum::LoadFrom(d_no->GetDirectory("numu_comp")) +
806  *Spectrum::LoadFrom(d_no->GetDirectory("antinumu_comp"))).ToTH1(kNDPOT, kNumuBackgroundColor, 1);
807  auto hnc_no = (*Spectrum::LoadFrom(d_no->GetDirectory("nc_tot_comp"))).ToTH1(kNDPOT, kNCBackgroundColor, 1);
808  auto htot_no = (TH1*)hnue_no->Clone();
809  htot_no->Add(hnumu_no);htot_no->Add(hnc_no);
810 
811  auto hnue_sh = (*Spectrum::LoadFrom(d_sh->GetDirectory("nue_comp")) +
812  *Spectrum::LoadFrom(d_sh->GetDirectory("antinue_comp"))).ToTH1(kNDPOT, kBeamNueBackgroundColor, 2);
813  auto hnumu_sh = (*Spectrum::LoadFrom(d_sh->GetDirectory("numu_comp")) +
814  *Spectrum::LoadFrom(d_sh->GetDirectory("antinumu_comp"))).ToTH1(kNDPOT, kNumuBackgroundColor, 2);
815  auto hnc_sh = (*Spectrum::LoadFrom(d_sh->GetDirectory("nc_tot_comp"))).ToTH1(kNDPOT, kNCBackgroundColor, 2);
816  auto htot_sh = (TH1*)hnue_sh->Clone();
817  htot_sh->Add(hnumu_sh);htot_sh->Add(hnc_sh);
818 
819  htot_no->SetLineColor(kTotalMCColor+1);
820  htot_sh->SetLineColor(kTotalMCColor+2);
821 
822  htot_no->SetMaximum(FindLimitY({htot_no,hdata}));
823  htot_no->SetTitle(plottitle);
824  hnumu_no->SetMaximum(FindLimitY({hnumu_no,hnc_no,hnue_no,hnumu_sh,hnc_sh,hnue_sh}));
825  hnumu_no->SetTitle(plottitle);
826 
827  CenterTitles(htot_no);
828  CenterTitles(hnue_no);
829  TGaxis::SetMaxDigits(3);
830  htot_no->GetXaxis()->SetDecimals();
831 
832  auto leg1=new TLegend(0.63,0.67,0.85,0.87);
833  leg1->SetFillStyle(0);
834  leg1->AddEntry(hdata, "ND Data","epl");
835  leg1->AddEntry(htot_no,"Nominal","l");
836  leg1->AddEntry(htot_sh,"Shift","l");
837 
838  PlotNDDataTotalMCComparison(hdata,{htot_no,htot_sh},leg1,pidtag,out_name,printratioplots,pidaxis);
839 
840  htot_sh->SetLineColor(kTotalMCColor+1);
841  std::vector<TH1*>vnom={hnumu_no,hnc_no,hnue_no};
842  std::vector<TH1*>vshi={hnumu_sh,hnc_sh,hnue_sh};
843 
844  auto leg2 = DefaultNueLegend(false,false,false);
845  TH1* htemp = new TH1F("", "", 1, 0, 1);
846  htemp->SetLineColor(kGray+1);
847  htemp->SetLineStyle(2);
848  leg2->AddEntry(htemp,"Shifted","l");
849 
850  PlotMCComponentsComparison(vnom,vshi,leg2,pidtag,out_name,printratioplots,false,pidaxis);
851 
852  if(printtable){
853  std::cout << "%" <<std::string(60,'-') << std::endl;
854  std:: cout << "% - " << (plottitle.Remove(plottitle.First(";"))).Data() << " - ND - " << pidtag << " - " << std::endl;
855  ComparisonTable ({htot_no,hnue_no,hnc_no,hnumu_no},
856  {htot_sh,hnue_sh,hnc_sh,hnumu_sh},
857  {"Total MC","Beam $\\nu_{e}$","NC", "$\\nu_{\\mu}$ CC"},
858  (plottitle.Remove(plottitle.First(";")))+pidtag
859  );
860  std::cout << "%" <<std::string(60,'-') << std::endl;
861  }
862 
863 
864  }
865 
866 
867  void ComparePredictions(IPrediction * pred_no, IPrediction * pred_sh, TString plottitle, TString out_name, TString pidtag, bool printtable=true, bool pidaxis = false){
868 
869  TH1* hNueSig_no = pred_no->PredictComponent(&calc, ana::Flavors::kNuMuToNuE,
871  TH1* hNueBeam_no = pred_no->PredictComponent(&calc, ana::Flavors::kNuEToNuE,
873  TH1* hNumu_no = pred_no->PredictComponent(&calc, ana::Flavors::kAllNuMu,
875  TH1* hNC_no = pred_no->PredictComponent(&calc, ana::Flavors::kAll,
877  TH1* hTau_no = pred_no->PredictComponent(&calc, ana::Flavors::kAllNuTau,
879 
880  TH1* hNueSig_sh = pred_sh->PredictComponent(&calc, ana::Flavors::kNuMuToNuE,
882  TH1* hNueBeam_sh = pred_sh->PredictComponent(&calc, ana::Flavors::kNuEToNuE,
884  TH1* hNumu_sh = pred_sh->PredictComponent(&calc, ana::Flavors::kAllNuMu,
886  TH1* hNC_sh = pred_sh->PredictComponent(&calc, ana::Flavors::kAll,
888  TH1* hTau_sh = pred_sh->PredictComponent(&calc, ana::Flavors::kAllNuTau,
890 
891 // std::cout<<"Making plots"<<std::endl;
892  hNueBeam_no->SetLineColor(kBeamNueBackgroundColor);
893  hNumu_no->SetLineColor(kNumuBackgroundColor);
894  hNC_no->SetLineColor(kNCBackgroundColor);
895  hNueSig_no->SetLineColor(kNueSignalColor);
896 
897  hNueBeam_sh->SetLineColor(kBeamNueBackgroundColor);
898  hNumu_sh->SetLineColor(kNumuBackgroundColor);
899  hNC_sh->SetLineColor(kNCBackgroundColor);
900  hNueSig_sh->SetLineColor(kNueSignalColor);
901 
902  hNueBeam_sh->SetLineStyle(2);
903  hNumu_sh->SetLineStyle(2);
904  hNC_sh->SetLineStyle(2);
905  hNueSig_sh->SetLineStyle(2);
906 
907  hNueSig_no->SetMaximum(1.3*max(hNueSig_no->GetMaximum(),hNueSig_sh->GetMaximum()));
908  hNueSig_no->SetMaximum(FindLimitY({hNueSig_no,hNueBeam_no,hNC_no,hNumu_no,
909  hNueSig_sh,hNueBeam_sh,hNC_sh,hNumu_sh}));
910  hNueSig_no->SetTitle(plottitle);
911  CenterTitles(hNueSig_no);
912  TGaxis::SetMaxDigits(3);
913 
914  std::vector<TH1*>vnom={hNueSig_no,hNueBeam_no,hNC_no,hNumu_no};
915  std::vector<TH1*>vshi={hNueSig_sh,hNueBeam_sh,hNC_sh,hNumu_sh};
916  auto leg = DefaultNueLegend(true,false,false);
917  TH1* htemp = new TH1F("", "", 1, 0, 1);
918  htemp->SetLineColor(kGray+1);
919  htemp->SetLineStyle(2);
920  leg->AddEntry(htemp,"Shifted","l");
921 
922  PlotMCComponentsComparison(vnom,vshi,leg,pidtag,out_name,printratioplots,false,pidaxis);
923 
924  if(printtable){
925  auto hTotBkg_no = (TH1*) hNueBeam_no->Clone();
926  hTotBkg_no->Add(hNC_no); hTotBkg_no->Add(hNumu_no); hTotBkg_no->Add(hTau_no);
927 
928  auto hTotBkg_sh = (TH1*) hNueBeam_sh->Clone();
929  hTotBkg_sh->Add(hNC_sh); hTotBkg_sh->Add(hNumu_sh); hTotBkg_sh->Add(hTau_sh);
930 
931  std::cout << "% " << std::string(60,'-') << std::endl;
932  std::cout << "% - " << (plottitle.Remove(plottitle.First(";"))).Data() << " - " << pidtag << " - " << std::endl;
933  ComparisonTable ({hNueSig_no, hTotBkg_no, hNueBeam_no,hNC_no,hNumu_no,hTau_no},
934  {hNueSig_sh, hTotBkg_sh, hNueBeam_sh,hNC_sh,hNumu_sh,hTau_sh},
935  {"$\\nu_e$ signal","Total Beam Bkg","Beam $\\nu_{e}$CC","NC", "$\\nu_{\\mu}$CC", "$\\nu_{\\tau}$CC"},
936  (plottitle.Remove(plottitle.First(";")))+pidtag
937  );
938  std::cout << "%" << std::string(60,'-') << std::endl;
939  }
940 
941  }
942 
943  void FormatErrorBand(TH1* hplus, TH1* hminus, bool signal=false, bool fixbins=false){
944 
945  int mccolor = signal ? kViolet -9 : kTotalMCErrorBandColor;
946  hplus ->SetLineColor(mccolor);
947  hminus->SetLineColor(mccolor);
948  hplus ->SetFillColor(mccolor);
949  hminus->SetFillColor(10);
950  hplus ->SetMarkerColor(kRed-7); hminus->SetMarkerColor(kBlue-7);
951 
952  if (fixbins){
953  for(int i=1;i<=hplus->GetNbinsX();i++){
954  double tplus = hplus->GetBinContent(i);
955  double tminus = hminus->GetBinContent(i);
956  hplus->SetBinContent(i,max(tplus, tminus));
957  hminus->SetBinContent(i,min(tplus,tminus));
958  }
959  }
960 
961  }
962 
963  void PlotMCComponentsErrorBand(std::vector<TH1*>vnom,std::vector<TH1*>vshi,
964  TLegend * leg, TString pidtag="",
965  TString out_name="plot_nd",bool ratioplot=false,
966  bool ratioerrors=false, bool pidaxis = 0)
967  {
968  auto hplu_plot = (TH1*) vshi[0]->Clone();
969  auto hmin_plot = (TH1*) vshi[1]->Clone();
970  FormatErrorBand(hplu_plot,hmin_plot,out_name.Contains("_sig"),true);
971  std::vector<TH1*>vshi_plot={hplu_plot,hmin_plot};
972 
973  if(pidaxis) FixLegend(leg,"ExPID");
974  if(!ratioplot){
975  auto c2 = new TCanvas ("c2","c2");
976  vnom[0]->SetMaximum(1.2*vnom[0]->GetMaximum());
977  vnom[0]->Draw("hist");
978  for(auto v:vshi_plot) v->Draw("hist same");
979  for(auto v:vnom) v->Draw("hist same");
980  if (pidaxis){
981  PIDDivisions();
982  PIDBinLabels((std::string)pidtag);
983  PID2DAxis(vnom[0]);
984  }
985  PIDTag(pidtag);
986  leg->Draw("same");
987  c2->SaveAs("plots/"+out_name+"_mccomp.pdf");
988  c2->SaveAs("plots/png_"+out_name+"_mccomp.png");
989  c2->SaveAs("plots/rootfiles/"+out_name+"_mccomp.root");
990  c2->Close();
991  }
992  else{
993  std::vector<TH1*> vratio;
994  for(unsigned int ii=0;ii<vshi.size();ii++){
995  auto htemp = (TH1*) vshi[ii]->Clone(ana::UniqueName().c_str());
996  htemp->Divide(vshi[ii],vnom[0]);
997  htemp->SetMarkerStyle(34);
998  htemp->SetMarkerSize(1.2);
999  htemp->GetYaxis()->SetTitle("Shift / Nom");
1000  vratio.push_back(htemp);
1001  }
1002  vratio[0]->GetYaxis()->SetRangeUser(0.5,1.5);
1003  vnom.insert(vnom.begin()+1,vshi_plot.begin(),vshi_plot.end());
1004  std::vector<TString> topOption (vnom.size(),"hist");
1005  std::vector<TString> bottomOption (vratio.size(),ratioerrors ?"p":"hist p");
1006  auto cratio2= RatioPlot(vnom, topOption,vratio,bottomOption,pidtag,pidaxis);
1007  cratio2->cd();
1008  leg->Draw("same");
1009  PIDTag(pidtag);
1010  cratio2->SaveAs("plots/"+out_name+"_mccomp_ratio.pdf");
1011  cratio2->SaveAs("plots/png_"+out_name+"_mccomp_ratio.png");
1012  cratio2->SaveAs("plots/rootfiles/"+out_name+"_mccomp_ratio.root");
1013  }
1014  }
1015 
1016 ////////////////////////////////////////
1017  struct NueMCComponents{
1018  TH1* sig, *bkg,*beam,*numu,*nc,*tau;
1020  };
1021 
1022  struct NueMCComponents GetNDMCComponents(TDirectory * d_no, double kNDPOT,int linestyle =1){
1024  ret.beam = (*Spectrum::LoadFrom(d_no->GetDirectory("nue_comp")) +
1025  *Spectrum::LoadFrom(d_no->GetDirectory("antinue_comp"))).ToTH1(kNDPOT, kBeamNueBackgroundColor, linestyle);
1026  ret.numu = (*Spectrum::LoadFrom(d_no->GetDirectory("numu_comp")) +
1027  *Spectrum::LoadFrom(d_no->GetDirectory("antinumu_comp"))).ToTH1(kNDPOT, kNumuBackgroundColor, linestyle);
1028  ret.nc = (*Spectrum::LoadFrom(d_no->GetDirectory("nc_tot_comp"))).ToTH1(kNDPOT, kNCBackgroundColor, linestyle);
1029  ret.bkg = (TH1*)ret.beam->Clone();
1030  ret.bkg->Add(ret.numu);ret.bkg->Add(ret.nc);
1031  ret.sig=0; ret.tau=0;
1032  ret.nuetonumu=0; ret.nuetonutau=0;
1033  ret.numutonumu=0; ret.numutonutau=0;
1034  ret.bkg->SetLineColor(kTotalMCColor);
1035  return ret;
1036 
1037  }
1038 ////////////////////////////////////////
1039  void CompareNDDataMC (TDirectory * d_no, TDirectory* d_pl,TDirectory*d_mi,TString plottitle,TString out_name, TString pidtag, bool printtable=true, bool pidaxis = 0){
1040  auto data = Spectrum::LoadFrom(d_no->GetDirectory("data_comp"));
1041  auto kNDPOT = data->POT();
1042  auto hdata = data->ToTH1(kNDPOT);
1043  hdata ->SetMarkerStyle(kFullCircle);
1044 
1045  auto hnom = GetNDMCComponents(d_no,kNDPOT);
1046  auto hplu = GetNDMCComponents(d_pl,kNDPOT);
1047  auto hmin = GetNDMCComponents(d_mi,kNDPOT);
1048 
1049  FormatErrorBand(hplu.bkg,hmin.bkg);
1050 
1051  hnom.bkg->SetMaximum(FindLimitY({hnom.bkg,hdata,hplu.bkg,hmin.bkg}));
1052  hnom.bkg->SetTitle(plottitle);
1053 
1054  CenterTitles(hnom.bkg);
1055  TGaxis::SetMaxDigits(3);
1056  hnom.bkg->GetXaxis()->SetDecimals();
1057 
1058  auto leg1=new TLegend(kLegendDefaultX1,kLegendDefaultY1,kLegendDefaultX2,kLegendDefaultY2);
1059  leg1->SetFillStyle(0);
1060  leg1->AddEntry(hdata, "ND Data","epl");
1061  leg1->AddEntry(hnom.bkg,"Nominal","l");
1062  leg1->AddEntry(hplu.bkg,"Syst.","f");
1063 
1064  PlotNDDataTotalMCComparison(hdata,{hnom.bkg,hplu.bkg,hmin.bkg},
1065  leg1,pidtag,out_name,printratioplots,pidaxis);
1066 
1067 
1068  auto leg2 = DefaultNueLegend(false,true,false);
1069  leg2->AddEntry(hplu.bkg,"Syst.","f");
1070  PlotMCComponentsErrorBand({hnom.bkg,hnom.numu,hnom.nc,hnom.beam},{hplu.bkg,hmin.bkg},
1071  leg2,pidtag,out_name,printratioplots,false,pidaxis);
1072 
1073 
1074  if(printtable){
1075  std::cout << "% " << std::string(60,'-') << std::endl;
1076  std::cout << "% - " << (plottitle.Remove(plottitle.First(";"))).Data() << " - ND - " << pidtag << " - " << std::endl;
1077  ComparisonTable ({hnom.bkg,hnom.beam,hnom.nc,hnom.numu},
1078  {hplu.bkg,hplu.beam,hplu.nc,hplu.numu},
1079  {hmin.bkg,hmin.beam,hmin.nc,hmin.numu},
1080  {"Total MC","Beam $\\nu_{e}$","NC", "$\\nu_{\\mu}$ CC"},
1081  (plottitle.Remove(plottitle.First(";")))+pidtag
1082  );
1083  std::cout << "% " << std::string(60,'-') << std::endl;
1084  }
1085 
1086 
1087  }
1088 
1089  struct NueMCComponents GetFDMCComponents(IPrediction * pred_no,int linestyle=1){
1091  ret.sig = pred_no->PredictComponent(&calc, ana::Flavors::kNuMuToNuE,
1093  ret.beam = pred_no->PredictComponent(&calc, ana::Flavors::kNuEToNuE,
1095  ret.numu = pred_no->PredictComponent(&calc, ana::Flavors::kAllNuMu,
1097  ret.nc = pred_no->PredictComponent(&calc, ana::Flavors::kAll,
1099  ret.tau = pred_no->PredictComponent(&calc, ana::Flavors::kAllNuTau,
1101  ret.bkg = (TH1*) ret.beam->Clone();
1102  ret.bkg->Add(ret.nc); ret.bkg->Add(ret.numu); ret.bkg->Add(ret.tau);
1103  ret.bkg->SetLineColor(kTotalMCColor);
1104 
1105  ret.nuetonumu = pred_no->PredictComponent(&calc, ana::Flavors::kNuEToNuMu,
1107  ret.numutonumu = pred_no->PredictComponent(&calc, ana::Flavors::kNuMuToNuMu,
1109  ret.nuetonutau = pred_no->PredictComponent(&calc, ana::Flavors::kNuEToNuTau,
1111  ret.numutonutau = pred_no->PredictComponent(&calc, ana::Flavors::kNuMuToNuTau,
1113  return ret;
1114  }
1115 
1116  TString Latexify (TString s){
1117  for (int i=0; i<s.Length(); i++){
1118  if (s[i]=='_'){
1119  s.Insert(i,"\\");
1120  i++;
1121  }
1122  }
1123  return s;
1124  }
1125 
1126 
1127  void ComparePredictions(IPrediction * pred_no, IPrediction * pred_pl, IPrediction * pred_mi,TString plottitle, TString out_name, TString pidtag, bool printtable=true, bool pidaxis = false){
1128 
1129  auto hnom = GetFDMCComponents(pred_no);
1130  auto hplu = GetFDMCComponents(pred_pl);
1131  auto hmin = GetFDMCComponents(pred_mi);
1132 
1133  TGaxis::SetMaxDigits(3);
1134 
1135  //FD Bkg Error band plot
1136  FormatErrorBand(hplu.bkg,hmin.bkg,false);
1137  hnom.bkg->SetMaximum(FindLimitY({hnom.bkg,hplu.bkg,hmin.bkg}));
1138  hnom.bkg->SetTitle(plottitle);
1139  CenterTitles(hnom.bkg);
1140  auto leg = DefaultNueLegend(false,true,false);
1141  leg->AddEntry(hplu.bkg,"Syst.","f");
1142 
1143  PlotMCComponentsErrorBand({hnom.bkg,hnom.beam,hnom.nc,hnom.numu},{hplu.bkg,hmin.bkg},
1144  leg,pidtag,out_name+"_bkg",printratioplots,false,pidaxis);
1145 
1146  //FD Sig Error band plot
1147  FormatErrorBand(hplu.sig,hmin.sig,true);
1148  hnom.sig->SetMaximum(FindLimitY({hnom.sig,hplu.sig,hmin.sig}));
1149  hnom.sig->SetTitle(plottitle);
1150  CenterTitles(hnom.sig);
1151  auto leg2 = new TLegend(kLegendDefaultX1,kLegendDefaultY1,kLegendDefaultX2,kLegendDefaultY2);
1152  leg2->AddEntry(hnom.sig,"#nu_{e} Signal","l");
1153  leg2->AddEntry(hplu.sig,"Syst.","f");
1154 
1155  PlotMCComponentsErrorBand({hnom.sig},{hplu.sig,hmin.sig},leg2,pidtag,out_name+"_sig",printratioplots,false,pidaxis);
1156 
1157  //Full comparison table
1158  if(printtable){
1159  std::cout << "% "<< std::string(60,'-') << std::endl;
1160  std::cout << "% - " << (plottitle.Remove(plottitle.First(";"))).Data() << " - " << pidtag << " - " << std::endl;
1161  ComparisonTable ({hnom.sig, hnom.bkg, hnom.beam, hnom.nc, hnom.numu, hnom.tau},
1162  {hplu.sig, hplu.bkg, hplu.beam, hplu.nc, hplu.numu, hplu.tau},
1163  {hmin.sig, hmin.bkg, hmin.beam, hmin.nc, hmin.numu, hmin.tau},
1164  {"$\\nu_e$ signal","Total Beam Bkg","Beam $\\nu_{e}$CC","NC", "$\\nu_{\\mu}$CC", "$\\nu_{\\tau}$CC"},
1165  out_name
1166  );
1167  std::cout << "% "<< std::string(60,'-') << std::endl;
1168  }
1169 
1170  }
1171 
1172  void MakeNueSystematicsFile(IPrediction * pred_no, IPrediction * pred_sh, TString pid, TString sigma, TFile* nue_syst_file){
1173  auto hnom = GetFDMCComponents(pred_no);
1174  auto hshi = GetFDMCComponents(pred_sh);
1175 
1176  std::vector<TString> channels = {"numutonumucc" ,"numutonuecc",
1177  "numutonutaucc","nuetonumucc",
1178  "nuetonuecc" ,"nuetonutaucc" ,
1179  "nc"};
1180  std::vector<TH1*> vshi = {hshi.numutonumu, hshi.sig,
1181  hshi.numutonutau, hshi.nuetonumu,
1182  hshi.beam, hshi.nuetonutau,
1183  hshi.nc};
1184  std::vector<TH1*> vnom = {hnom.numutonumu, hnom.sig,
1185  hnom.numutonutau, hnom.nuetonumu,
1186  hnom.beam, hnom.nuetonutau,
1187  hnom.nc};
1188  if(hshi.sig->GetNbinsX()!=30) {
1189  std::cerr << "\nWARNING\n"
1190  <<"MakeNueSystematicsFile doesn't know what to do with this binning\n"
1191  <<"Nothing will be saved\n";
1192  return;
1193  };
1194 
1195  const ana::Binning* bins = 0;
1196  if(pid == "LID") bins = &ana::kLIDNLBinning;
1197  if(pid == "LEM") bins = &ana::kLEMNLBinning;
1198  if(pid == "CVN") bins = &ana::kCVNNLBinning;
1199  assert(bins);
1200 
1201  Double_t pidEdges[]={bins->Edges()[0],bins->Edges()[1],bins->Edges()[2],bins->Edges()[3]};
1202  for(unsigned int i=0;i<vshi.size();i++){
1203  vshi[i]->Divide(vnom[i]);
1204  TString hname = "fdratio_"+pid+"_"+channels[i]+"_"+sigma;
1205  TString htitle = "FD "+channels[i]+";Reconstructed Energy (GeV);"+pid;
1206 
1207  TH2D* h_epid = new TH2D(hname,htitle,3,pidEdges,10,0,5);
1208  for(int xx=1;xx<=3;xx++){
1209  for(int yy=1;yy<=10;yy++){
1210  double ratio = vshi[i]->GetBinContent(yy+10*(xx-1));
1211  // double shift = ratio>0? ratio -1:0;
1212  if(ratio>0)h_epid->SetBinContent(xx,yy,ratio-1);
1213  }
1214  }
1215  nue_syst_file->cd();
1216  h_epid->Write();
1217  }
1218 
1219  }
1220 
1221  void PrintRawEventCounts(TDirectory * dpred, TString title){
1222 
1223  std::vector<int> ndcounts;
1224  std::vector<int> fdcounts;
1225  std::vector<TString> labels;
1226  //Data
1227  labels.push_back("Data - $\\nu_\\mu $");
1228  ndcounts.push_back(((TH1*)dpred->Get("extrap/MEextrap/Decomp/data_comp/hist"))->GetEntries());
1229  fdcounts.push_back(0);
1230  labels.push_back("Data - $\\nu_e $");
1231  ndcounts.push_back(((TH1*)dpred->Get("extrap/EEextrap/Decomp/data_comp/hist"))->GetEntries());
1232  fdcounts.push_back(0);
1233 
1234  //Signal
1235  labels.push_back("Signal - $\\nu_\\mu \\to \\nu_e $");
1236  labels.push_back("Signal - $\\bar{\\nu}_\\mu \\to \\bar{\\nu}_e $");
1237  for(TString xp:{"ME","MEAnti"}){
1238  ndcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/RecoToTrueND/hist"))->GetEntries());
1239  fdcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/TrueToRecoFD/hist"))->GetEntries());
1240  }
1241 
1242  //Major bkg
1243  labels.push_back("Bkg. - $\\nu_e \\to \\nu_e $");
1244  labels.push_back("Bkg. - $\\nu_\\mu \\to \\nu_\\mu $");
1245  labels.push_back("Bkg. - $NC \\to NC $");
1246  for(TString xp:{"EE","MM","NC"}){
1247  ndcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/RecoND/hist"))->GetEntries());
1248  fdcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/TrueToRecoFD/hist"))->GetEntries());
1249  }
1250 
1251  //Minor bkg
1252  labels.push_back("Bkg. - $\\nu_e \\to \\nu_\\tau $");
1253  labels.push_back("Bkg. - $\\nu_\\mu \\to \\nu_\\tau $");
1254  for(TString xp: {"ET","MT"}){
1255  ndcounts.push_back(0);
1256  fdcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/RecoFD/hist"))->GetEntries());
1257  }
1258  labels.push_back("Bkg. - Other");
1259  double fdcount_minor=0;
1260  for(TString xp: {"EM","EEAnti","EMAnti","ETAnti","MMAnti","MTAnti"}){
1261  fdcount_minor+=(((TH1*)dpred->Get("extrap/"+xp+"extrap/RecoFD/hist"))->GetEntries());
1262  }
1263  ndcounts.push_back(0);
1264  fdcounts.push_back(fdcount_minor);
1265 
1266  //Print
1267  std::cout << "% " << std::string(60,'-') << std::endl;
1268  std::cout << title << "\n";
1269  std::cout << "\\begin{tabular}{ l | r | r }\n"
1270  << std::setw(45) << " " << " & "
1271  << std::setw(10) << "ND" << " & "
1272  << std::setw(10) << "FD"
1273  << " \\\\ \\hline \n";
1274  for(unsigned int i=0;i<labels.size();i++){
1275  std::cout << std::setw(45) << labels[i] << " & "
1276  << std::setw(10)<< ndcounts[i] << " & "
1277  << std::setw(10)<< fdcounts[i] << " \\\\ ";
1278  if(i==1 || i==3 || i==6) std::cout << "\\hline \n";
1279  else std::cout << "\n";
1280  }
1281  std::cout << "\\end{tabular} \n";
1282  std::cout << "% " << std::string(60,'-') << std::endl;
1283  }
1284 
1285 
1286 
1287 }
1288 
1289 
1290 
T max(const caf::Proxy< T > &a, T b)
struct NueMCComponents GetNDMCComponents(TDirectory *d_no, double kNDPOT, int linestyle=1)
double FindLimitY(std::vector< TH1 * > histos, bool doMax=true)
std::string MakeLatexCommandName(const std::string &instring)
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Double_t xx
Definition: macro.C:12
Oscillation analysis framework, runs over CAF files outside of ART.
TString Latexify(TString s)
void PIDBinLabels(std::string pid)
set< int >::iterator it
Float_t y1[n_points_granero]
Definition: compare.C:5
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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:553
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
const Binning kCVNNLBinning
Nonlinear binning for CVN.
Definition: Binning.h:121
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
(&#39;beam &#39;)
Definition: IPrediction.h:15
void ComparisonTable(std::vector< TH1 * > mcnom, std::vector< TH1 * > mcshift, std::vector< TString > labels)
const char * p
Definition: xmltok.h:285
static const double kFDPOT
void PIDBinLabels(bool shortaxis)
struct NueMCComponents GetFDMCComponents(IPrediction *pred_no, int linestyle=1)
OStream cerr
Definition: OStream.cxx:7
const double kLegendDefaultX2
TLegend * leg1
Definition: plot_hist.C:105
TH1F * hplus
Definition: Xsec_final.C:75
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1128
TGraphAsymmErrors * GraphWithPoissonErrors(const TH1 *h, bool noErrorsXaxis, bool drawEmptyBins)
Calculate statistical errors appropriate for small Poisson numbers.
Definition: Plots.cxx:842
const Color_t kTotalMCErrorBandColor
Definition: Style.h:18
TH1F * hsig
Definition: Xdiff_gwt.C:59
cout<< t1-> GetEntries()
Definition: plottest35.C:29
const Color_t kNumuBackgroundColor
Definition: Style.h:31
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 Binning bins
Definition: sa_fd_arrays.h:120
void PIDBinLabelsShortAxis()
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)
const Binning pidBins
Definition: NumuCC_CPiBin.h:5
const XML_Char const XML_Char * data
Definition: expat.h:268
void MakeNueSystematicsFile(IPrediction *pred_no, IPrediction *pred_sh, TString pid, TString sigma, TFile *nue_syst_file)
c2
Definition: demo5.py:33
void Prelim()
osc::OscCalculatorDumb calc
const XML_Char * s
Definition: expat.h:262
Charged-current interactions.
Definition: IPrediction.h:39
std::map< ToFCounter, std::vector< unsigned int > > channels
void RatioPlot(const Spectrum &data, const Spectrum &expected, const Spectrum &fit, double miny, double maxy)
Plot data/expected, compared with fit/expected.
Definition: Plots.cxx:215
void CompareNDDataMC(TDirectory *d_no, TDirectory *d_sh, TString plottitle, TString out_name, TString pidtag, bool printtable=true, bool pidaxis=0)
void FixLegend(TLegend *leg, TString opt="default")
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
const Color_t kBeamNueBackgroundColor
Definition: Style.h:24
void PIDDivisionsShortAxis()
double Ebins[11]
Definition: Xdiff_gwt.C:8
const double kLegendDefaultY1
const double kLegendDefaultX1
void PIDTag(TString pid)
TH1 * RebinToShortAxis(TH1 *h)
const double j
Definition: BetheBloch.cxx:29
double he
Definition: runWimpSim.h:113
TH2D * histo
void PID2DShortAxis(TH1 *hist)
void FormatErrorBand(TH1 *hplus, TH1 *hminus, bool signal=false, bool fixbins=false)
float bin[41]
Definition: plottest35.C:14
TLatex * prelim
Definition: Xsec_final.C:133
const std::vector< double > & Edges() const
Definition: Binning.h:30
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
TString MakeLatexCommandName(TString mystring)
const bool printratioplots
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
TLegend * CustomLegend(std::vector< TH1 * >h, std::vector< bool >isdata, std::vector< TString >title, double x1=kLegendDefaultX1, double y1=kLegendDefaultY1, double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
const Binning bins
Definition: NumuCC_CPiBin.h:8
TH1F * h1
const double kLegendDefaultY2
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
ifstream in
Definition: comparison.C:7
TChain * le
Definition: leana.C:12
const Binning kLIDNLBinning
Nonlinear binning for LID.
Definition: Binning.h:114
TH1F * hminus
Definition: Xsec_final.C:76
const Binning kLEMNLBinning
Nonlinear binning for LEM.
Definition: Binning.h:107
Neutral-current interactions.
Definition: IPrediction.h:40
const double kSecondAnaRunMatchedMCPOTScale
Definition: Exposures.h:96
void PlotDataMC(std::vector< TH1 * >vnom, std::vector< bool >isdata, TLegend *leg, TString pidtag="", TString htag="", TString out_name="plot_FD", bool ratioplot=false, bool ratioerrors=false, bool pidaxis=0)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
c1
Definition: demo5.py:24
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
Double_t ymin
Definition: plot.C:24
void CenterTitles(TH1 *histo)
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:17
void PIDDivisions()
void PID2DAxis(TH1 *axes, bool third=false)
void HTag(TString pid)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir)
Definition: Spectrum.cxx:1055
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
void ComparePredictions(IPrediction *pred_no, IPrediction *pred_sh, TString plottitle, TString out_name, TString pidtag, bool printtable=true, bool pidaxis=false)
TCanvas * ExPIDPlot(std::vector< TH1 * > topHistos, std::vector< TString > topOption)
const Color_t kNCBackgroundColor
Definition: Style.h:22
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
void PID2DAxis(TH1 *hist, bool shortaxis)
void PrintRawEventCounts(TDirectory *dpred, TString title)
TLegend * DefaultNueLegend(bool sig=true, bool tot=true, bool data=false, double x1=kLegendDefaultX1, double y1=kLegendDefaultY1, double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
h
Definition: demo3.py:41
void PIDDivisions(bool shortaxis)