temp_draw_plots_util.h

Go to the documentation of this file.

#pragma once

#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Style.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "OscLib/OscCalcDumb.h"

#include "TStyle.h"
#include "TFile.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TGaxis.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLine.h"
#include "TSystem.h"

#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>

osc::OscCalcDumb calc;
static const double kFDPOT= ana::kSecondAnaRunMatchedMCPOTScale;
const bool printratioplots=true;

// Put a "NOvA Preliminary" tag in the corner
void Prelim()
{
  TLatex* prelim = new TLatex(.9, .95, "NOvA Preliminary");
  prelim->SetTextColor(kBlue);
  prelim->SetNDC();
  prelim->SetTextSize(2/36.);
  prelim->SetTextAlign(32);
  prelim->Draw();
}

void CenterTitles(TH1* histo)
{
  histo->GetXaxis()->CenterTitle();
  histo->GetYaxis()->CenterTitle();
  histo->GetZaxis()->CenterTitle();
}

void PIDTag(TString pid)
{
  TLatex* tag = new TLatex(.83, .92, pid);
  tag->SetTextColor(kAzure+3);
  tag->SetNDC();
  tag->SetTextSize(0.038);
  tag->SetTextFont(42);
  tag->SetTextAlign(12);
  tag->Draw();
}

void HTag(TString pid)
{
  TLatex* tag = new TLatex(.18, .825, pid);
  if ( pid == "IH") tag->SetTextColor(kOrange+2);
  if ( pid == "NH") tag->SetTextColor(kAzure+3);
  tag->SetNDC();
  tag->SetTextSize(0.043);
  tag->SetTextFont(42);
  tag->SetTextAlign(12);
  tag->Draw();
}

// NB, loses all the style info, so call early
TH1* RebinToShortAxis(TH1* h)
{
  TH1* ret = new TH1D(ana::UniqueName().c_str(), h->GetTitle(),
                      18, 0, 18);
  ret->GetXaxis()->SetTitle(h->GetXaxis()->GetTitle());
  ret->GetYaxis()->SetTitle(h->GetYaxis()->GetTitle());

  std::vector<std::pair<int, int>> binMap = {{3, 2},
                                             {4, 3},
                                             {5, 4},
                                             {6, 5}, // end of 1st spectrum
                                             {13, 8},
                                             {14, 9},
                                             {15, 10},
                                             {16, 11}, // end of 2nd
                                             {23, 14},
                                             {24, 15},
                                             {25, 16},
                                             {26, 17}};

  for(auto it: binMap) ret->SetBinContent(it.second, h->GetBinContent(it.first));

  return ret;
}


  TString MakeLatexCommandName(TString mystring){
    //There must be a better way
    std::vector<TString> in = {"0","1","2","3","4","5",
                               "6","7","8","9","_"," ","(",")"}; 
    std::vector<TString> out = {"zero","one","two","three","four","five",
                                "six","seven","eight","nine","","","",""};
    for(unsigned int i=0;i<in.size();i++)
      mystring.ReplaceAll(in[i],out[i]);
    return mystring;
    
  }

void PIDDivisions(){
    // Separators between the different PID ranges
    // TGraph not TLine because TGraph is properly clipped by frame
    TGraph* l1 = new TGraph;
    l1->SetPoint(0, 10, 0);
    l1->SetPoint(1, 10, 50000);
    TGraph* l2 = new TGraph;
    l2->SetPoint(0, 20, 0);
    l2->SetPoint(1, 20, 50000);
    for(TGraph* l: {l1, l2}){
      l->SetLineWidth(1);
      l->SetLineColor(kGray+3);
      l->Draw("l same");
    }
}

void PIDBinLabels( std::string pid){
    // Put the labels near the top of the plot
  gPad->Update();  
  const double y = gPad->GetUymax() * 0.955;

    const ana::Binning* bins = 0;
    if(pid == "LID") bins = &ana::kLIDNLBinning;
    if(pid == "LEM") bins = &ana::kLEMNLBinning;
    if(pid == "CVN") bins = &ana::kCVNNLBinning;
    assert(bins);

    for(int i = 0; i < 3; ++i){
//      TLatex* ltx = new TLatex(3+6*i, y, TString::Format("%g  <  %s  <  %g",
      TLatex* ltx = new TLatex(5+10*i, y, TString::Format("%g  <  %s  <  %g",
                                                          bins->Edges()[i],
                                                          pid.c_str(),
                                                          bins->Edges()[i+1]));
      ltx->SetTextSize(0.032);
      ltx->SetTextColor(kGray+3);
      ltx->SetTextAlign(22);
      ltx->Draw();
    }

}

void PID2DAxis(TH1* axes, bool third = false)
{

  double xmax = axes->GetXaxis()->GetXmax();  
  double xmin = axes->GetXaxis()->GetXmin();  
  double xrange = xmax-xmin;

  double ymin = axes->GetMinimum();
  TGaxis* gax1 = new TGaxis(0, ymin, 8, ymin, 0, 4, 4, "");
  TGaxis* gax2 = new TGaxis(10, ymin, 18, ymin, 0, 4, 4, "");
  TGaxis* gax3 = new TGaxis(20, ymin, 30, ymin, 0, 5, 5, "");
  
  float size = axes->GetXaxis()->GetLabelSize();
  float offs = axes->GetXaxis()->GetLabelOffset();
  int   font = axes->GetXaxis()->GetLabelFont();
  // And then hide its title
  axes->GetXaxis()->SetLabelSize(0);
  // and make all the tick marks the same size
  axes->GetXaxis()->SetNdivisions(3000);
  axes->GetXaxis()->SetTitle("Reconstructed Neutrino Energy (GeV)");
//  axes->GetYaxis()->SetTitle("Events / 6.05 x 10^{20} POT-equiv");
  axes->GetYaxis()->SetTitleOffset(0.8);
  axes->SetTitle("");
  CenterTitles(axes);

  if ( third ){
    gax1->SetLabelSize(size);
    gax1->SetLabelOffset(offs);
    gax1->SetLabelFont(font);
    gax1->Draw();
  }
  else
    for(TGaxis* ax: {gax1, gax2, gax3}){
      // Copy settings from pre-existing axis
      ax->SetLabelSize(size-0.001);
      ax->SetLabelOffset(offs);
      ax->SetLabelFont(font);
      ax->Draw();
    }

  TLine* xline = new TLine();
  xline->SetLineWidth(2);
  xline->DrawLine(xmin, 0, xmax, 0);
  gPad->RedrawAxis();
}
//----------------------------------------------------------------------
void PIDDivisionsShortAxis()
{
  // Separators between the different PID ranges
  // TGraph not TLine because TGraph is properly clipped by frame
  TGraph* l1 = new TGraph;
  l1->SetPoint(0, 6, 0);
  l1->SetPoint(1, 6, 50000);
  TGraph* l2 = new TGraph;
  l2->SetPoint(0, 12, 0);
  l2->SetPoint(1, 12, 50000);
  for(TGraph* l: {l1, l2}){
    //      l->SetLineStyle(3);
    l->SetLineWidth(1);
    l->SetLineColor(kGray+3);
    l->Draw("l same");
  }

  // Hardcoding maxy is obviously unsustainable
  TGaxis* gax1 = new TGaxis(6, 0, 6, 18, 0, 18, 504, "U+-");
  TGaxis* gax2 = new TGaxis(12, 0, 12, 18, 0, 18, 504, "U+-");
}
//----------------------------------------------------------------------
void PIDBinLabelsShortAxis()
{
  const std::string fPID = "CVN"; // HACK

  // Put the labels near the top of the plot
  std::string pid = (std::string)fPID;
  gPad->Update();
  const double y = gPad->GetUymax() * .94;

  const ana::Binning* bins = 0;
  if(fPID == "LID") bins = &ana::kLIDNLBinning;
  if(fPID == "LEM") bins = &ana::kLEMNLBinning;
  if(fPID == "CVN") bins = &ana::kCVNNLBinning;
  assert(bins);

  for(int i = 0; i < 3; ++i){
    TLatex* ltx = new TLatex(3+6*i, y, TString::Format("%g < %s < %g",
                                                       bins->Edges()[i],
                                                       pid.c_str(),
                                                       bins->Edges()[i+1]));
    ltx->SetTextSize(.04);//0.035);
    ltx->SetTextAlign(22);
    ltx->Draw();
  }
}
void PID2DShortAxis(TH1* hist)
{
  double ymin = hist->GetMinimum();
  // Paint our own x-axis fragments
  TGaxis* gax1 = new TGaxis(1, ymin, 5, ymin, 1, 3, 3, "");
  TGaxis* gax2 = new TGaxis(7, ymin, 11, ymin, 1, 3, 3, "");
  TGaxis* gax3 = new TGaxis(13, ymin, 17, ymin, 1, 3, 3, "");

  float tick = hist->GetXaxis()->GetTickLength();
  float size = hist->GetXaxis()->GetLabelSize();
  float offs = hist->GetXaxis()->GetLabelOffset();
  int   font = hist->GetXaxis()->GetLabelFont();
  // And then hide its title
  hist->GetXaxis()->SetLabelSize(0);
  // and make all the tick marks the same size
  hist->GetXaxis()->SetNdivisions(1800);
  hist->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
  hist->GetYaxis()->SetTitle("Events / 0.5 GeV Bin");
  hist->GetYaxis()->SetTitleOffset(0.8);
  hist->GetXaxis()->CenterTitle();
  hist->GetYaxis()->CenterTitle();
  hist->SetTitle("");

  if (false){// fThird ){
    gax1->SetTickSize(tick);
    gax1->SetLabelSize(size);
    gax1->SetLabelOffset(offs);
    gax1->SetLabelFont(font);
    gax1->Draw();
  }
  else{
    for(TGaxis* ax: {gax1, gax2, gax3}){
      // Copy settings from pre-existing axis
      ax->SetTickSize(tick);
      ax->SetLabelSize(size);
      ax->SetLabelOffset(offs);
      ax->SetLabelFont(font);
      ax->Draw();
    }
  }
}

double FindLimitY (std::vector<TH1*> histos, bool doMax = true)
{
  double max = histos[0]->GetMaximum();
  double min = histos[0]->GetMinimum();

  for(unsigned int i=0; i<histos.size(); i++){
    if ( histos[i]->GetMaximum() > max ) max = histos[i]->GetMaximum();
    if ( histos[i]->GetMinimum() < min ) min = histos[i]->GetMinimum();
  }

  if ( doMax ) return max;
  else         return min;

}

namespace ana
{

  void ComparisonTable(std::vector<TH1*> mcnom, std::vector<TH1*> mcshift, std::vector<TString> labels
                       ,TString ltxcommand = ""){
    unsigned int nentries = mcnom.size();
    if(mcshift.size()!=nentries || labels.size()!=nentries){std::cerr << "incompatible sizes\n"; return;}
    
    if(ltxcommand != "") std::cout << "\n\\newcommand{\\" << MakeLatexCommandName(ltxcommand) << "}\n";
    std::cout << "{\\begin{tabular}{ l | r r r } \n"
              << std::setw(20) <<" "      << " & "
              << std::setw(10) << "Nominal"  << " & "
              << std::setw(10) << "Shift" << " & "
              << std::setw(12) << "\\% Diff." 
              <<"\t \\\\ \\hline \n";
    for(unsigned int i=0;i<nentries;i++){
      std::cout.unsetf(std::ios_base::floatfield);
      std::cout << std::setprecision(3) 
                << std::setw(20) << labels[i]              << " & " 
                << std::setw(10) << mcnom[i]  ->Integral() << " & " 
                << std::setw(10) << mcshift[i]->Integral() << " & "  
                << std::setw(10) << std::setprecision(1) << std::fixed 
                << ( mcnom[i]->Integral()>0 ? (((mcshift[i]->Integral()/mcnom[i]->Integral()) -1.)*100.):0.)
                << "\\% \t \\\\ \n";
    }
    std::cout << "\\end{tabular}\n}\n";
  }

  
  void ComparisonTable(std::vector<TH1*> mcnom, std::vector<TH1*> mcplus,  std::vector<TH1*> mcminus,
                       std::vector<TString> labels,TString ltxcommand = ""){
    unsigned int nentries = mcnom.size();
    if(mcplus.size()!=nentries || mcminus.size()!=nentries || labels.size()!=nentries){std::cerr << "incompatible sizes\n"; return;}
    if(ltxcommand != "") std::cout << "\\newcommand{\\" << MakeLatexCommandName(ltxcommand) << "}\n";
    std::cout << "{ \\begin{tabular}{ l | r r r | r r } \n"
              << std::setw(20) <<" "      << " & "
              << std::setw(10) << "Nominal"  << " & "
              << std::setw(10) << "Shift (+)" << " & "
              << std::setw(10) << "Shift (-)" << " & "
              << std::setw(12) << "\\% Diff. (+)" << " & "
              << std::setw(12) << "\\% Diff. (-)" 
              <<"\t \\\\ \\hline \n";
    for(unsigned int i=0;i<nentries;i++){
      std::cout.unsetf(std::ios_base::floatfield);
      std::cout << std::setprecision(3) 
                << std::setw(20) << labels[i]              << " & " 
                << std::setw(10) << mcnom[i]  ->Integral() << " & " 
                << std::setw(10) << mcplus[i] ->Integral() << " & "  
                << std::setw(10) << mcminus[i]->Integral() << " & "  
                << std::setw(10) << std::setprecision(1) << std::fixed 
                << ( mcnom[i]->Integral()>0 ? (((mcplus [i]->Integral()/mcnom[i]->Integral()) -1.)*100.):0.)
                << "\\% & "
                << std::setw(10) 
                << ( mcnom[i]->Integral()>0 ? (((mcminus[i]->Integral()/mcnom[i]->Integral()) -1.)*100.):0.)
                << "\\% \t \\\\ \n";
    }
    std::cout << "\\end{tabular} } \n";
  }

  ////////////////////////////////////////

  const double kLegendDefaultX1 = 0.63;
  const double kLegendDefaultX2 = 0.85;
  const double kLegendDefaultY1 = 0.68;
  const double kLegendDefaultY2 = 0.88;

  void FixLegend(TLegend * leg, TString opt="default"){
    leg->SetFillStyle(0);
    if(opt=="default"){
      leg->SetX1(kLegendDefaultX1);      leg->SetX2(kLegendDefaultX2);
      leg->SetY1(kLegendDefaultY1);      leg->SetY2(kLegendDefaultY2);
    }
    if(opt=="ExPID"){
      leg->SetX1(0.405);      leg->SetX2(0.565);
      leg->SetY1(0.59);      leg->SetY2(0.82);
      leg->SetBorderSize(0);
      leg->SetTextSize(0.037);
    }
    if(opt=="ExPIDPads"){
      leg->SetX1(0.55);      leg->SetX2(0.89);
      leg->SetY1(0.75);      leg->SetY2(0.89);
      leg->SetTextSize(0.045);
    }
  }
 ////////////////////////////////////////
  TLegend * DefaultNueLegend (bool sig=true, bool tot=true,bool data=false,
                              double x1=kLegendDefaultX1, double y1=kLegendDefaultY1,
                              double x2=kLegendDefaultX2, double y2=kLegendDefaultY2)
  {
    TLegend* leg = new TLegend(x1,y1,x2,y2);
    leg->SetFillStyle(0);
    TH1* hdata = new TH1F("", "", 1, 0, 1);
    TH1* hsig = new TH1F("", "", 1, 0, 1);
    TH1* hNC = new TH1F("", "", 1, 0, 1);
    TH1* hCC = new TH1F("", "", 1, 0, 1);
    TH1* hBeam = new TH1F("", "", 1, 0, 1);
    TH1* htot = new TH1F("", "", 1, 0, 1);
    hdata->SetMarkerStyle(kFullCircle);
    hsig->SetLineColor(kNueSignalColor);
    hNC->SetLineColor(kNCBackgroundColor);
    hCC->SetLineColor(kNumuBackgroundColor);
    hBeam->SetLineColor(kBeamNueBackgroundColor);
    htot->SetLineColor(kTotalMCColor);
    if(sig)  leg->AddEntry(hsig,"#nu_{e} Signal","l");
    if(data) leg->AddEntry(hdata,"ND Data","epl");
    if(tot)leg->AddEntry(htot,"Total Bkg","l");
    leg->AddEntry(hBeam,"Beam #nu_{e}","l");
    leg->AddEntry(hNC,"NC","l");
    leg->AddEntry(hCC,"#nu_{#mu} CC","l");
    return leg;
    
  }
  ////////////////////////////////////////
  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)
  {
    TLegend* leg = new TLegend(x1,y1,x2,y2);
    leg->SetFillStyle(0);
    for ( unsigned int v = 0; v<h.size(); ++v){
      if (isdata[v]) TLegendEntry* le = leg->AddEntry(h[v],title[v],"epl");
      else{
        if (h[v]->GetLineColor() == kRed)  TLegendEntry* le = leg->AddEntry(h[v],title[v],"l");
        else TLegendEntry* le = leg->AddEntry(h[v],title[v],"f");
      }       
    }
    return leg;
    
  }
 ////////////////////////////////////////
  TCanvas * ExPIDPlot (std::vector<TH1*>    topHistos,
                       std::vector<TString> topOption)
  {
    TCanvas *cpid = new TCanvas("cpid", "cpid", 410, 850);
    float Ebins = 10;
    int pidBins = 3;

    TPad* p[pidBins];
    p[0] = new TPad("", "", 0, 0, 1, 1);
    p[1] = new TPad("", "", 0, 0, 1, 1);
    p[2] = new TPad("", "", 0, 0, 1, 1);
    p[0]->SetTopMargin(0.1);
    p[0]->SetBottomMargin(.634);
    p[1]->SetTopMargin(.369);
    p[1]->SetBottomMargin(.369);
    p[2]->SetTopMargin(.634);
    p[2]->SetBottomMargin(0.1);
    p[0]->SetFillStyle(0);
    p[1]->SetFillStyle(0);
    p[2]->SetFillStyle(0);
    p[0]->Draw();
    p[1]->Draw();
    p[2]->Draw();

    TH1* he[pidBins];
    TH1* h[topHistos.size()][pidBins];
    for(int k=pidBins-1; k>=0; --k){
      he[k] = (TH1*) topHistos[0]->Clone();
      he[k]->GetYaxis()->SetRangeUser(FindLimitY(topHistos,0) ,1.2*FindLimitY(topHistos,1));
      he[k]->GetXaxis()->SetRangeUser(Ebins*k,Ebins*(k+1));
      if ( k>0 ) he[k]->GetXaxis()->SetLabelSize(0);
      p[pidBins-1-k]->cd();
      he[k]->Draw(topOption[0]);
      he[k]->GetYaxis()->SetTitleSize(26);
      he[k]->GetYaxis()->SetTitleFont(43);
      he[k]->GetYaxis()->SetTitleOffset(1.25);
      he[k]->GetYaxis()->SetLabelFont(43);
      he[k]->GetYaxis()->SetLabelSize(14);
      he[k]->GetXaxis()->SetLabelOffset(-0.008);
      he[k]->GetXaxis()->SetTitleOffset(0.5);
      p[pidBins-1-k]->SetGridy();
    }
    for(unsigned int ii=0;ii<topHistos.size();ii++){
      for(int j=pidBins-1; j>=0; --j){
        h[ii][j] = (TH1*) topHistos[ii]->Clone();
        h[ii][j]->GetXaxis()->SetRangeUser(Ebins*j,Ebins*(j+1));
        p[pidBins-1-j]->cd();
        h[ii][j]->Draw("same "+topOption[ii]);
      }
    }
    for(int h=pidBins-1; h>=0; --h){
      he[h]->SetStats(0);
      p[pidBins-1-h]->cd();
      he[h]->Draw("same "+topOption[0]);
    }
    p[0]->RedrawAxis(); //needed for white histograms e.g. beam
    p[1]->RedrawAxis();
    p[2]->RedrawAxis();
    return cpid;
  }
  ////////////////////////////////////////
  TCanvas * RatioPlot (std::vector<TH1*>    topHistos,
                       std::vector<TString> topOption, 
                       std::vector<TH1*>    bottomHistos,
                       std::vector<TString> bottomOption, 
                       TString pidtag,
                       bool pidaxis=false)
  {
    TCanvas *c = new TCanvas("c", "canvas", 600+(150*pidaxis), 700-(50*pidaxis));
    auto h1 = (TH1*) topHistos[0]->Clone();
    auto h3 = (TH1*) bottomHistos[0]->Clone();

    CenterTitles(h1);
    CenterTitles(h3);

    TPad* p1 = new TPad("", "", 0, 0, 1, 1);
    TPad* p2 = new TPad("", "", 0, 0, 1, 1);
    p1->SetBottomMargin(.3);
    p2->SetTopMargin(.7);
    p1->SetFillStyle(0);
    p2->SetFillStyle(0);
    p1->Draw();
    p2->Draw();

    p1->cd();               
    h1->SetStats(0);          
    h1->Draw(topOption[0]);
    for(unsigned int ii=0;ii<topHistos.size();ii++){
      topHistos[ii]->Draw("same "+topOption[ii]);}
    h1->Draw("same "+topOption[0]);
    p1->RedrawAxis(); //white histograms cover axes 
    h1->GetYaxis()->SetRangeUser(FindLimitY(topHistos,0) ,1.2*FindLimitY(topHistos,1));
    h1->GetYaxis()->SetTitleSize(26);
    h1->GetYaxis()->SetTitleFont(43);
    h1->GetYaxis()->SetTitleOffset(1.2);
    h1->GetYaxis()->SetLabelFont(43);
    h1->GetYaxis()->SetLabelSize(18);
    h1->GetXaxis()->SetLabelSize(0);
    h1->GetXaxis()->SetTitleSize(0);
    if (pidaxis){
      PIDDivisions();
      PIDBinLabels((std::string)pidtag);
    }

    p2->cd(); 
    h3->SetTitle("");
    h3->GetYaxis()->SetRangeUser(max(0.,2-FindLimitY(bottomHistos,1)) ,1.01*FindLimitY(bottomHistos,1));
    h3->GetYaxis()->SetTitleSize(26);
    h3->GetYaxis()->SetTitleFont(43);
    h3->GetYaxis()->SetTitleOffset(1.2);
    h3->GetYaxis()->SetLabelFont(43); 
    h3->GetYaxis()->SetLabelSize(15);
    h3->GetYaxis()->SetDecimals();
    
    h3->GetXaxis()->SetTitle(h1->GetXaxis()->GetTitle());
    h3->GetXaxis()->SetTitleSize(26);
    h3->GetXaxis()->SetTitleFont(43);
     h3->GetXaxis()->SetLabelFont(43); 
    h3->GetXaxis()->SetLabelSize(18);
    h3->GetXaxis()->SetLabelOffset(0.005);
    h3->Draw(bottomOption[0]);
    for(unsigned int ii=0;ii<bottomHistos.size();ii++){
      bottomHistos[ii]->Draw("same "+bottomOption[ii]);}
    h3->Draw("same "+bottomOption[0]);
    
    auto lone = new TLine();
    lone->SetLineStyle(3);
    p2->Update();
    lone->DrawLine(p2->GetUxmin(),1.0,p2->GetUxmax(),1.0);
    lone->DrawLine(p2->GetUxmin(),1.2,p2->GetUxmax(),1.2);
    lone->DrawLine(p2->GetUxmin(),0.8,p2->GetUxmax(),0.8);
    float p2Y = p2->GetUymax()-p2->GetUymin();
    h3->GetYaxis()->SetRangeUser(p2->GetUymin()-(0.01*p2Y),p2->GetUymax()+(0.01*p2Y));

    if (pidaxis){
      PIDDivisions();
      PID2DAxis(h3);
    }
    return c;
  }
  ////////////////////////////////////////
  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)
  {
    gStyle->SetPadLeftMargin(0.15);
    gStyle->SetPadTopMargin(0.08);

    // Fill empty bins with 0 if in range and <0 if out of range
    for ( unsigned int bin = 1; bin<=vnom[0]->GetNbinsX(); ++bin){
      bool filledbin=false;
      for ( unsigned int v = 0; v<vnom.size(); ++v){
        if ( !isdata[v] ) 
          if (vnom[v]->GetBinContent(bin) > 0 ) filledbin=true;
      }
      for ( unsigned int v = 0; v<vnom.size(); ++v){
        if (vnom[v]->GetBinContent(bin) == 0 ){
          if ( !filledbin ) vnom[v]->SetBinContent(bin,-0.01);
          else vnom[v]->SetBinContent(bin,0.001);
        }
        vnom[v]->GetYaxis()->SetTitleOffset(1.05);
      }
    }

    // Draw marker w/error or line when appropriate
    int datahist=-1;
    for ( unsigned int v = 0; v<vnom.size(); ++v) 
      if (vnom[v]->GetLineColor() == kBlack ) datahist=v;

    // Option to make a ratio plot given
    if(!ratioplot){
      auto c2 = new TCanvas ("c2","c2",715,500);
      for ( unsigned int v = 0; v<vnom.size(); ++v){
        TGraphAsymmErrors* gr = GraphWithPoissonErrors(vnom[v], 0, 1);
        if ( isdata[v] ){
          if (v == datahist ) gr->Draw("ep same");
          else vnom[v]->Draw("p same");
        }
        else             vnom[v]->Draw("hist same");
      }
      gPad->RedrawAxis();

      // Adjust for ExPID axis
      if (pidaxis){
        FixLegend(leg,"ExPID");
        NuePlotStyle style(pidtag,0);
        style.PIDDivisions(true);
        style.PID2DAxis(vnom[0],true);
        style.PIDBinLabels(true);
      }

      // pot label
      TLatex* ltx = new TLatex(.175, 0.74, 
                               TString::Format("%g#times10^{20} POT equiv.", 6.05).Data());
      ltx->SetNDC();
      ltx->SetTextAlign(11);
      ltx->SetTextSize(0.0365);
//      ltx->Draw();
      //    HTag(htag);
      leg->Draw("same");

      // Save with & without prelim tag
      gPad->SetFillStyle(0);
      c2->SetFillStyle(0);
      c2->SaveAs("plots/"+out_name+"_datamc.pdf");
      c2->SaveAs("plots/"+out_name+"_datamc.eps");
      c2->SaveAs("plots/"+out_name+"_datamc.C");
      c2->SaveAs("plots/"+out_name+"_datamc.root");
//      Prelim();
      //    c2->SaveAs("plots/"+out_name+"_pre_datamc.pdf");
      //c2->SaveAs("plots/"+out_name+"_pre_datamc.eps");
      //c2->SaveAs("plots/"+out_name+"_pre_datamc.C");
      //c2->SaveAs("plots/"+out_name+"_pre_datamc.root");
    }
    else{
      std::vector<TH1*> vratio;
      for(unsigned int ii=1;ii<vnom.size();ii++){
        auto htemp = (TH1*) vnom[ii]->Clone(ana::UniqueName().c_str());
        htemp->Divide(vnom[datahist],vnom[ii]);
        htemp->SetMarkerColor(vnom[ii]->GetLineColor());
        htemp->SetMarkerStyle(34);
        htemp->SetMarkerSize(1.2);
        htemp->GetYaxis()->SetTitle("Data / MC");
        htemp->GetYaxis()->SetRangeUser(0.5,1.5);
        vratio.push_back(htemp);
      }
      std::vector<TString> topOption  = {" ","hist"};
      std::vector<TString> bottomOption (vratio.size(),ratioerrors ?"p":"hist p");
      auto cratio2= RatioPlot(vnom, topOption,vratio,bottomOption,pidtag,pidaxis);
      cratio2->cd();
      leg->Draw("same");
      //HTag(pidtag);
      cratio2->SaveAs("plots/"+out_name+"_datamc_ratio.pdf");
      cratio2->SaveAs("plots/png_"+out_name+"_datamc_ratio.png");
      cratio2->SaveAs("plots/rootfiles/"+out_name+"_datamc_ratio.root");
    }
  }

  ////////////////////////////////////////
  void PlotNDDataTotalMCComparison(TH1* hdata,std::vector<TH1*> htots,
                                   TLegend * leg,TString pidtag="",
                                   TString out_name="plot_nd",
                                   bool ratioplot=false, bool pidaxis=false)
  {
    if(pidaxis) FixLegend(leg,"ExPID");
    std::vector<TString> h_opt;
    for(auto h:htots) h_opt.push_back("hist");
    h_opt.push_back("p");
    
    if(!ratioplot){    
      auto c1 = new TCanvas ("c1","c1");
      htots[0]->SetMaximum(1.3*htots[0]->GetMaximum());
      for(auto h:htots) h->Draw("hist same");
      htots[0]->Draw("hist same"); //need to redraw for error bands
      hdata->Draw("p same");
      PIDTag(pidtag);
      if (pidaxis){
        PIDDivisions();
        PIDBinLabels((std::string)pidtag);
        PID2DAxis(htots[0]);
      }
      leg->Draw("same");
      c1->SaveAs("plots/"+out_name+"_datamc.pdf");
      c1->SaveAs("plots/png_"+out_name+"_datamc.png");
      c1->SaveAs("plots/rootfiles/"+out_name+"_datamc.root");
    }
    else{
      std::vector<TH1*> ratios;
      std::vector<TString> r_opt;
      for(auto h:htots){
        auto ratio = (TH1*) h->Clone();
        ratio->Divide(hdata,h);
        ratio->GetYaxis()->SetRangeUser(0.5,1.5);
        ratios.push_back(ratio);
        r_opt.push_back("e1");
      } 
      ratios[0]->GetYaxis()->SetTitle("Data / MC");
      ratios[0]->GetYaxis()->SetRangeUser(0.5,1.5);
      htots.push_back(hdata);
      auto cratio1 =RatioPlot(htots,h_opt,ratios,r_opt,pidtag,pidaxis);
      cratio1->cd();
      leg->Draw("same");
      PIDTag(pidtag);
      cratio1->SaveAs("plots/"+out_name+"_datamc_ratio.pdf");
      cratio1->SaveAs("plots/png_"+out_name+"_datamc_ratio.png");
      cratio1->SaveAs("plots/rootfiles/"+out_name+"_datamc_ratio.root");
    }
    if(pidaxis){
      htots.push_back(hdata);
      auto cpads =ExPIDPlot(htots,h_opt);
      cpads->cd();
      FixLegend(leg,"ExPIDPads");
      leg->Draw("same");
      PIDTag(pidtag);
//      cpads->SaveAs("plots/"+out_name+"_datamc_pads.pdf");
      //    cpads->SaveAs("plots/png_"+out_name+"_datamc_pads.png");
      //cpads->SaveAs("plots/rootfiles/"+out_name+"_datamc_pads.root");
    }
  }

  ////////////////////////////////////////
  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)
  {

    if(pidaxis) FixLegend(leg,"ExPID");
    if(!ratioplot){
      auto c2 = new TCanvas ("c2","c2");
      for(auto v:vnom) v->Draw("hist same");
      for(auto v:vshi) v->Draw("hist same");
      if (pidaxis){
        PIDDivisions();
        PIDBinLabels((std::string)pidtag);
        PID2DAxis(vnom[0]);
      }
      PIDTag(pidtag);
      leg->Draw("same");
      c2->SaveAs("plots/"+out_name+"_mccomp.pdf");
      c2->SaveAs("plots/png_"+out_name+"_mccomp.png");
      c2->SaveAs("plots/rootfiles/"+out_name+"_mccomp.root");
      c2->Close();
    }
    else{
      std::vector<TH1*> vratio;
      for(unsigned int ii=0;ii<vshi.size();ii++){
        auto htemp = (TH1*) vnom[ii]->Clone(ana::UniqueName().c_str());
        htemp->Divide(vshi[ii],vnom[ii]);
        htemp->SetMarkerColor(vshi[ii]->GetLineColor());
        htemp->SetMarkerStyle(34);
        htemp->SetMarkerSize(1.2);
        htemp->GetYaxis()->SetTitle("Shift / Nom");
        htemp->GetYaxis()->SetRangeUser(0.5,1.5);
        vratio.push_back(htemp);
      }
      vnom.insert(vnom.end(),vshi.begin(),vshi.end());
      std::vector<TString> topOption    (vnom.size(),"hist");
      std::vector<TString> bottomOption (vratio.size(),ratioerrors ?"p":"hist p");
      auto cratio2= RatioPlot(vnom, topOption,vratio,bottomOption,pidtag,pidaxis);
      cratio2->cd();
      leg->Draw("same");
      PIDTag(pidtag);
      cratio2->SaveAs("plots/"+out_name+"_mccomp_ratio.pdf");
      cratio2->SaveAs("plots/png_"+out_name+"_mccomp_ratio.png");
      cratio2->SaveAs("plots/rootfiles/"+out_name+"_mccomp_ratio.root");
    }
    if(pidaxis){
      vnom.insert(vnom.end(),vshi.begin(),vshi.end());
      std::vector<TString> topOption    (vnom.size(),"hist");
      auto cpads =ExPIDPlot(vnom,topOption);
      cpads->cd();
      FixLegend(leg,"ExPIDPads");
      leg->Draw("same");
      PIDTag(pidtag);
//      cpads->SaveAs("plots/"+out_name+"_mccomp_pads.pdf");
      //    cpads->SaveAs("plots/png_"+out_name+"_mccomp_pads.png");
      //cpads->SaveAs("plots/rootfiles/"+out_name+"_mccomp_pads.root");

    }
  }

  void CompareNDDataMC (TDirectory * d_no, TDirectory* d_sh,TString plottitle,TString out_name, TString pidtag, bool printtable=true, bool pidaxis = 0){
    auto data = Spectrum::LoadFrom(d_no, "data_comp");
    auto kNDPOT = data->POT();
    auto hdata = data->ToTH1(kNDPOT);
    hdata ->SetMarkerStyle(kFullCircle);
    
    auto hnue_no  = (*Spectrum::LoadFrom(d_no, "nue_comp") +
                     *Spectrum::LoadFrom(d_no, "antinue_comp")).ToTH1(kNDPOT, kBeamNueBackgroundColor, 1);
    auto hnumu_no = (*Spectrum::LoadFrom(d_no, "numu_comp") +
                     *Spectrum::LoadFrom(d_no, "antinumu_comp")).ToTH1(kNDPOT, kNumuBackgroundColor, 1);
    auto hnc_no   = (*Spectrum::LoadFrom(d_no, "nc_tot_comp")).ToTH1(kNDPOT, kNCBackgroundColor, 1);
    auto htot_no = (TH1*)hnue_no->Clone();
    htot_no->Add(hnumu_no);htot_no->Add(hnc_no);
    
    auto hnue_sh  = (*Spectrum::LoadFrom(d_sh, "nue_comp") +
                     *Spectrum::LoadFrom(d_sh, "antinue_comp")).ToTH1(kNDPOT, kBeamNueBackgroundColor, 2);
    auto hnumu_sh = (*Spectrum::LoadFrom(d_sh, "numu_comp") +
                     *Spectrum::LoadFrom(d_sh, "antinumu_comp")).ToTH1(kNDPOT, kNumuBackgroundColor, 2);
    auto hnc_sh   = (*Spectrum::LoadFrom(d_sh, "nc_tot_comp")).ToTH1(kNDPOT, kNCBackgroundColor, 2);
    auto htot_sh = (TH1*)hnue_sh->Clone();
    htot_sh->Add(hnumu_sh);htot_sh->Add(hnc_sh);
    
    htot_no->SetLineColor(kTotalMCColor+1);
    htot_sh->SetLineColor(kTotalMCColor+2);
    
    htot_no->SetMaximum(FindLimitY({htot_no,hdata}));
    htot_no->SetTitle(plottitle);
    hnumu_no->SetMaximum(FindLimitY({hnumu_no,hnc_no,hnue_no,hnumu_sh,hnc_sh,hnue_sh}));
    hnumu_no->SetTitle(plottitle);
 
    CenterTitles(htot_no);
    CenterTitles(hnue_no);
    TGaxis::SetMaxDigits(3);
    htot_no->GetXaxis()->SetDecimals();
    
    auto leg1=new TLegend(0.63,0.67,0.85,0.87);
    leg1->SetFillStyle(0);
    leg1->AddEntry(hdata,   "ND Data","epl");
    leg1->AddEntry(htot_no,"Nominal","l");
    leg1->AddEntry(htot_sh,"Shift","l");
    
    PlotNDDataTotalMCComparison(hdata,{htot_no,htot_sh},leg1,pidtag,out_name,printratioplots,pidaxis);
    
    htot_sh->SetLineColor(kTotalMCColor+1);
    std::vector<TH1*>vnom={hnumu_no,hnc_no,hnue_no};
    std::vector<TH1*>vshi={hnumu_sh,hnc_sh,hnue_sh};
    
    auto leg2 = DefaultNueLegend(false,false,false);
    TH1* htemp = new TH1F("", "", 1, 0, 1);
    htemp->SetLineColor(kGray+1);
    htemp->SetLineStyle(2);
    leg2->AddEntry(htemp,"Shifted","l");

    PlotMCComponentsComparison(vnom,vshi,leg2,pidtag,out_name,printratioplots,false,pidaxis);
    
    if(printtable){
      std::cout << "%" <<std::string(60,'-') << std::endl;
      std:: cout << "% - " << (plottitle.Remove(plottitle.First(";"))).Data() << " - ND - " << pidtag << " - " << std::endl;
      ComparisonTable ({htot_no,hnue_no,hnc_no,hnumu_no},
                       {htot_sh,hnue_sh,hnc_sh,hnumu_sh},
                       {"Total MC","Beam $\\nu_{e}$","NC", "$\\nu_{\\mu}$ CC"},
                       (plottitle.Remove(plottitle.First(";")))+pidtag
                       );
      std::cout << "%" <<std::string(60,'-') << std::endl;
    }
    
    
  }
  
  
  void ComparePredictions(IPrediction * pred_no, IPrediction * pred_sh, TString plottitle, TString out_name, TString pidtag, bool printtable=true, bool pidaxis = false){
    
    TH1* hNueSig_no  = pred_no->PredictComponent(&calc, ana::Flavors::kNuMuToNuE,
                                                 ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT);
    TH1* hNueBeam_no = pred_no->PredictComponent(&calc, ana::Flavors::kNuEToNuE,
                                                 ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT);
    TH1* hNumu_no    = pred_no->PredictComponent(&calc, ana::Flavors::kAllNuMu,
                                                 ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT);
    TH1* hNC_no      = pred_no->PredictComponent(&calc, ana::Flavors::kAll,
                                                 ana::Current::kNC, ana::Sign::kBoth).ToTH1(kFDPOT);
    TH1* hTau_no     = pred_no->PredictComponent(&calc, ana::Flavors::kAllNuTau,
                                                 ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT);

    TH1* hNueSig_sh  = pred_sh->PredictComponent(&calc, ana::Flavors::kNuMuToNuE,
                                                 ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT);
    TH1* hNueBeam_sh = pred_sh->PredictComponent(&calc, ana::Flavors::kNuEToNuE,
                                                 ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT);
    TH1* hNumu_sh    = pred_sh->PredictComponent(&calc, ana::Flavors::kAllNuMu,
                                                 ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT);
    TH1* hNC_sh      = pred_sh->PredictComponent(&calc, ana::Flavors::kAll,
                                                 ana::Current::kNC, ana::Sign::kBoth).ToTH1(kFDPOT);
    TH1* hTau_sh     = pred_sh->PredictComponent(&calc, ana::Flavors::kAllNuTau,
                                                 ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT);
    
//    std::cout<<"Making plots"<<std::endl;
    hNueBeam_no->SetLineColor(kBeamNueBackgroundColor);
    hNumu_no->SetLineColor(kNumuBackgroundColor);
    hNC_no->SetLineColor(kNCBackgroundColor);
    hNueSig_no->SetLineColor(kNueSignalColor);
    
    hNueBeam_sh->SetLineColor(kBeamNueBackgroundColor);
    hNumu_sh->SetLineColor(kNumuBackgroundColor);
    hNC_sh->SetLineColor(kNCBackgroundColor);
    hNueSig_sh->SetLineColor(kNueSignalColor);
    
    hNueBeam_sh->SetLineStyle(2);
    hNumu_sh->SetLineStyle(2);
    hNC_sh->SetLineStyle(2);
    hNueSig_sh->SetLineStyle(2);
    
    hNueSig_no->SetMaximum(1.3*max(hNueSig_no->GetMaximum(),hNueSig_sh->GetMaximum()));
    hNueSig_no->SetMaximum(FindLimitY({hNueSig_no,hNueBeam_no,hNC_no,hNumu_no,
                                       hNueSig_sh,hNueBeam_sh,hNC_sh,hNumu_sh}));
    hNueSig_no->SetTitle(plottitle);
    CenterTitles(hNueSig_no);
    TGaxis::SetMaxDigits(3);
    
    std::vector<TH1*>vnom={hNueSig_no,hNueBeam_no,hNC_no,hNumu_no};
    std::vector<TH1*>vshi={hNueSig_sh,hNueBeam_sh,hNC_sh,hNumu_sh};
    auto leg = DefaultNueLegend(true,false,false);
    TH1* htemp = new TH1F("", "", 1, 0, 1);
    htemp->SetLineColor(kGray+1);
    htemp->SetLineStyle(2);
    leg->AddEntry(htemp,"Shifted","l");

    PlotMCComponentsComparison(vnom,vshi,leg,pidtag,out_name,printratioplots,false,pidaxis);
    
    if(printtable){
      auto hTotBkg_no = (TH1*) hNueBeam_no->Clone();
      hTotBkg_no->Add(hNC_no);     hTotBkg_no->Add(hNumu_no);     hTotBkg_no->Add(hTau_no);
      
      auto hTotBkg_sh = (TH1*) hNueBeam_sh->Clone();
      hTotBkg_sh->Add(hNC_sh);     hTotBkg_sh->Add(hNumu_sh);     hTotBkg_sh->Add(hTau_sh);
      
      std::cout << "% " << std::string(60,'-') << std::endl;
      std::cout << "% - " << (plottitle.Remove(plottitle.First(";"))).Data() << " - " << pidtag << " - " << std::endl;
      ComparisonTable ({hNueSig_no, hTotBkg_no, hNueBeam_no,hNC_no,hNumu_no,hTau_no},
                       {hNueSig_sh, hTotBkg_sh, hNueBeam_sh,hNC_sh,hNumu_sh,hTau_sh},
                       {"$\\nu_e$ signal","Total Beam Bkg","Beam $\\nu_{e}$CC","NC", "$\\nu_{\\mu}$CC", "$\\nu_{\\tau}$CC"},
                       (plottitle.Remove(plottitle.First(";")))+pidtag
                       );
      std::cout << "%" << std::string(60,'-') << std::endl;
    }    
    
  }
  
  void FormatErrorBand(TH1* hplus, TH1* hminus, bool signal=false, bool fixbins=false){
    
    int mccolor = signal ? kViolet -9 : kTotalMCErrorBandColor;
    hplus ->SetLineColor(mccolor);
    hminus->SetLineColor(mccolor);                                                                                                                                                        
    hplus ->SetFillColor(mccolor);                                                                                                                                                        
    hminus->SetFillColor(10);                                                                                                                                                                            
    hplus ->SetMarkerColor(kRed-7);                                                                                                                                                                  hminus->SetMarkerColor(kBlue-7);
    
    if (fixbins){
      for(int i=1;i<=hplus->GetNbinsX();i++){
        double tplus = hplus->GetBinContent(i);
        double tminus = hminus->GetBinContent(i);
        hplus->SetBinContent(i,max(tplus, tminus));
        hminus->SetBinContent(i,min(tplus,tminus));
      }
    }
    
  }
  
  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)
  {
    auto hplu_plot = (TH1*) vshi[0]->Clone();
    auto hmin_plot = (TH1*) vshi[1]->Clone();
    FormatErrorBand(hplu_plot,hmin_plot,out_name.Contains("_sig"),true);
    std::vector<TH1*>vshi_plot={hplu_plot,hmin_plot};

    if(pidaxis) FixLegend(leg,"ExPID");
    if(!ratioplot){
      auto c2 = new TCanvas ("c2","c2");
      vnom[0]->SetMaximum(1.2*vnom[0]->GetMaximum());
      vnom[0]->Draw("hist");
      for(auto v:vshi_plot) v->Draw("hist same");
      for(auto v:vnom) v->Draw("hist same");
      if (pidaxis){
        PIDDivisions();
        PIDBinLabels((std::string)pidtag);
        PID2DAxis(vnom[0]);
      }
      PIDTag(pidtag);
      leg->Draw("same");
      c2->SaveAs("plots/"+out_name+"_mccomp.pdf");
      c2->SaveAs("plots/png_"+out_name+"_mccomp.png");
      c2->SaveAs("plots/rootfiles/"+out_name+"_mccomp.root");
      c2->Close();
    }
    else{
      std::vector<TH1*> vratio;
      for(unsigned int ii=0;ii<vshi.size();ii++){
        auto htemp = (TH1*) vshi[ii]->Clone(ana::UniqueName().c_str());
        htemp->Divide(vshi[ii],vnom[0]);
        htemp->SetMarkerStyle(34);
        htemp->SetMarkerSize(1.2);
        htemp->GetYaxis()->SetTitle("Shift / Nom");
        vratio.push_back(htemp);
      }
      vratio[0]->GetYaxis()->SetRangeUser(0.5,1.5);
      vnom.insert(vnom.begin()+1,vshi_plot.begin(),vshi_plot.end());
      std::vector<TString> topOption    (vnom.size(),"hist");
      std::vector<TString> bottomOption (vratio.size(),ratioerrors ?"p":"hist p");
      auto cratio2= RatioPlot(vnom, topOption,vratio,bottomOption,pidtag,pidaxis);
      cratio2->cd();
      leg->Draw("same");
      PIDTag(pidtag);
      cratio2->SaveAs("plots/"+out_name+"_mccomp_ratio.pdf");
      cratio2->SaveAs("plots/png_"+out_name+"_mccomp_ratio.png");
      cratio2->SaveAs("plots/rootfiles/"+out_name+"_mccomp_ratio.root");
    }
  }

////////////////////////////////////////
  struct NueMCComponents{
    TH1* sig, *bkg,*beam,*numu,*nc,*tau;
    TH1* nuetonumu,*numutonumu,*nuetonutau,*numutonutau;
  };
  
  struct NueMCComponents GetNDMCComponents(TDirectory * d_no, double kNDPOT,int linestyle =1){
    NueMCComponents ret;
    ret.beam = (*Spectrum::LoadFrom(d_no, "nue_comp") +
                *Spectrum::LoadFrom(d_no, "antinue_comp")).ToTH1(kNDPOT, kBeamNueBackgroundColor, linestyle);
    ret.numu = (*Spectrum::LoadFrom(d_no, "numu_comp") +
                *Spectrum::LoadFrom(d_no, "antinumu_comp")).ToTH1(kNDPOT, kNumuBackgroundColor, linestyle);
    ret.nc   = (*Spectrum::LoadFrom(d_no, "nc_tot_comp")).ToTH1(kNDPOT, kNCBackgroundColor, linestyle);
    ret.bkg = (TH1*)ret.beam->Clone();
    ret.bkg->Add(ret.numu);ret.bkg->Add(ret.nc);
    ret.sig=0; ret.tau=0;
    ret.nuetonumu=0; ret.nuetonutau=0;
    ret.numutonumu=0; ret.numutonutau=0;
    ret.bkg->SetLineColor(kTotalMCColor);
    return ret;
    
  }
////////////////////////////////////////
  void CompareNDDataMC (TDirectory * d_no, TDirectory* d_pl,TDirectory*d_mi,TString plottitle,TString out_name, TString pidtag, bool printtable=true, bool pidaxis = 0){
    auto data = Spectrum::LoadFrom(d_no, "data_comp");
    auto kNDPOT = data->POT();
    auto hdata = data->ToTH1(kNDPOT);
    hdata ->SetMarkerStyle(kFullCircle);
    
    auto hnom = GetNDMCComponents(d_no,kNDPOT);
    auto hplu = GetNDMCComponents(d_pl,kNDPOT);
    auto hmin = GetNDMCComponents(d_mi,kNDPOT);
    
    FormatErrorBand(hplu.bkg,hmin.bkg);

    hnom.bkg->SetMaximum(FindLimitY({hnom.bkg,hdata,hplu.bkg,hmin.bkg}));
    hnom.bkg->SetTitle(plottitle);
    
    CenterTitles(hnom.bkg);
    TGaxis::SetMaxDigits(3);
    hnom.bkg->GetXaxis()->SetDecimals();
   
    auto leg1=new TLegend(kLegendDefaultX1,kLegendDefaultY1,kLegendDefaultX2,kLegendDefaultY2);
    leg1->SetFillStyle(0);
    leg1->AddEntry(hdata,  "ND Data","epl");
    leg1->AddEntry(hnom.bkg,"Nominal","l");
    leg1->AddEntry(hplu.bkg,"Syst.","f");

    PlotNDDataTotalMCComparison(hdata,{hnom.bkg,hplu.bkg,hmin.bkg},
                                leg1,pidtag,out_name,printratioplots,pidaxis);
    

    auto leg2 = DefaultNueLegend(false,true,false);
    leg2->AddEntry(hplu.bkg,"Syst.","f");
    PlotMCComponentsErrorBand({hnom.bkg,hnom.numu,hnom.nc,hnom.beam},{hplu.bkg,hmin.bkg},
                              leg2,pidtag,out_name,printratioplots,false,pidaxis);

    
    if(printtable){
      std::cout << "% " << std::string(60,'-') << std::endl;
      std::cout << "% - " << (plottitle.Remove(plottitle.First(";"))).Data() << " - ND - " << pidtag << " - " << std::endl;
      ComparisonTable ({hnom.bkg,hnom.beam,hnom.nc,hnom.numu},
                       {hplu.bkg,hplu.beam,hplu.nc,hplu.numu},
                       {hmin.bkg,hmin.beam,hmin.nc,hmin.numu},
                       {"Total MC","Beam $\\nu_{e}$","NC", "$\\nu_{\\mu}$ CC"},
                       (plottitle.Remove(plottitle.First(";")))+pidtag
        );
      std::cout << "% " << std::string(60,'-') << std::endl;
    }
    
    
  }
  
  struct NueMCComponents GetFDMCComponents(IPrediction * pred_no,int linestyle=1){
    NueMCComponents ret;
    ret.sig  = pred_no->PredictComponent(&calc, ana::Flavors::kNuMuToNuE,
                                         ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT,kNueSignalColor,linestyle);
    ret.beam = pred_no->PredictComponent(&calc, ana::Flavors::kNuEToNuE,
                                         ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT,kBeamNueBackgroundColor,linestyle);
    ret.numu = pred_no->PredictComponent(&calc, ana::Flavors::kAllNuMu,
                                         ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT,kNumuBackgroundColor,linestyle);
    ret.nc   = pred_no->PredictComponent(&calc, ana::Flavors::kAll,
                                         ana::Current::kNC, ana::Sign::kBoth).ToTH1(kFDPOT,kNCBackgroundColor,linestyle);
    ret.tau   = pred_no->PredictComponent(&calc, ana::Flavors::kAllNuTau,
                                          ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT,kNumuBackgroundColor,linestyle);
    ret.bkg = (TH1*) ret.beam->Clone();
    ret.bkg->Add(ret.nc);     ret.bkg->Add(ret.numu);     ret.bkg->Add(ret.tau);
    ret.bkg->SetLineColor(kTotalMCColor);

    ret.nuetonumu   = pred_no->PredictComponent(&calc, ana::Flavors::kNuEToNuMu,
                                                ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT,kNumuBackgroundColor,linestyle);
    ret.numutonumu  = pred_no->PredictComponent(&calc, ana::Flavors::kNuMuToNuMu,
                                                ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT,kNumuBackgroundColor,linestyle);
    ret.nuetonutau  = pred_no->PredictComponent(&calc, ana::Flavors::kNuEToNuTau,
                                                ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT,kNumuBackgroundColor,linestyle);
    ret.numutonutau = pred_no->PredictComponent(&calc, ana::Flavors::kNuMuToNuTau,
                                                ana::Current::kCC, ana::Sign::kBoth).ToTH1(kFDPOT,kNumuBackgroundColor,linestyle);
    return ret;
  }  

  TString Latexify (TString s){
    for (int i=0; i<s.Length(); i++){
      if (s[i]=='_'){
        s.Insert(i,"\\");
        i++;
      }
    }
    return s;
  }
  
  
  void ComparePredictions(IPrediction * pred_no, IPrediction * pred_pl, IPrediction * pred_mi,TString plottitle, TString out_name, TString pidtag, bool printtable=true, bool pidaxis = false){
    
    auto hnom = GetFDMCComponents(pred_no);  
    auto hplu = GetFDMCComponents(pred_pl);  
    auto hmin = GetFDMCComponents(pred_mi);  

    TGaxis::SetMaxDigits(3);

    //FD Bkg Error band plot
    FormatErrorBand(hplu.bkg,hmin.bkg,false);
    hnom.bkg->SetMaximum(FindLimitY({hnom.bkg,hplu.bkg,hmin.bkg}));
    hnom.bkg->SetTitle(plottitle);
    CenterTitles(hnom.bkg);
    auto leg = DefaultNueLegend(false,true,false);
    leg->AddEntry(hplu.bkg,"Syst.","f");

    PlotMCComponentsErrorBand({hnom.bkg,hnom.beam,hnom.nc,hnom.numu},{hplu.bkg,hmin.bkg},
                              leg,pidtag,out_name+"_bkg",printratioplots,false,pidaxis);
    
    //FD Sig Error band plot
    FormatErrorBand(hplu.sig,hmin.sig,true);
    hnom.sig->SetMaximum(FindLimitY({hnom.sig,hplu.sig,hmin.sig}));
    hnom.sig->SetTitle(plottitle);
    CenterTitles(hnom.sig);
    auto leg2 = new TLegend(kLegendDefaultX1,kLegendDefaultY1,kLegendDefaultX2,kLegendDefaultY2);
    leg2->AddEntry(hnom.sig,"#nu_{e} Signal","l");
    leg2->AddEntry(hplu.sig,"Syst.","f");

    PlotMCComponentsErrorBand({hnom.sig},{hplu.sig,hmin.sig},leg2,pidtag,out_name+"_sig",printratioplots,false,pidaxis);

    //Full comparison table
    if(printtable){
      std::cout << "% "<< std::string(60,'-') << std::endl;
      std::cout << "% - " << (plottitle.Remove(plottitle.First(";"))).Data() << " - " << pidtag << " - " << std::endl;
      ComparisonTable ({hnom.sig, hnom.bkg, hnom.beam, hnom.nc, hnom.numu, hnom.tau},
                       {hplu.sig, hplu.bkg, hplu.beam, hplu.nc, hplu.numu, hplu.tau},
                       {hmin.sig, hmin.bkg, hmin.beam, hmin.nc, hmin.numu, hmin.tau},
                       {"$\\nu_e$ signal","Total Beam Bkg","Beam $\\nu_{e}$CC","NC", "$\\nu_{\\mu}$CC", "$\\nu_{\\tau}$CC"},
                       out_name
                       );
      std::cout << "% "<< std::string(60,'-') << std::endl;
    }    
    
  }

  void MakeNueSystematicsFile(IPrediction * pred_no, IPrediction * pred_sh, TString pid, TString sigma, TFile* nue_syst_file){
    auto hnom = GetFDMCComponents(pred_no);
    auto hshi = GetFDMCComponents(pred_sh);

    std::vector<TString> channels  = {"numutonumucc" ,"numutonuecc",
                                    "numutonutaucc","nuetonumucc",
                                    "nuetonuecc"   ,"nuetonutaucc" , 
                                    "nc"};
    std::vector<TH1*> vshi  =  {hshi.numutonumu, hshi.sig,
                                hshi.numutonutau, hshi.nuetonumu,
                                hshi.beam, hshi.nuetonutau,
                                hshi.nc};
    std::vector<TH1*> vnom = {hnom.numutonumu, hnom.sig,
                              hnom.numutonutau, hnom.nuetonumu,
                              hnom.beam, hnom.nuetonutau,
                              hnom.nc};
    if(hshi.sig->GetNbinsX()!=30) {
      std::cerr << "\nWARNING\n"
                <<"MakeNueSystematicsFile doesn't know what to do with this binning\n"
                <<"Nothing will be saved\n";
      return;
    };
    
    const ana::Binning* bins = 0;
    if(pid == "LID") bins = &ana::kLIDNLBinning;
    if(pid == "LEM") bins = &ana::kLEMNLBinning;
    if(pid == "CVN") bins = &ana::kCVNNLBinning;
    assert(bins);    
    
    Double_t pidEdges[]={bins->Edges()[0],bins->Edges()[1],bins->Edges()[2],bins->Edges()[3]};
    for(unsigned int i=0;i<vshi.size();i++){
      vshi[i]->Divide(vnom[i]);
      TString hname = "fdratio_"+pid+"_"+channels[i]+"_"+sigma;
      TString htitle = "FD "+channels[i]+";Reconstructed Energy (GeV);"+pid;

      TH2D* h_epid = new TH2D(hname,htitle,3,pidEdges,10,0,5);
      for(int xx=1;xx<=3;xx++){
        for(int yy=1;yy<=10;yy++){
          double ratio = vshi[i]->GetBinContent(yy+10*(xx-1));
          //      double shift = ratio>0? ratio -1:0;
          if(ratio>0)h_epid->SetBinContent(xx,yy,ratio-1);
        }
      }
      nue_syst_file->cd();
      h_epid->Write();
    }
        
  }

  void PrintRawEventCounts(TDirectory * dpred, TString title){
    
    std::vector<int> ndcounts;
    std::vector<int> fdcounts;
    std::vector<TString> labels;
    //Data
    labels.push_back("Data - $\\nu_\\mu $");
    ndcounts.push_back(((TH1*)dpred->Get("extrap/MEextrap/Decomp/data_comp/hist"))->GetEntries());
    fdcounts.push_back(0);
    labels.push_back("Data - $\\nu_e $");
    ndcounts.push_back(((TH1*)dpred->Get("extrap/EEextrap/Decomp/data_comp/hist"))->GetEntries());
    fdcounts.push_back(0);
   
    //Signal 
    labels.push_back("Signal - $\\nu_\\mu  \\to \\nu_e $");
    labels.push_back("Signal - $\\bar{\\nu}_\\mu  \\to \\bar{\\nu}_e $");
    for(TString xp:{"ME","MEAnti"}){
      ndcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/RecoToTrueND/hist"))->GetEntries());
      fdcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/TrueToRecoFD/hist"))->GetEntries());
    }

    //Major bkg 
    labels.push_back("Bkg. - $\\nu_e  \\to \\nu_e $");
    labels.push_back("Bkg. - $\\nu_\\mu  \\to \\nu_\\mu $");
    labels.push_back("Bkg. - $NC \\to NC $");
    for(TString xp:{"EE","MM","NC"}){
      ndcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/RecoND/hist"))->GetEntries());
      fdcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/TrueToRecoFD/hist"))->GetEntries());
    }

    //Minor bkg
    labels.push_back("Bkg. - $\\nu_e  \\to \\nu_\\tau $");
    labels.push_back("Bkg. - $\\nu_\\mu  \\to \\nu_\\tau $");
    for(TString xp: {"ET","MT"}){
      ndcounts.push_back(0); 
      fdcounts.push_back(((TH1*)dpred->Get("extrap/"+xp+"extrap/RecoFD/hist"))->GetEntries());
    }
    labels.push_back("Bkg. - Other");
    double fdcount_minor=0;
    for(TString xp: {"EM","EEAnti","EMAnti","ETAnti","MMAnti","MTAnti"}){
      fdcount_minor+=(((TH1*)dpred->Get("extrap/"+xp+"extrap/RecoFD/hist"))->GetEntries());
    }
    ndcounts.push_back(0);
    fdcounts.push_back(fdcount_minor);

    //Print
    std::cout << "% " << std::string(60,'-') << std::endl;
    std::cout << title << "\n";
    std::cout << "\\begin{tabular}{ l | r | r }\n"
              << std::setw(45) << " " << " & "
              << std::setw(10) << "ND" << " & "
              << std::setw(10) << "FD" 
              << " \\\\ \\hline \n";
    for(unsigned int i=0;i<labels.size();i++){
      std::cout << std::setw(45) << labels[i] << " & "
                << std::setw(10)<< ndcounts[i] << " & "
                << std::setw(10)<< fdcounts[i] << " \\\\ ";
      if(i==1 || i==3 || i==6) std::cout << "\\hline \n";
      else std::cout << "\n";
    }
    std::cout << "\\end{tabular} \n";
    std::cout << "% " << std::string(60,'-') << std::endl;
  }



}