draw_plots_util.h

Go to the documentation of this file.

#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Plots.h"
#include "CAFAna/Analysis/Calcs.h"
#include "3FlavorAna/Plotting/NuePlotStyle.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Core/HistAxis.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "CAFAna/Prediction/PredictionInterp.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "3FlavorAna/Prediction/PredictionExtendToPeripheral.h"
#include "3FlavorAna/Prediction/PredictionSystJoint2018.h"
#include "CAFAna/Vars/FitVars.h"
#include "CAFAna/Decomp/IDecomp.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Extrap/IExtrap.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "CAFAna/Extrap/ModularExtrapComponent.h"

#include "OscLib/OscCalcDumb.h"
#include "OscLib/IOscCalc.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 "TPaveText.h"
#include "TLegend.h"
#include "TLine.h"
#include "TSystem.h"
#include "TError.h"

#include <iostream>
#include <fstream>
#include <iomanip>
#include <cxxabi.h>

#include "tex_tables_util.h"

static bool isFHC=true; // Global variable for FHC/RHC

//static const double kFDPOT= ana::kAna2017POT;
const bool printratioplots=true;

enum AxisType {kNue1bin, kNue3bin, kNumuOpt, kNumuOther};

bool usedumb = false;
bool mergePeripheral=false;

osc::IOscCalc * GetCalculator (bool usingdumb){
  double dCP_Pi = 1.48;
  double ssth23 = 0.558;
  double dmsq32 = 2.44e-3;

  osc::IOscCalc * calc;
  if(usingdumb){
    osc::OscCalcDumb dumb;
    calc = dumb.Copy();
  }
  else{
    auto calc2 = DefaultOscCalc();
    kFitDeltaInPiUnits.SetValue(calc2, dCP_Pi);
    kFitSinSqTheta23.SetValue(calc2, ssth23);
    kFitDmSq32.SetValue(calc2, dmsq32);
    calc = calc2->Copy();
  }
  return calc;
}


//////////////////////////////////////////////////////////////////////
// Histogram
//////////////////////////////////////////////////////////////////////

TH1* PrettyRatio(TH1* shi, TH1* nom, int color,
                 double alpha=1., TString titley="Shift / Nom")
{
  auto htemp = reinterpret_cast<TH1*>(nom->Clone(ana::UniqueName().c_str()));
  htemp->Divide(shi,nom);
  htemp->SetMarkerColorAlpha(color,alpha);
  htemp->SetLineColorAlpha(color,alpha);
  htemp->SetMarkerStyle(34);
  htemp->SetMarkerSize(1.4);
  htemp->GetYaxis()->SetTitle(titley);
  return htemp;
}

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

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;
}

//////////////////////////////////////////////////////////////////////
void FormatErrorBand(TH1* hplus, TH1* hminus,
                     bool signal=false, bool fixbins=false)
{
  int mccolor = signal ? kViolet -9 : kTotalMCErrorBandColor;
  std::cout << "hplus: " << hplus->Integral() << std::endl;
  std::cout << "hminus: " << hminus->Integral() << std::endl;
  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));
    }
  }
}

//////////////////////////////////////////////////////////////////////
TH1* HistNumuOpt (TH1* orig){
  auto ret = (TH1*) orig->Clone(UniqueName().c_str());
  ret->Scale(0.1,"width");
  ret->GetYaxis()->SetTitle("Events / 0.1 GeV");
  ret->SetBit(kCanDelete);
  ret->SetDirectory(gROOT);
  return ret;
}

//////////////////////////////////////////////////////////////////////
// Adornments
//////////////////////////////////////////////////////////////////////

void PrettyTag(TString pid, int color, double xndc, double yndc){
  TLatex* tag = new TLatex(xndc, yndc, pid);
  tag->SetTextColor(color);
  tag->SetNDC();
  tag->SetTextSize(0.04);
  tag->SetTextFont(42);
  tag->SetTextAlign(22);
  tag->Draw();
  tag->SetBit(kCanDelete);
}

void PIDTag(TString pid) { PrettyTag (pid, kAzure+3, .83, .92); }
void HTag  (TString hie) { PrettyTag (hie, (hie=="NH")?kAzure+3:kOrange+2, .18, .825);}

//////////////////////////////////////////////////////////////////////
// Legends
//////////////////////////////////////////////////////////////////////

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

TLegend * CustomLegend(
  std::vector<TH1*>h, std::vector<TString> titles, std::vector<TString> option,
  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)
    leg->AddEntry(h[v], titles[v], option[v]);
  leg->SetBit(kCanDelete);
  return leg;
}

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

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.355); leg->SetX2(0.515);
    leg->SetY1(0.65); leg->SetY2(0.82);
    leg->SetBorderSize(0);
    leg->SetTextSize(0.03);
  }
  if(opt=="ExPIDNearDet")
  {
    leg->SetX1(0.695); leg->SetX2(0.855);
    leg->SetY1(0.65); leg->SetY2(0.82);
    leg->SetBorderSize(0);
    leg->SetTextSize(0.03);
  }
  if(opt=="ExPIDPads")
  {
    leg->SetX1(0.55); leg->SetX2(0.89);
    leg->SetY1(0.75); leg->SetY2(0.89);
    leg->SetTextSize(0.045);
  }
}
//////////////////////////////////////////////////////////////////////

TGraph * DummyGraph(int color=kBlack, int style =1, int linewidth=2, int markerstyle = kFullCircle)
{
  auto gr = new TGraph();
  gr->SetLineColor(color);
  gr->SetLineStyle(style);
  gr->SetLineWidth(linewidth);
  gr->SetMarkerStyle(markerstyle);
  gr->SetMarkerColor(color);
  return gr;
}
//////////////////////////////////////////////////////////////////////

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);
  if(sig)  leg->AddEntry( DummyGraph(kNueSignalColor), "#nu_{e} Signal","l");
  if(data) leg->AddEntry( DummyGraph(),                "ND Data","epl");
  if(tot)  leg->AddEntry( DummyGraph(kTotalMCColor),   "Total Bkg","l");
  leg->AddEntry( DummyGraph(kBeamNueBackgroundColor),  "Beam #nu_{e}","l");
  leg->AddEntry( DummyGraph(kNCBackgroundColor),       "NC","l");
  leg->AddEntry( DummyGraph(kNumuBackgroundColor),     "#nu_{#mu} CC","l");
  leg->SetBit(kCanDelete);
  return leg;

}
//////////////////////////////////////////////////////////////////////

TLegend * DefaultNumuLegend(
  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);
  if(sig)
    leg->AddEntry(DummyGraph(kNueSignalColor), "#nu_{#mu} Signal", "l");
  if(data)
    leg->AddEntry(DummyGraph(), "ND Data", "epl");
  if(tot)
    leg->AddEntry(DummyGraph(kTotalMCColor), "Total Bkg", "l");
  if(isNumuAna)
    leg->AddEntry(DummyGraph(kNumuBackgroundColor), "#nu_{#mu} app", "l");
  else
    leg->AddEntry(DummyGraph(kNumuBackgroundColor), "#nu_{#mu} CC", "l");
  leg->AddEntry(DummyGraph(kNCBackgroundColor), "NC","l");
  leg->AddEntry(DummyGraph(kBeamNueBackgroundColor), "#nu_{e} CC","l");
  leg->AddEntry(DummyGraph(kNuTauBackgroundColor), "#nu_{#tau} CC","l");
  leg->SetBit(kCanDelete);
  return leg;
}
//////////////////////////////////////////////////////////////////////

void PrintCanvasAll (TCanvas * c, const TString outname, const TString outdir = "plots/",
         const std::vector <TString> exts = {"pdf","png","root"} )
{
  for(auto const & ext:exts){
    gSystem -> MakeDirectory( outdir );
    gSystem -> MakeDirectory( outdir+"rootfiles/" );
    if(ext=="png")       c->SaveAs(outdir+"png_"+outname + "." + ext);
    else if(ext=="root") c->SaveAs(outdir+"rootfiles/"+outname + "." + ext);
    else                 c->SaveAs(outdir+outname + "." + ext);
  }
  delete c;
}

//////////////////////////////////////////////////////////////////////
// Nue 3 bins functions
//////////////////////////////////////////////////////////////////////

// NB, loses all the style info, so call early
/*TH1* RebinToShortAxis(TH1* h)
{
  TH1* ret = new TH1D(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;
}
*/
//////////////////////////////////////////////////////////////////////
void NuePID3bin(
  TH1* h, bool coreOnlyND=false, bool binLabels=true,
  const double textSize=0.04, const int textColor=kGray+3,
  const int lineColor=kBlack, const int lineStyle=3, const double yNDC=0.96)
{
  if (coreOnlyND)
    Info("NuePID3bin", "specified core only/ND.");

  Nue2018ThreeBinDivisions(coreOnlyND, lineColor, lineStyle);
  if (binLabels)
  {
    Nue2018ThreeBinLabels(yNDC, textSize, textColor, coreOnlyND);
  }
  else
  {
    if (textSize != 0.04)
      Warning("NuePID3bin",
              "binLabels == false --> specified text size has no effect");
    if (textColor != kGray+3)
      Warning("NuePID3bin",
              "binLabels == false --> specified text color has no effect");
    if (yNDC != 0.96)
      Warning("NuePID3bin",
              "binLabels == false --> specified yNDC has no effect");
  }
}
//////////////////////////////////////////////////////////////////////
// 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(24);
//     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;
// }

//////////////////////////////////////////////////////////////////////
// Pretty Canvas
//////////////////////////////////////////////////////////////////////
TCanvas * RatioPlot(
  std::vector<TH1*> topHistos, std::vector<TString> topOption,
  std::vector<TH1*> bottomHistos, std::vector<TString> bottomOption,
  TString pidtag, AxisType pidaxis=kNue1bin, const bool nearDet=false)
{
  gStyle->SetTitleStyle(0);
  TCanvas *c = new TCanvas("c", "canvas", 600, 700);
  // if(pidaxis == kNue3bin) c->SetCanvasSize(750,650);
  //define top and bottom pads
  TPad* p1 = new TPad("", "", 0, 0, 1, 1);
  TPad* p2 = new TPad("", "", 0, 0, 1, 1);
  p1->SetBottomMargin(.3);
  p2->SetTopMargin(.7);

  for(auto p:{p1,p2}){
    p->SetFillStyle(0);
    p->Draw();
  }

  // format histograms
  auto h1 = (TH1*) topHistos[0]->Clone();
  auto h3 = (TH1*) bottomHistos[0]->Clone();

  for(auto & h:{h1,h3}){
    h->SetStats(0);

    h->GetYaxis()->SetTitleSize(26);
    h->GetYaxis()->SetTitleFont(43);
    h->GetYaxis()->SetTitleOffset(1.3);
    h->GetYaxis()->SetLabelFont(43);

    h->GetXaxis()->SetTitleOffset(1.2);
    h->GetXaxis()->SetTitleSize(28);
    h->GetXaxis()->SetTitleFont(43);
    h->GetXaxis()->SetLabelFont(43);
  }

  if (pidaxis == kNue3bin)
    h1->GetYaxis()->SetRangeUser(
      FindLimitY(topHistos,0), 1.5*FindLimitY(topHistos, 1));
  else
    h1->GetYaxis()->SetRangeUser(
      FindLimitY(topHistos,0), 1.2*FindLimitY(topHistos, 1));
  h1->GetYaxis()->SetLabelSize(18);

  h1->GetXaxis()->SetLabelSize(0);
  h1->GetXaxis()->SetTitleSize(0);

  h3->SetTitle("");
  h3->GetYaxis()->SetRangeUser(max(0.,2-FindLimitY(bottomHistos,1)),
             1.01*FindLimitY(bottomHistos,1));
  h3->GetYaxis()->SetLabelSize(15);
  h3->GetYaxis()->SetDecimals();

  h3->GetXaxis()->SetTitle(h1->GetXaxis()->GetTitle());
  h3->GetXaxis()->SetLabelSize(18);
  h3->GetXaxis()->SetLabelOffset(0.005);

  p1->cd();

  h1->Draw(topOption[0]);

  for(unsigned int ii = 1; ii < topHistos.size(); ++ii) {
    topHistos[ii]->Draw("same "+topOption[ii]);}
  h1->Draw("same "+topOption[0]);

  p1->RedrawAxis(); //white histograms cover axes

  if(pidaxis==kNue3bin)
    NuePID3bin(h1, nearDet);

  p2->cd();
  h3->Draw(bottomOption[0]);
  for(unsigned int ii = 1; ii < bottomHistos.size(); ++ii){
    bottomHistos[ii]->Draw("same " + bottomOption[ii]);}
  h3->Draw("same " + bottomOption[0]);

  auto lone = new TLine();
  lone->SetLineStyle(3);
  lone->SetLineColor(kGray+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));
  p2->RedrawAxis();
  if(pidaxis==kNue3bin)
    NuePID3bin(h3, nearDet, false); // hist, coreOnlyND, binLabels

  return c;
}
//////////////////////////////////////////////////////////////////////
// Pretty plot styles: Data MC, MC components, error bands
//////////////////////////////////////////////////////////////////////

/*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, AxisType pidaxis = kNue1bin)
{
  gStyle->SetPadLeftMargin(0.15);
  gStyle->SetPadTopMargin(0.08);

  // Fill empty bins with 0 if in range and <0 if out of range
  for ( 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 ( uint 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==kNue3bin)
    {
      FixLegend(leg, "ExPID");
      Nue2018ThreeBinDivisions();
      Nue2018ThreeBinLabels();
    }

    // pot label
    PrettyTag(TString::Format("%g#times10^{20} POT equiv.", 8.85).Data(),
        kBlack, .175, .74);
    HTag(htag);
    leg->Draw("same");

    gPad->SetFillStyle(0);
    c2->SetFillStyle(0);
    PrintCanvasAll (c2,out_name+"_datamc","plots/",{"pdf","eps","png","C","root"});
    PrintLatexFigure(out_name+"_datamc");
    //removed: prelim tag
  }
  else{
    std::vector<TH1*> vratio;
    for(unsigned int ii=1;ii<vnom.size();ii++){
      auto htemp = PrettyRatio (vnom[datahist], vnom[ii], vnom[ii]->GetLineColor(), 1.,
        "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);
    PrintCanvasAll(cratio2,out_name+"_datamc_ratio");
    PrintLatexFigure("/png_"+out_name+"_datamc_ratio");
  }
}
*/
//////////////////////////////////////////////////////////////////////

void PlotNDDataTotalMCComparison(
  TH1* hdata, std::vector<TH1*> htots, TLegend * leg,
  TString pidtag="", TString out_name="plot_nd",
  bool ratioplot=false, AxisType pidaxis = kNue1bin)
{
  Info("PlotNDDataTotalMCComparison",
       TString::Format("pidaxis = %d", pidaxis));

  if(pidaxis==kNue3bin)
    FixLegend(leg, "ExPIDNearDet");
  std::vector<TString> h_opt;
  for(auto const & 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 const & 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==kNue3bin)
      NuePID3bin(htots[0], true); // hist, coreOnlyND, binLabels
    leg->Draw("same");
    PrintCanvasAll(c1, out_name+"_datamc");
    PrintLatexFigure("/png_"+out_name+"_datamc");
  }
  else
  {
    std::vector<TH1*> ratios;
    std::vector<TString> r_opt;
    auto idx = 0;
    for(auto const & h:htots)
    {
      auto ratio = PrettyRatio(hdata, h, h->GetLineColor(),
                               idx>0 ? 0.4:1., "Data / MC");
      ratio->GetYaxis()->SetRangeUser(0.5,1.5);
      ratios.push_back(ratio);
      r_opt.push_back("e");
      idx++;
    }
    htots.push_back(hdata);
    auto cratio1 = RatioPlot(htots, h_opt, ratios,
                             r_opt, pidtag, pidaxis, true);
    cratio1->cd();
    leg->Draw("same");
    PIDTag(pidtag);

    PrintCanvasAll(cratio1,out_name+"_datamc_ratio");
    PrintLatexFigure("/png_"+out_name+"_datamc_ratio");
  }
}
//////////////////////////////////////////////////////////////////////
void PlotMCComponentsComparison(
  std::vector<TH1*>vnom, std::vector<TH1*>vshi,
  TLegend * leg, TString pidtag="",
  TString out_name="plot_nd"/* Why ND? */, bool ratioplot=false,
  bool ratioerrors=false, AxisType pidaxis=kNue1bin)
{
  // Is it NearDet?
  Info("PlotMCComponentsErrorBand", "out_name = " + out_name);
  const auto isNearDet = out_name.Contains("ND", TString::kIgnoreCase);
  if (isNearDet)
    Info("PlotMCComponentsErrorBand", "isNearDet = true");
  else
    Info("PlotMCComponentsErrorBand", "isNearDet = false");

  if(pidaxis==kNue3bin)
  {
    if (isNearDet)
      FixLegend(leg, "ExPIDNearDet");
    else
      FixLegend(leg, "ExPID");
  }


  if(!ratioplot)
  {
    auto c2 = new TCanvas(UniqueName().c_str(), "c2");
    for(auto const & v:vnom)
      v->Draw("hist same");
    for(auto const & v:vshi)
      v->Draw("hist same");
    if(pidaxis==kNue3bin)
      NuePID3bin(vnom[0], isNearDet);
    PIDTag(pidtag);
    leg->Draw("same");
    PrintCanvasAll(c2, out_name+"_mccomp");
    PrintLatexFigure("/png_"+out_name+"_mccomp");
    c2->Close();
  }
  else
  {
    std::vector<TH1*> vratio;
    for(unsigned int ii=0;  ii<vshi.size();ii++)
    {
      auto htemp = PrettyRatio(vshi[ii], vnom[ii],
                               vnom[ii]->GetLineColor(), ii>0 ? 0.4:1.);
      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, isNearDet);
    cratio2->cd();
    leg->Draw("same");
    PIDTag(pidtag);
    PrintCanvasAll(cratio2, out_name+"_mccomp_ratio");
    PrintLatexFigure("/png_"+out_name+"_mccomp_ratio");
  }
}

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

void PlotMCComponentsErrorBand(
  std::vector<TH1*>vnom, std::vector<TH1*>vshi,
  TLegend * leg, TString pidtag="", TString out_name="plot_nd",
  bool ratioplot=false, bool ratioerrors=false, AxisType pidaxis=kNue1bin)
{
  // Is it NearDet?
  Info("PlotMCComponentsErrorBand", "out_name = " + out_name);
  const auto isNearDet = out_name.Contains("ND", TString::kIgnoreCase);
  if (isNearDet)
    Info("PlotMCComponentsErrorBand", "isNearDet = true");
  else
    Info("PlotMCComponentsErrorBand", "isNearDet = false");

  auto hplu_plot = reinterpret_cast<TH1*>(vshi[0]->Clone());
  auto hmin_plot = reinterpret_cast<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==kNue3bin)
  {
    if (isNearDet)
      FixLegend(leg, "ExPIDNearDet");
    else
      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==kNue3bin)
      NuePID3bin(vnom[0], isNearDet);
    PIDTag(pidtag);
    leg->Draw("same");
    PrintCanvasAll(c2, out_name+"_mccomp");
    PrintLatexFigure("/png_"+out_name+"_mccomp");
    c2->Close();
  }
  else
  {
    std::vector<TH1*> vratio;
    for(unsigned int ii=0; ii<vshi.size(); ii++)
    {
      auto htemp = PrettyRatio(vshi[ii], vnom[0], vshi[ii]->GetMarkerColor());
      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, isNearDet);
    cratio2->cd();
    leg->Draw("same");
    PIDTag(pidtag);
    PrintCanvasAll(cratio2, out_name+"_mccomp_ratio");
    PrintLatexFigure("/png_"+out_name+"_mccomp_ratio");
  }
}

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



//////////////////////////////////////////////////////////////////////
// Define container for histograms + format
//////////////////////////////////////////////////////////////////////

namespace ana
{

  struct DataMCComponents
  {
    TH1* data;
    TH1* mcTotal;
    TH1* sig, *wrongSign, *bkg, *beam, *numu, *nc, *tau;
    TH1* nuetonumu,*numutonumu,*nuetonutau,*numutonutau;
  };

  void DefaultFormatNue (DataMCComponents & comp, int linestyle = 1){
    std::vector <std::pair <TH1*, int>> pairs = {
      {comp.sig, kNueSignalColor},
      {comp.wrongSign, kNueSignalWSColor},
      {comp.beam, kBeamNueBackgroundColor},
      {comp.nc, kNCBackgroundColor},
      {comp.numu, kNumuBackgroundColor},
      {comp.bkg, kTotalMCColor},
      {comp.mcTotal, kTotalMCColor},
      {comp.nuetonumu, kNumuBackgroundColor},
      {comp.numutonumu, kNumuBackgroundColor},
      {comp.nuetonutau, kGray},
      {comp.numutonutau, kGray}
    };
    for (auto const & mypair:pairs)
    {
      if(mypair.first == 0)
      {
        Info("DefaultFormatNue", "Skipping formatting for component...");
        continue;
      }
      mypair.first->SetLineColor(mypair.second);
      mypair.first->SetLineStyle(linestyle);
      ana::CenterTitles(mypair.first);
    }
    if(comp.data) comp.data->SetMarkerStyle(kFullCircle);
    return;
  }
//////////////////////////////////////////////////////////////////////

  const DataMCComponents GetNDCompsFromDecomp(const IDecomp* decomp)
  {
    // Get components from decomp
    Spectrum nueComp = decomp->MC_NueComponent();
    Info("GetNDComponents", "Loaded NueComponent...");
    Spectrum antinueComp = decomp->MC_AntiNueComponent();
    Info("GetNDComponents", "Loaded AnitNueComponent...");
    Spectrum numuComp = decomp->MC_NumuComponent();
    Info("GetNDComponents", "Loaded NumuComponent...");
    Spectrum antinumuComp = decomp->MC_AntiNumuComponent();
    Info("GetNDComponents", "Loaded AntiNumuComponent...");
    Spectrum ncTotComp = decomp->MC_NCTotalComponent();
    Info("GetNDComponents", "Loaded NCTotalComponent...");
    Spectrum dataComp = decomp->Data_Component();
    Info("GetNDComponents", "Loaded Data_Component...");

    auto beam = nueComp + antinueComp;
    auto numu = numuComp + antinumuComp;

    // Get correct POT for horn current
    double ndPOT;

    if (isFHC)
      ndPOT = kAna2018FHCPOT;
    else
      ndPOT = kAna2018RHCPOT;

    DataMCComponents ret;
    ret.beam = beam.ToTH1(ndPOT);
    std::cout << "beam component: " << ret.beam->Integral()
      << ", in " << ret.beam->GetNbinsX() << " bins" << std::endl;
    ret.mcTotal = reinterpret_cast<TH1*>(ret.beam->Clone());

    ret.numu = numu.ToTH1(ndPOT);
    std::cout << "numu component: " << ret.numu->Integral()
      << ", in " << ret.numu->GetNbinsX() << " bins" << std::endl;
    ret.mcTotal->Add(ret.numu);

    ret.nc = ncTotComp.ToTH1(ndPOT);
    std::cout << "nc component: " << ret.nc->Integral()
      << ", in " << ret.nc->GetNbinsX() << " bins" << std::endl;
    ret.mcTotal->Add(ret.nc);
    std::cout << "Total MC: " << ret.mcTotal->Integral()
      << ", in " << ret.mcTotal->GetNbinsX() << " bins" << std::endl;

    ret.data = dataComp.ToTH1(ndPOT);
    std::cout << "data component: " << ret.data->Integral()
      << ", in " << ret.beam->GetNbinsX() << " bins" << std::endl;

    ret.bkg = reinterpret_cast<TH1*>(ret.nc->Clone());
    ret.sig = 0;
    ret.wrongSign = 0;
    ret.tau = 0;
    ret.nuetonumu = 0;
    ret.nuetonutau = 0;
    ret.numutonumu = 0;
    ret.numutonutau = 0;

    std::cout << "Loaded ND components" << std::endl;

    // Add formatting
    // TODO: it's quite unintuitive that all formatting goes through
    // DefaultFormatNue, even for Numu!
    DefaultFormatNue(ret);

    return ret;
  }

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

  struct DataMCComponents GetFDMCComponents(
    osc::IOscCalc * calc, IPrediction * pred_no, TString
    output_name="nue", int linestyle=1, bool bkgdetails=false)
  {
    DataMCComponents ret;
    // Set some defaults for the components.
    Flavors::Flavors_t MySig, MyBeam, MyNuMu, MyNC, MyTau;
    MySig = MyBeam = MyNuMu = MyNC = MyTau = Flavors::kAll;
    // Change the components depending on what I'm looking at.
    if (!isNumuAna)
    { // Nu e
      MySig  = Flavors::kNuMuToNuE;
      MyBeam = Flavors::kNuEToNuE;
      MyNuMu = Flavors::kAllNuMu;
      MyNC   = Flavors::kAll;
      MyTau  = Flavors::kAllNuTau;
    }
    else
    { // Nu mu
      MySig  = Flavors::kNuMuToNuMu;
      MyBeam = Flavors::kAllNuE;
      MyNuMu = Flavors::kNuEToNuMu;
      MyNC   = Flavors::kAll;
      MyTau  = Flavors::kAllNuTau;
    }

    // Specify sign depending on polarity
    Sign::Sign_t sign, wrongSign;
    sign = wrongSign = Sign::kBoth;
    if (!isNumuAna) {
      if (isFHC)
  {
    sign = Sign::kNu;
    wrongSign = Sign::kAntiNu;
  }
      else
  {
    sign = Sign::kAntiNu;
    wrongSign = Sign::kNu;
  }
    }
    const double pot = (isFHC) ? kAna2018FHCPOT:kAna2018RHCPOT;

    // Predict components
    ret.sig = pred_no->PredictComponent(
      calc, MySig, Current::kCC, sign).ToTH1(pot);
    ret.wrongSign = pred_no->PredictComponent(
      calc, MySig, Current::kCC, wrongSign).ToTH1(pot);
    ret.beam = pred_no->PredictComponent(
      calc, MyBeam, Current::kCC, Sign::kBoth).ToTH1(pot);
    ret.numu = pred_no->PredictComponent(
      calc, MyNuMu, Current::kCC, Sign::kBoth).ToTH1(pot);
    ret.nc = pred_no->PredictComponent(
      calc, MyNC,   Current::kNC, Sign::kBoth).ToTH1(pot);
    ret.tau = pred_no->PredictComponent(
      calc, MyTau,  Current::kCC, Sign::kBoth).ToTH1(pot);

    ret.bkg = (pred_no->Predict(calc) - pred_no->PredictComponent(
      calc, MySig, Current::kCC, sign)).ToTH1(pot);

    ret.mcTotal = reinterpret_cast<TH1*>(ret.bkg->Clone());
    ret.mcTotal->Add(ret.sig);

    if(bkgdetails)
    {
      ret.nuetonumu = pred_no->PredictComponent(
        calc, Flavors::kNuEToNuMu, Current::kCC, Sign::kBoth).ToTH1(pot);
      ret.numutonumu  = pred_no->PredictComponent(
        calc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kBoth).ToTH1(pot);
      ret.nuetonutau  = pred_no->PredictComponent(
        calc, Flavors::kNuEToNuTau, Current::kCC, Sign::kBoth).ToTH1(pot);
      ret.numutonutau = pred_no->PredictComponent(
        calc, Flavors::kNuMuToNuTau, Current::kCC, Sign::kBoth).ToTH1(pot);
    }
    else
    {
      ret.nuetonumu = 0; ret.numutonumu = 0;
      ret.nuetonutau = 0; ret.numutonutau = 0;
    }

    ret.data = 0;
    DefaultFormatNue(ret, linestyle);

    return ret;
  }

//////////////////////////////////////////////////////////////////////
// Compare ND Data MC
//////////////////////////////////////////////////////////////////////

  void CompareNDDataOneMC(DataMCComponents h_no, DataMCComponents h_sh,
                          TString plottitle, TString out_name, TString pidtag,
                          AxisType pidaxis, bool printtable=true)
  {
    h_sh.bkg->SetLineColor(kTotalMCErrorBandColor);
    h_sh.bkg->SetLineStyle(1);

    h_no.bkg->SetMaximum(FindLimitY({h_no.mcTotal, h_no.data}));
    h_no.bkg->SetTitle(plottitle);

    TGaxis::SetMaxDigits(3);
    h_no.bkg->GetXaxis()->SetDecimals();

    auto leg1 = CustomLegend(
      {h_no.data, h_no.mcTotal, h_sh.mcTotal},
      {"ND Data","Nominal","Shift"},{"epl","l","l"});

    auto leg2 = DefaultNueLegend(false, true, false);

    TH1* mydata = h_no.data;
    std::vector <TH1*> comp_plot = {h_no.mcTotal, h_sh.mcTotal};
    std::vector <TH1*> vnom_plot = {
      h_no.mcTotal, h_no.beam, h_no.nc, h_no.numu};
    std::vector <TH1*> vshi_plot = {h_sh.mcTotal};

    std::vector <TH1*> vnom = {h_no.mcTotal, h_no.beam, h_no.nc, h_no.numu};
    std::vector <TH1*> vshi = {h_sh.mcTotal, h_sh.beam, h_sh.nc, h_sh.numu};
    std::vector<TString> labels = {
      "Total MC", "Beam $\\nu_{e}$", "NC", "$\\nu_{\\mu}$ CC"};

    // Formatting for shift -->  do this before re-binning for Numu
    h_sh.mcTotal->SetLineStyle(2);
    h_sh.mcTotal->SetLineColor(kTotalMCColor+1);

    if(pidaxis == kNumuOpt)
    {
      leg2 = DefaultNumuLegend(false, true, false);
      mydata = HistNumuOpt(h_no.data);
      comp_plot = {HistNumuOpt(h_no.mcTotal), HistNumuOpt(h_sh.mcTotal)};
      vnom_plot = {HistNumuOpt(h_no.mcTotal), HistNumuOpt(h_no.numu),
                   HistNumuOpt(h_no.nc), HistNumuOpt(h_no.beam)};
      vshi_plot = {HistNumuOpt(h_sh.mcTotal)};

      vnom = {HistNumuOpt(h_no.mcTotal), HistNumuOpt(h_no.numu),
              HistNumuOpt(h_no.nc), HistNumuOpt(h_no.beam)};
      vshi = {HistNumuOpt(h_sh.mcTotal), HistNumuOpt(h_sh.numu),
              HistNumuOpt(h_sh.nc), HistNumuOpt(h_sh.beam)};
      labels = {"Total MC", "$\\nu_{\\mu}$ CC", "NC", "$\\nu_{e}$ CC"};
    }

    PlotNDDataTotalMCComparison(
      mydata, comp_plot, leg1, pidtag, out_name, printratioplots, pidaxis);

    leg2->AddEntry(h_sh.mcTotal, "Shifted", "l");

    PlotMCComponentsComparison(
      vnom_plot, vshi_plot, leg2, pidtag,
      out_name, printratioplots, false, pidaxis);

    if(printtable)
    {
      ComparisonTableOne(vnom, vshi, labels, out_name);

      int nbins = 1;
      if (pidaxis== kNue3bin)
        nbins = 2;  // No peripheral bin for ND
      ComparisonTableOneNbins(vnom, vshi, labels, out_name, nbins);
    }
  }
//////////////////////////////////////////////////////////////////////

  void CompareNDDataTwoMC(
    DataMCComponents hnom, DataMCComponents hplu, DataMCComponents hmin,
    TString plottitle, TString out_name, TString pidtag,
    AxisType pidaxis, bool printtable=true)
  {
    Info("CompareNDDataTwoMC",
         TString::Format("pidaxis = %d", pidaxis));
    // std::cout << "DEBUG[CompareNDDataTwoMC]: hnom.bkg "
    //   << hnom.bkg->Integral() << std::endl;
    // std::cout << "DEBUG[CompareNDDataTwoMC]: hplu.bkg "
    //   << hplu.bkg->Integral() << std::endl;
    // std::cout << "DEBUG[CompareNDDataTwoMC]: hmin.bkg "
    //   << hmin.bkg->Integral() << std::endl;
    //
    // std::cout << "DEBUG[CompareNDDataTwoMC]: hnom.data "
    //   << hnom.data->Integral() << std::endl;
    // std::cout << "DEBUG[CompareNDDataTwoMC]: hplu.data "
    //   << hplu.data->Integral() << std::endl;
    // std::cout << "DEBUG[CompareNDDataTwoMC]: hmin.data "
    //   << hmin.data->Integral() << std::endl;

    // Call here to get correct formatting in legends
    // Need to call again with fixbins as true later
    FormatErrorBand(hplu.mcTotal, hmin.mcTotal);

    hnom.mcTotal->SetMaximum(FindLimitY(
      {hnom.mcTotal, hnom.data, hplu.mcTotal, hmin.mcTotal}));
    hnom.mcTotal->SetTitle(plottitle);

    TGaxis::SetMaxDigits(3);
    hnom.mcTotal->GetXaxis()->SetDecimals();

    auto leg1 = CustomLegend(
      {hnom.data, hnom.mcTotal, hplu.mcTotal},
      {"ND Data","Nominal MC","Syst."}, {"epl","l","f"});
    auto leg2 = DefaultNueLegend(false, true, false);

    TH1 * mydata = hnom.data;
    std::vector <TH1*> comp_plot = {hnom.mcTotal, hplu.mcTotal, hmin.mcTotal};
    std::vector <TH1*> vnom_plot = {
      hnom.mcTotal, hnom.beam, hnom.nc, hnom.numu};
    std::vector <TH1*> vshi_plot = {hplu.mcTotal, hmin.mcTotal};

    std::vector <TH1*> vnom = {hnom.mcTotal, hnom.beam, hnom.nc, hnom.numu};
    std::vector <TH1*> vplu = {hplu.mcTotal, hplu.beam, hplu.nc, hplu.numu};
    std::vector <TH1*> vmin = {hmin.mcTotal, hmin.beam, hmin.nc, hmin.numu};
    std::vector<TString> labels = {
      "Total MC", "Beam $\\nu_{e}$", "NC", "$\\nu_{\\mu}$ CC"};

    if(pidaxis == kNumuOpt)
    {
      leg2 = DefaultNumuLegend(false, true, false);
      mydata = HistNumuOpt(hnom.data);
      comp_plot = {HistNumuOpt(hnom.mcTotal),
                   HistNumuOpt(hplu.mcTotal), HistNumuOpt(hmin.mcTotal)};
      vnom_plot = {
                   HistNumuOpt(hnom.mcTotal),
                   HistNumuOpt(hnom.numu),
                   HistNumuOpt(hnom.nc),
                   HistNumuOpt(hnom.beam)
                 };
      vshi_plot = {HistNumuOpt(hplu.mcTotal), HistNumuOpt(hmin.mcTotal)};

      // vnom = {hnom.mcTotal, hnom.numu, hnom.nc, hnom.beam};
      // vplu = {hplu.mcTotal, hplu.numu, hplu.nc, hplu.beam};
      // vmin = {hmin.mcTotal, hmin.numu, hmin.nc, hmin.beam};
      vnom = {HistNumuOpt(hnom.mcTotal), HistNumuOpt(hnom.numu),
              HistNumuOpt(hnom.nc), HistNumuOpt(hnom.beam)};
      vplu = {HistNumuOpt(hplu.mcTotal), HistNumuOpt(hplu.numu),
              HistNumuOpt(hplu.nc), HistNumuOpt(hplu.beam)};
      vmin = {HistNumuOpt(hmin.mcTotal), HistNumuOpt(hmin.numu),
              HistNumuOpt(hmin.nc), HistNumuOpt(hmin.beam)};
      labels = {"Total MC", "$\\nu_{\\mu}$ CC", "NC", "$\\nu_{e}$ CC"};
    }

    // Format error band on data/MC plot.
    // Do this here because function is shared with single shift case.
    FormatErrorBand(comp_plot[1], comp_plot[2],
                    out_name.Contains("_sig"), true);

    PlotNDDataTotalMCComparison(
      mydata, comp_plot, leg1, pidtag,
      out_name, printratioplots, pidaxis);

    leg2->AddEntry(hplu.mcTotal, "Syst.", "f");

    PlotMCComponentsErrorBand(
      vnom_plot, vshi_plot, leg2, pidtag,
      out_name, printratioplots, false, pidaxis);
    // TODO: seg fault when ratioplot is false
    // --> do we ever not want to plot the ratio?

    if(printtable)
    {
      ComparisonTable(vnom, vplu, vmin, labels, out_name);
      int nbins = 1;
      if (pidaxis== kNue3bin)
        nbins = 2;  // no peripheral bin for ND
      ComparisonTableNbins(vnom, vplu, vmin, labels, out_name, nbins);
    }
  }

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

  enum EModExtrapComps
  {
    kEEextrap, kEEAntiextrap,
    kMMextrap, kMMAntiextrap,
    kMEextrap, kMEAntiextrap,
    kEMextrap, kEMAntiextrap,
    // NC
    kNCTotalextrap, kNCextrap, kNCAntiextrap,
    // End NC
    kMTextrap, kMTAntiextrap,
    kETextrap, kETAntiextrap,
  };

  const IDecomp* GetDecomp(IPrediction* prediction,
                           EModExtrapComps modExtrapComp)
  {
    // First cast IPrediction as a PredictionExtrap object and check it extist
    std::cout << "Checking if IPrediction exists" << std::endl;
    if (!prediction)
      abort();

    PredictionExtrap* predExtrap;

    // If numu then it's simmple;
    if (isNumuAna)
    {
      std::cout << abi::__cxa_demangle(typeid(*prediction).name(), 0, 0, 0) << std::endl;
      predExtrap = dynamic_cast<PredictionExtrap*>(prediction);
    }
    else
    {
      // If nue then we have to do some wrangling to get the correct object
      std::cout << abi::__cxa_demangle(typeid(*prediction).name(), 0, 0, 0) << std::endl;
      // This will be a PredictionExtendToPeripheral so cast as that
      PredictionExtendToPeripheral* predExtend = dynamic_cast<PredictionExtendToPeripheral*>(prediction);
      std::cout << abi::__cxa_demangle(typeid(*predExtend).name(), 0, 0, 0) << std::endl;
      std::cout << "Checking if PredictionExtendToPeripheral exists" << std::endl;
      if (!predExtend)
        abort();

      // Now we want to get the core part of this
      IPrediction* predCore = predExtend->GetCore();
      std::cout << "Checking if IPrediction (core) exists" << std::endl;
      if (!predCore)
        abort();

      // Now we want to cast this as a PredictionExtrap
      std::cout << abi::__cxa_demangle(typeid(*predCore).name(), 0, 0, 0) << std::endl;
      predExtrap = dynamic_cast<PredictionExtrap*>(predCore);
      std::cout << "Checking if PredictionExtrap exists" << std::endl;
      std::cout << abi::__cxa_demangle(typeid(*predExtrap).name(), 0, 0, 0) << std::endl;
    }
    if (!predExtrap)
      abort();
    // Get the Extrap object
    IExtrap* extrap = predExtrap->GetExtrap();

    // Cast as modular extrap
    ModularExtrap* modExtrap = dynamic_cast<ModularExtrap*>(extrap);

    // Get mod extrap components
    std::vector<ModularExtrapComponent*> modExtrapComps =
      modExtrap->GetModExtrapComponents();

    // Get decomp for EEextrap
    // Use vector index to select appropriate set of extrap components
    std::cout << "*************************************" << std::endl;
    std::cout << "Getting modular extrap component: "
      << modExtrapComp << std::endl;
    std::cout << "*************************************" << std::endl;
    const IDecomp* decomp = modExtrapComps[modExtrapComp]->GetDecomp();

    return decomp;
  }


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

  void CompareNDDataMCFromVector(
    PredictionSystJoint2018* predictionSyst, const ISyst* syst,
    EModExtrapComps modExtrapComp, TString plottitle, TString out_name,
    TString tag, AxisType pidaxis, bool printtable=true)
  {
    Info("CompareNDDataMCFromVector",
         TString::Format("pidaxis = %d", pidaxis));

    // Get nominal from PredictionSystJoint2018
    std::cout << std::endl;
    std::cout << "Getting Nominal prediction" << std::endl;
    auto pred_no = predictionSyst->GetNominalPrediction();
    std::cout << "In ComparePredictionsFromVector -- pred_no is " << pred_no << std::endl;
    // Get Nominal decomp
    // Option 1: through decomp
    const IDecomp* nomDecomp = GetDecomp(pred_no, modExtrapComp);
    DataMCComponents nomNDComps = GetNDCompsFromDecomp(nomDecomp);
    std::cout << "Got nomDecomp" << std::endl;
    std::cout << "Got nomNDComps" << std::endl;
    // Get shifted predictions from PredictionSystJoint2018
    const PredictionInterp::ShiftedPreds& shiftedPred =
      predictionSyst->GetShiftedPrediction(syst);

    std::cout << std::endl;
    std::cout << "ShiftedPreds name: " << shiftedPred.systName << std::endl;

    int numSh = shiftedPred.preds.size();
    std::cout << "Number of shifts: " << numSh << std::endl;
    // for (auto shift:shiftedPred.preds) { std::cout << "Func shift: " << shift << std::endl; }

    std::vector<IPrediction*> pred_sh = shiftedPred.preds;

    // create vecor of components
    std::vector<DataMCComponents> comps;
    for (auto pred:pred_sh)
    {
      // Option 1: through decomp
      const IDecomp* shiftedDecomp = GetDecomp(pred, modExtrapComp);
      std::cout << "Got shiftedDecomp" << std::endl;
      DataMCComponents shiftedComps = GetNDCompsFromDecomp(shiftedDecomp);
      std::cout << "Got shiftedComps" << std::endl;
      std::cout << std::endl;
      comps.push_back(shiftedComps);
    }

    std::cout << "comps has length: " << comps.size() << std::endl;

    if (numSh == 3)
    {
      // shift 0 = -1
      // shift 1 = Nominal
      // shift 2 = +1
      std::cout << "Entering shift = 3" << std::endl;
      CompareNDDataTwoMC(
        nomNDComps, comps[2], comps[0], plottitle,
        out_name+"_1sigma", "#pm1#sigma", pidaxis, printtable);
    }
    else if (numSh == 5)
    {
      // shift 1 = -1
      // shift 2 = Nominal
      // shift 3 = +1
      std::cout << "Entering shift = 5" << std::endl;
      CompareNDDataTwoMC(
        nomNDComps, comps[3], comps[1], plottitle,
        out_name+"_1sigma", "#pm1#sigma", pidaxis, printtable);
    }
    else if (numSh == 2)
    {
      // shift 0 = Nominal
      // shift 1 = +1
      std::cout << "Entering shift = 2" << std::endl;
      CompareNDDataOneMC(
        nomNDComps, comps[1], plottitle, out_name+"_1sigma",
        "#plus1#sigma", pidaxis, printtable);
    }
  }

//////////////////////////////////////////////////////////////////////
// Compare FD Predictions
//////////////////////////////////////////////////////////////////////

  void CompareOneShiftPred(IPrediction * pred_no, IPrediction * pred_sh,
         TString plottitle, TString out_name, TString pidtag,
         bool printtable=true, AxisType pidaxis = kNue1bin, bool PrintParams=false){

    auto calc = GetCalculator (usedumb);
    auto h_no = GetFDMCComponents (calc, pred_no, out_name);
    auto h_sh = GetFDMCComponents (calc, pred_sh, out_name, 2);

    h_no.sig->SetMaximum(FindLimitY({h_no.sig,h_sh.sig}));
    h_no.bkg->SetMaximum(FindLimitY({h_no.bkg,h_sh.bkg}));
    h_no.sig->SetTitle(plottitle);
    h_no.bkg->SetTitle(plottitle);
    TGaxis::SetMaxDigits(3);

    auto leg = DefaultNueLegend(false,true,false);
    if (isNumuAna) leg = DefaultNumuLegend(false,true,false);
    leg->AddEntry(DummyGraph(kGray+1,2),"Shifted","l");

    std::vector<TString> legnames = {"#nu_{e} Signal","Shifted"};
    if (isNumuAna) legnames = {"#nu_{#mu} Signal","Shifted"};
    auto leg2 = CustomLegend({h_no.sig, h_sh.sig},legnames,{"l","l"});

    if(pidaxis == kNumuOpt)
    { //fix the plot but don't mess up the tables

      PlotMCComponentsComparison({ HistNumuOpt(h_no.bkg), HistNumuOpt(h_no.beam), HistNumuOpt(h_no.nc), HistNumuOpt(h_no.numu) },
   { HistNumuOpt(h_sh.bkg), HistNumuOpt(h_sh.beam), HistNumuOpt(h_sh.nc), HistNumuOpt(h_sh.numu) },
   leg,pidtag,out_name+"_bkg",printratioplots,false,pidaxis);
      PlotMCComponentsComparison({ HistNumuOpt(h_no.sig) },
   { HistNumuOpt(h_sh.sig) },
   leg2,pidtag,out_name+"_sig",printratioplots,false,pidaxis);
    }
    else
    {
      PlotMCComponentsComparison({ h_no.bkg, h_no.beam, h_no.nc, h_no.numu },
         { h_sh.bkg, h_sh.beam, h_sh.nc, h_sh.numu },
         leg,pidtag, out_name+"_bkg", printratioplots, false, pidaxis);
      PlotMCComponentsComparison({ h_no.sig },
         { h_sh.sig },
         leg2,pidtag,out_name+"_sig",printratioplots,false,pidaxis);
    }

    if(PrintParams) PrintOscilationParameters (calc, usedumb);

    if(printtable){
      std::vector<TString> labels = {"$\\nu_e$ signal","Tot beam bkg","Beam $\\nu_{e}$ CC","NC", "$\\nu_{\\mu}$ CC", "$\\nu_{\\tau}$ CC"};

      if (isNumuAna) {
  labels = {"$\\nu_{\\mu}$ signal","Tot beam bkg","$\\nu_{e}$ CC","NC", "$\\nu_{\\mu} $App", "$\\nu_{\\tau}$ CC"};
      }
      ComparisonTableOne ({h_no.sig, h_no.bkg, h_no.beam, h_no.nc, h_no.numu, h_no.tau},
        {h_sh.sig, h_sh.bkg, h_sh.beam, h_sh.nc, h_sh.numu, h_sh.tau},
        labels,
        out_name
  );
      int nbins = 1;
      if (pidaxis== kNue3bin) nbins = 3;
      ComparisonTableOneNbins({h_no.sig, h_no.bkg, h_no.beam,h_no.nc,h_no.numu,h_no.tau},
            {h_sh.sig, h_sh.bkg, h_sh.beam,h_sh.nc,h_sh.numu,h_sh.tau},
            labels,
            out_name,
            nbins);
    }//end if table
  }
//////////////////////////////////////////////////////////////////////

  void CompareTwoShiftPred(
    IPrediction * pred_no, IPrediction * pred_pl, IPrediction * pred_mi,
    TString plottitle, TString out_name, TString pidtag,
    bool printtable=true, AxisType pidaxis = kNue1bin, bool PrintParams=false)
  {

    auto calc = GetCalculator(usedumb);

    auto hnom = GetFDMCComponents(calc, pred_no, out_name);
    auto hplu = GetFDMCComponents(calc, pred_pl, out_name);
    auto hmin = GetFDMCComponents(calc, pred_mi, out_name);

    FormatErrorBand(hplu.bkg, hmin.bkg, false);
    hnom.bkg->SetMaximum(FindLimitY({hnom.bkg, hplu.bkg, hmin.bkg}));
    hnom.bkg->SetTitle(plottitle);

    FormatErrorBand(hplu.sig, hmin.sig, true);
    hnom.sig->SetMaximum(FindLimitY({hnom.sig, hplu.sig, hmin.sig}));
    hnom.sig->SetTitle(plottitle);

    TGaxis::SetMaxDigits(3);

    TLegend * leg, *leg2;
    if( pidaxis == kNumuOpt || pidaxis == kNumuOther)
    {
      leg = DefaultNumuLegend(false, true, false);
      leg2 = CustomLegend({hnom.sig, hplu.sig},
                          {"#nu_{#mu} Signal","Shifted"},
                          {"l","f"});
    }
    else
    {
      leg = DefaultNueLegend(false, true, false);
      leg2 = CustomLegend({hnom.sig, hplu.sig},
                          {"#nu_{e} Signal","Shifted"},
                          {"l","f"});
    }
    leg->AddEntry(hplu.bkg,"Syst.","f");

    std::vector <TH1*> nom_bkg_plot = {hnom.bkg, hnom.beam, hnom.nc, hnom.numu, hnom.tau};
    std::vector <TH1*> shi_bkg_plot = {hplu.bkg, hmin.bkg};
    std::vector <TH1*> nom_sig_plot = {hnom.sig};
    std::vector <TH1*> shi_sig_plot = {hplu.sig, hmin.sig};

    std::vector <TH1*> vnom = {
      hnom.sig, hnom.bkg, hnom.beam, hnom.nc, hnom.numu, hnom.tau};
    std::vector <TH1*> vplu = {
      hplu.sig, hplu.bkg, hplu.beam, hplu.nc, hplu.numu, hplu.tau};
    std::vector <TH1*> vmin = {
      hmin.sig, hmin.bkg, hmin.beam, hmin.nc, hmin.numu, hmin.tau};

    std::vector <TString> labels = {
      "$\\nu_e$ signal", "Tot beam bkg", "Beam $\\nu_{e}$ CC",
      "NC", "$\\nu_{\\mu}$ CC", "$\\nu_{\\tau}$ CC"};
    if(pidaxis == kNumuOpt || pidaxis == kNumuOther)
    {
      labels = {"$\\nu_{\\mu}$ signal", "Tot beam bkg", "NC",
                "$\\nu_{\\mu}$ App", "$\\nu_{e}$ CC", "$\\nu_{\\tau}$ CC"};
      vnom = {hnom.sig, hnom.bkg, hnom.nc, hnom.numu, hnom.beam, hnom.tau};
      vplu = {hplu.sig, hplu.bkg, hplu.nc, hplu.numu, hnom.beam, hplu.tau};
      vmin = {hmin.sig, hmin.bkg, hmin.nc, hmin.numu, hnom.beam, hmin.tau};
    }
    if(pidaxis == kNumuOpt)
    {
      nom_bkg_plot = {HistNumuOpt(hnom.bkg), HistNumuOpt(hnom.beam), HistNumuOpt(hnom.nc), HistNumuOpt(hnom.numu), HistNumuOpt(hnom.tau)};
      shi_bkg_plot = {HistNumuOpt(hplu.bkg), HistNumuOpt(hmin.bkg)};
      nom_sig_plot = {HistNumuOpt(hnom.sig)};
      shi_sig_plot = {HistNumuOpt(hplu.sig), HistNumuOpt(hmin.sig)};
    }

    PlotMCComponentsErrorBand(nom_bkg_plot, shi_bkg_plot,
            leg,  pidtag, out_name + "_bkg", printratioplots, false, pidaxis);
    PlotMCComponentsErrorBand(nom_sig_plot, shi_sig_plot,
            leg2, pidtag, out_name + "_sig", printratioplots, false, pidaxis);

    if(PrintParams) PrintOscilationParameters (calc, usedumb);

    if(printtable){

      ComparisonTable (vnom, vplu, vmin, labels, out_name);

      int nbins = 1;
      if (pidaxis== kNue3bin) nbins = 3;
      ComparisonTableNbins (vnom, vplu, vmin, labels, out_name, nbins);
    }
  }

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

  void ComparePredictionsFromVector(
    PredictionSystJoint2018* predictionSyst, const ISyst* syst,
    TString plottitle, TString out_name,
    bool printtable=true, AxisType pidaxis=kNue3bin)
  {
    // Get nominal from PredictionSystJoint2018
    auto pred_no = predictionSyst->GetNominalPrediction();
    //std::cout << "In ComparePredictionsFromVector -- pred_no is " << pred_no << std::endl;

    // Get shifted predictions from PredictionSystJoint2018
    const PredictionInterp::ShiftedPreds& shiftedPred =
      predictionSyst->GetShiftedPrediction(syst);

    //std::cout << "ShiftedPreds name: " << shiftedPred.systName << std::endl;

    int numSh = shiftedPred.preds.size();
    //std::cout << "Number of shifts: " << numSh << std::endl;
    //for (auto shift:shiftedPred.preds) { std::cout << "Func shift: " << shift << std::endl; }

    std::vector<IPrediction*> pred_sh = shiftedPred.preds;

    if (numSh == 3)
    {
      // shift 0 = -1
      // shift 1 = Nominal
      // shift 2 = +1
      //std::cout << "Passing; " << pred_no << ", " << pred_sh[0] << ", " << pred_sh[2] << std::endl;
      CompareTwoShiftPred(
        pred_no, pred_sh[0], pred_sh[2], plottitle,
        out_name+"_1sigma", "#pm1#sigma", printtable, pidaxis);
    }
    else if (numSh == 5)
    {
      // shift 1 = -1
      // shift 2 = Nominal
      // shift 3 = +1
      //std::cout << "Passing; " << pred_no << ", " << pred_sh[1] << ", " << pred_sh[3] << std::endl;
      CompareTwoShiftPred(
        pred_no, pred_sh[1], pred_sh[3], plottitle,
        out_name+"_1sigma", "#pm1#sigma", printtable, pidaxis);
    }
    else if (numSh == 2)
    {
      // shift 0 = Nominal
      // shift 1 = +1
      //std::cout << "Passing; " << pred_no << ", " << pred_sh[1] << std::endl;
      CompareOneShiftPred(
        pred_no, pred_sh[1], plottitle, out_name+"_1sigma",
        "#plus1#sigma", printtable, pidaxis);
    }
  }
}