SystMaker.cxx

Go to the documentation of this file.

#include "NuXAna/Analysis/SystMaker.h"

#include "OscLib/IOscCalc.h"

#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Extrap/IExtrap.h"
#include "CAFAna/Analysis/Calcs.h"

#include "NuXAna/Prediction/PredictionGeneratorNuX.h"
#include "NuXAna/Prediction/FDPredictionSterile.h"
#include "NuXAna/Prediction/NDPredictionSterile.h"
#include "NuXAna/Prediction/PredictionSterile.h"

#include "TCanvas.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TObjString.h"
#include "TKey.h"

using std::cout;
using std::cerr;
using std::endl;
using std::string;
using std::vector;
using std::map;
using std::pair;
using std::make_pair;
using std::find;
using std::move;
using std::unique_ptr;
using std::make_unique;
using std::showpos;
using std::ostream;
using std::ostringstream;
using std::runtime_error;

using osc::IOscCalc;

using ana::covmx::Sample;

namespace ana {

  // -------------------------------------------------------------------------------------
  IPrediction* GetPrediction(SystPredType predType,
           const HistAxis* axis, const HistAxis* numuAxis,
           const Cut* FDCut, const Cut* NDCut, const Cut* numuCut,
           const SystShifts* shiftData, const SystShifts* shiftMC,
           const Var* weight,
           Loaders* loaders, Loaders* loadersND) {

    IPrediction* pred = nullptr;

    if (predType == SystPredType::kExtrap) {
      SterileGenerator genExtrap(*axis, *numuAxis,
               *FDCut, *NDCut, *numuCut,
               *shiftData, *weight);
      if (loadersND)
        pred = genExtrap.Generate(*loadersND, *loaders, *shiftMC).release();
      else
        pred = genExtrap.Generate(*loaders, *shiftMC).release();
    }

    else if (predType == SystPredType::kFD) {
      FDPredictionGenerator genFD(*axis, *FDCut, *shiftData, *weight);
      pred = genFD.Generate(*loaders, *shiftMC).release();
    }

    else if (predType == SystPredType::kND) {
      NDPredictionGenerator genND(*axis, *NDCut, *shiftData, *weight);
      pred = genND.Generate(*loaders, *shiftMC).release();
    }

    return pred;

  } // function GetPrediction

  // --------------------------------------------------------------------------
  TH1* GetSpectrum(IPrediction* pred, bool signal, const Sample sample,
    IOscCalc* calc) {

    // osc::IOscCalc* calc = DefaultSterileCalc(4);
    // calc->SetDm(4, 10.);
    // calc->SetAngle(2, 4, M_PI/4);
    // calc->SetAngle(3, 4, M_PI/2);

    // If NC selection, signal is NC and background is CC. Nice and easy.
    if (sample.selection == covmx::Selection::kNC) {
      return pred->PredictComponent(calc, Flavors::kAll,
        signal ? Current::kNC : Current::kCC, Sign::kBoth).ToTH1(sample.GetPOT());
    }

    // Get beam polarity and CC signal channel for this sample
    Sign::Sign_t sign = (sample.polarity==covmx::kFHC) ? Sign::kNu : Sign::kAntiNu;
    Flavors::Flavors_t flavor;
    if (sample.selection == covmx::kCCNumu)
      flavor = Flavors::kNuMuToNuMu;
    else if (sample.selection == covmx::kCCNue)
      flavor = Flavors::kNuMuToNuE;
    else assert(false and "Selection not recognised!");

    // Signal is the specified CC channel
    if (signal) {
      return pred->PredictComponent(calc, flavor, Current::kCC,
        sign).ToTH1(sample.GetPOT());
    }
    // Background is everything else
    else {
      // Start with correct sign backgrounds
      Spectrum ret = pred->PredictComponent(calc, Flavors::Flavors_t(~flavor),
        Current::kCC, sign);
      // CC wrong sign backgrounds
      ret += pred->PredictComponent(calc, Flavors::kAll,
        Current::kCC, Sign::Sign_t(~sign));
      // NC everything
      ret += pred->PredictComponent(calc, Flavors::kAll, Current::kNC,
        Sign::kBoth);
      return ret.ToTH1(sample.GetPOT());
    }

  } // function GetSpectrum

  // --------------------------------------------------------------------------
  string MakeLatexCommandName(const string& instring){

    TString mystring(instring.c_str());

    vector<TString> in = {"0","1","2","3","4","5",
             "6","7","8","9","_","#"," ",".","-","(",")"};
    vector<TString> out = {"zero","one","two","three","four","five",
              "six","seven","eight","nine","","","","","XX","",""};
    for(size_t i = 0; i < in.size(); ++i)
      mystring.ReplaceAll(in[i],out[i]);

    return string(mystring.Data());
  }

  // -------------------------------------------------------------------------------------
  /// SystMakerShiftSigma constructor
  SystMakerShiftSigma::SystMakerShiftSigma(string name, int sigma) :
    fName(name),
    fSigma(sigma) {
  }

  // -------------------------------------------------------------------------------------
  /// SystMakerShiftSigma destructor
  SystMakerShiftSigma::~SystMakerShiftSigma() {
    if (fPrediction)             delete fPrediction;
    if (fShiftedSpectrum.first)  delete fShiftedSpectrum.first;
    if (fShiftedSpectrum.second) delete fShiftedSpectrum.second;
    if (fShiftedRatio.first)     delete fShiftedRatio.first;
    if (fShiftedRatio.second)    delete fShiftedRatio.second;
  }

  // -------------------------------------------------------------------------------------
  /// Print the sigmas for a given shift
  void SystMakerShiftSigma::PrintSigma(ostream& os) {
    os << "    " << showpos << fSigma << ": "
       << fShift.first << " (signal), "
       << fShift.second << " (background)"
       << endl;
  }

  // -------------------------------------------------------------------------------------
  void SystMakerShiftSigma::SaveTo(TDirectory* dir) const {

    TDirectory* tmp = gDirectory;

    dir->cd();
    TObjString("SystMakerShiftSigma").Write("type");

    TObjString(fName.c_str()).Write("name");
    TVectorD sigma(1);
    sigma[0] = fSigma;
    sigma.Write("sigma");

    fPrediction->SaveTo(dir, "pred");
    TH1* shiftSig = fShiftedSpectrum.first;
    TH1* shiftBkg = fShiftedSpectrum.second;
    if (shiftSig)
      shiftSig->Write("shiftSig");
    if (shiftBkg)
      shiftBkg->Write("shiftBkg");
    TH1* shiftRatSig = fShiftedRatio.first;
    TH1* shiftRatBkg = fShiftedRatio.second;
    if (shiftRatSig)
      shiftRatSig->Write("shiftRatSig");
    if (shiftRatBkg)
      shiftRatBkg->Write("shiftRatBkg");

    TVectorD v(2);
    v[0] = fShift.first;
    v[1] = fShift.second;
    v.Write("shift");

    tmp->cd();

  }

  // -------------------------------------------------------------------------------------
  SystMakerShiftSigma* SystMakerShiftSigma::LoadFrom(TDirectory* dir) {

    DontAddDirectory guard;

    TObjString* tag = (TObjString*)dir->Get("type");
    assert(tag);
    assert(tag->GetString() == "SystMakerShiftSigma");

    TObjString* name = (TObjString*)dir->Get("name");;
    const TVectorD sigma = *(TVectorD*)dir->Get("sigma");
    SystMakerShiftSigma* systMakerShiftSigma
      = new SystMakerShiftSigma(name->GetString().Data(), sigma[0]);

    tag = (TObjString*)dir->Get("pred/type");
    assert(tag);
    if (tag->GetString() == "NDPredictionSterile")
      systMakerShiftSigma->fPrediction = NDPredictionSterile::LoadFrom(dir, "pred").release();
    else if (tag->GetString() == "FDPredictionSterile")
      systMakerShiftSigma->fPrediction = FDPredictionSterile::LoadFrom(dir, "pred").release();
    else if (tag->GetString() == "PredictionSterile")
      systMakerShiftSigma->fPrediction = PredictionSterile::LoadFrom(dir, "pred").release();
    else assert(false && "Prediction type not recognised!");

    TH1* shiftSig = (TH1*)dir->Get("shiftSig");
    TH1* shiftBkg = (TH1*)dir->Get("shiftBkg");
    if (shiftSig and shiftBkg)
      systMakerShiftSigma->fShiftedSpectrum = make_pair(shiftSig, shiftBkg);
    TH1* shiftRatSig = (TH1*)dir->Get("shiftRatSig");
    TH1* shiftRatBkg = (TH1*)dir->Get("shiftRatBkg");
    if (shiftRatSig and shiftRatBkg)
      systMakerShiftSigma->fShiftedRatio = make_pair(shiftRatSig, shiftRatBkg);

    const TVectorD v = *(TVectorD*)dir->Get("shift");
    systMakerShiftSigma->fShift = make_pair(v[0], v[1]);

    return systMakerShiftSigma;

  }

  // -------------------------------------------------------------------------------------
  SystMakerShift::SystMakerShift(string name, string label) :
    fName(name),
    fLabel(label)
  {}

  // -------------------------------------------------------------------------------------
  /// Default destructor
  SystMakerShift::~SystMakerShift()
  {
    for (auto sigma : fSigmas) {
      if (sigma.second) delete sigma.second;
    }
  }

  // -------------------------------------------------------------------------------------
  vector<int> SystMakerShift::GetSigmas() {
    vector<int> sigmas;
    for (const auto& sigma : fSigmas)
      sigmas.push_back(sigma.first);
    return sigmas;
  }

  // -------------------------------------------------------------------------------------
  void SystMakerShift::PrintShift(ostream& os) {
    os << endl << "  " << fLabel << " shifts:" << endl;
    for (const auto& sigma : fSigmas)
      sigma.second->PrintSigma(os);
    return;
  }

  // -------------------------------------------------------------------------------------
  bool SystMakerShift::CheckShifts() {
    // Require -1, +1 shifts
    bool good = true;
    if (!fSigmas.count(-1) or !fSigmas.count(+1))
      good = false;
    return good;
  }

  // -------------------------------------------------------------------------------------
  TCanvas* SystMakerShift::DrawSig(TH1* nominal) {

    // Drawing tools ------------
    TCanvas* canv = new TCanvas("canv", "", 800, 800);
    TPad* padSpecSigTotal  = new TPad("padSpecSigTotal",  "", 0., 0.375, 1., 1.);
    TPad* padRatioSigTotal = new TPad("padRatioSigTotal", "", 0., 0.,    1., 0.375);
    TLegend* leg = new TLegend(0.5, 0.5, 0.88, 0.88);
    map<int, TLegendEntry*> legEntries;

    // Draw nominal ------------
    TH1* hNom = (TH1*)nominal->Clone("hNom");
    for (size_t bin = 1; bin <= size_t(hNom->GetNbinsX()); ++bin)
      hNom->SetBinContent(bin, hNom->GetBinContent(bin)/hNom->GetBinWidth(bin));
    hNom->GetYaxis()->SetTitle("Events / GeV");
    hNom->SetLineStyle(kSolid);
    hNom->SetLineColor(kBlack);
    hNom->SetLineWidth(2);
    padSpecSigTotal->cd();
    hNom->Draw("hist");

    // ratio
    TH1* nomRatio = (TH1*)nominal->Clone("nomRatio");
    nomRatio->GetYaxis()->SetTitle("Shifted / Nominal");
    for (Int_t bin = 0; bin <= nomRatio->GetNbinsX(); ++bin)
      nomRatio->SetBinContent(bin, 1.);
    padRatioSigTotal->cd();
    nomRatio->Draw("hist");
    legEntries[0] = new TLegendEntry(nominal, "Nominal", "l");
    double ymaxRatio = 0;

    // Draw shifts ------------
    int it = 0;
    for (map<int, SystMakerShiftSigma*>::const_iterator sigma = fSigmas.begin(); sigma != fSigmas.end(); ++sigma, ++it) {
      if (sigma->second->fSigma == 0)
        continue;
      // Total ------------
      TH1* sigmaSig = (TH1*)sigma->second->fShiftedSpectrum.first->Clone();
      for (size_t bin = 1; bin <= size_t(sigmaSig->GetNbinsX()); ++bin)
        sigmaSig->SetBinContent(bin, sigmaSig->GetBinContent(bin)/sigmaSig->GetBinWidth(bin));
      sigmaSig->GetYaxis()->SetTitle("Events / GeV");
      sigmaSig->SetLineStyle(kDashed);
      sigmaSig->SetLineColor(abs(sigma->second->fSigma)+1);
      sigmaSig->SetLineWidth(2);
      padSpecSigTotal->cd();
      sigmaSig->DrawClone("hist same");
      TH1* sigmaRatioSig = (TH1*)sigmaSig->Clone("sigmaRatioSig");
      sigmaRatioSig->Divide(hNom);
      sigmaRatioSig->SetBinContent(0, 1);
      for (size_t i = 1; i <= (size_t)sigmaRatioSig->GetNbinsX(); ++i) {
        double val = sigmaRatioSig->GetBinContent(i);
        if (val == 0) {
          sigmaRatioSig->SetBinContent(i, 1);
          val = 1;
        }
        if (fabs(val-1) > ymaxRatio)
          ymaxRatio = fabs(val-1);
      }
      padRatioSigTotal->cd();
      sigmaRatioSig->DrawClone("hist same");
      legEntries[abs(sigma->second->fSigma)]
        = new TLegendEntry(sigmaSig->Clone(), Form("+/-%d#sigma", sigma->second->fSigma), "l");
    }

    // Draw legend ------------
    for (const auto& legEntry : legEntries)
      leg->AddEntry(legEntry.second->GetObject(),
        legEntry.second->GetLabel(),
        legEntry.second->GetOption());
    padSpecSigTotal->cd();
    leg->Draw();

    nomRatio->SetMinimum(1-(1.1*ymaxRatio));
    nomRatio->SetMaximum(1+(1.1*ymaxRatio));

    // Draw all together ------------
    canv->cd();
    canv->Clear();
    padSpecSigTotal->Draw();
    padRatioSigTotal->Draw();

    return canv;

  } // function SystMakerShift::DrawSig

  // -------------------------------------------------------------------------------------
  TCanvas* SystMakerShift::DrawBkg(TH1* nominal) {

    // Drawing tools ------------
    TCanvas* canv = new TCanvas("canv", "", 800, 800);
    TPad* padSpecBkgTotal  = new TPad("padSpecBkgTotal",  "", 0., 0.375, 1., 1.);
    TPad* padRatioBkgTotal = new TPad("padRatioBkgTotal", "", 0., 0.,    1., 0.375);
    TLegend* leg = new TLegend(0.5, 0.5, 0.88, 0.88);
    map<int, TLegendEntry*> legEntries;

    // Draw nominal ------------
    TH1* hNom = (TH1*)nominal->Clone("hNom");
    for (size_t bin = 1; bin <= size_t(hNom->GetNbinsX()); ++bin)
      hNom->SetBinContent(bin, hNom->GetBinContent(bin)/hNom->GetBinWidth(bin));
    hNom->GetYaxis()->SetTitle("Events / GeV");
    nominal->SetLineStyle(kSolid);
    nominal->SetLineColor(kBlack);
    nominal->SetLineWidth(2);
    padSpecBkgTotal->cd();
    nominal->Draw("hist");

    // ratio
    TH1* nomRatio = (TH1*)nominal->Clone("nomRatio");
    nomRatio->GetYaxis()->SetTitle("Shifted / Nominal");
    for (Int_t bin = 0; bin <= nomRatio->GetNbinsX(); ++bin)
      nomRatio->SetBinContent(bin, 1.);
    padRatioBkgTotal->cd();
    nomRatio->Draw("hist");
    legEntries[0] = new TLegendEntry(nominal, "Nominal", "l");
    double ymaxRatio = 0;

    // Draw shifts ------------
    int it = 0;
    for (map<int, SystMakerShiftSigma*>::const_iterator sigma = fSigmas.begin();
      sigma != fSigmas.end(); ++sigma, ++it) {
      if (sigma->second->fSigma == 0)
        continue;
      // Total ------------
      TH1* sigmaBkg = (TH1*)sigma->second->fShiftedSpectrum.second->Clone();
      for (size_t bin = 1; bin <= size_t(sigmaBkg->GetNbinsX()); ++bin)
        sigmaBkg->SetBinContent(bin, sigmaBkg->GetBinContent(bin)/sigmaBkg->GetBinWidth(bin));
      sigmaBkg->GetYaxis()->SetTitle("Events / GeV");
      sigmaBkg->SetLineStyle(kDashed);
      sigmaBkg->SetLineColor(abs(sigma->second->fSigma)+1);
      sigmaBkg->SetLineWidth(2);
      padSpecBkgTotal->cd();
      sigmaBkg->DrawClone("hist same");
      TH1* sigmaRatioBkg = (TH1*)sigmaBkg->Clone("sigmaRatioBkg");
      sigmaRatioBkg->Divide(hNom);
      sigmaRatioBkg->SetBinContent(0, 1);
      for (size_t i = 1; i <= (size_t)sigmaRatioBkg->GetNbinsX(); ++i) {
        double val = sigmaRatioBkg->GetBinContent(i);
        if (val == 0) {
          sigmaRatioBkg->SetBinContent(i, 1);
          val = 1;
        }
        if (fabs(val-1) > ymaxRatio)
          ymaxRatio = fabs(val-1);
      }
      padRatioBkgTotal->cd();
      sigmaRatioBkg->DrawClone("hist same");
      legEntries[abs(sigma->second->fSigma)]
        = new TLegendEntry(sigmaBkg->Clone(), Form("+/-%d#sigma",
          sigma->second->fSigma), "l");
    }

    // Draw legend ------------
    for (const auto& legEntry : legEntries)
      leg->AddEntry(legEntry.second->GetObject(),
        legEntry.second->GetLabel(),
        legEntry.second->GetOption());
    padSpecBkgTotal->cd();
    leg->Draw();

    nomRatio->SetMinimum(1-(1.1*ymaxRatio));
    nomRatio->SetMaximum(1+(1.1*ymaxRatio));

    // Draw all together ------------
    canv->cd();
    canv->Clear();
    padSpecBkgTotal->Draw();
    padRatioBkgTotal->Draw();

    return canv;

  } // function SystMakerShift::DrawBkg

  // -------------------------------------------------------------------------------------
  void SystMakerShift::DrawShift(TDirectory* outDir, pair<TH1*, TH1*>& nominal) {

    // Out directories ------------
    string dirSigShift = Form("Signal/Shifts/%s", fName.c_str());
    string dirBkgShift = Form("Background/Shifts/%s", fName.c_str());
    string dirSigShift_sigmas = dirSigShift+"/Sigmas/";
    string dirBkgShiftSigmas = dirBkgShift+"/Sigmas/";
    outDir->mkdir(dirSigShift.c_str());
    outDir->mkdir(dirBkgShift.c_str());
    outDir->mkdir(dirSigShift_sigmas.c_str());
    outDir->mkdir(dirBkgShiftSigmas.c_str());

    // Drawing tools ------------
    TCanvas* canv = new TCanvas("canv", "", 800, 800);
    TPad* padSpecSig  = new TPad("padSpecSig",  "", 0., 0.375, 1., 1.);
    TPad* padRatioSig = new TPad("padRatioSig", "", 0., 0.,    1., 0.375);
    TPad* padSpec_bg  = new TPad("padSpec_bg",  "", 0., 0.375, 1., 1.);
    TPad* padRatio_bg = new TPad("padRatio_bg", "", 0., 0.,    1., 0.375);
    TPad* padSpecSigTotal  = new TPad("padSpecSigTotal",  "", 0., 0.375, 1., 1.);
    TPad* padRatioSigTotal = new TPad("padRatioSigTotal", "", 0., 0.,    1., 0.375);
    TPad* padSpecBkgTotal  = new TPad("padSpecBkgTotal",  "", 0., 0.375, 1., 1.);
    TPad* padRatioBkgTotal = new TPad("padRatioBkgTotal", "", 0., 0.,    1., 0.375);
    TLegend* leg = new TLegend(0.5, 0.5, 0.88, 0.88);
    map<int, TLegendEntry*> legEntries;

    // Draw nominal ------------
    // signal
    TH1* nomSig = (TH1*)nominal.first->Clone("nomSig");
    for (size_t bin = 1; bin <= size_t(nomSig->GetNbinsX()); ++bin)
      nomSig->SetBinContent(bin, nomSig->GetBinContent(bin)/nomSig->GetBinWidth(bin));
    nomSig->GetYaxis()->SetTitle("Events / GeV");
    nomSig->SetLineStyle(kSolid);
    nomSig->SetLineColor(kBlack);
    nomSig->SetLineWidth(2);
    padSpecSigTotal->cd();
    nomSig->Draw("hist");
    // background
    TH1* nomBkg  = (TH1*)nominal.second->Clone("nomBkg");
    for (size_t bin = 1; bin <= size_t(nomBkg->GetNbinsX()); ++bin)
      nomBkg->SetBinContent(bin, nomBkg->GetBinContent(bin)/nomBkg->GetBinWidth(bin));
    nomBkg->GetYaxis()->SetTitle("Events / GeV");
    nomBkg->SetLineStyle(kSolid);
    nomBkg->SetLineColor(kBlack);
    nomBkg->SetLineWidth(2);
    padSpecBkgTotal->cd();
    nomBkg->Draw("hist");
    // ratio
    TH1* nomRatio = (TH1*)nomSig->Clone("nomRatio");
    nomRatio->GetYaxis()->SetTitle("Shifted / Nominal");
    for (Int_t bin = 1; bin <= nomRatio->GetNbinsX(); ++bin)
      nomRatio->SetBinContent(bin, 1.);
    padRatioSigTotal->cd();
    nomRatio->Draw("hist");
    padRatioBkgTotal->cd();
    nomRatio->Draw("hist");
    legEntries[0] = new TLegendEntry(nomSig, "Nominal", "l");

    // Draw shifts ------------
    int it = 0;
    for (map<int, SystMakerShiftSigma*>::const_iterator sigma = fSigmas.begin(); sigma != fSigmas.end(); ++sigma, ++it) {
      if (sigma->second->fSigma == 0)
        continue;
      // Total ------------
      // signal
      TH1* sigmaSig = (TH1*)sigma->second->fShiftedSpectrum.first->Clone();
      for (size_t bin = 1; bin <= size_t(sigmaSig->GetNbinsX()); ++bin)
        sigmaSig->SetBinContent(bin, sigmaSig->GetBinContent(bin)/sigmaSig->GetBinWidth(bin));
      sigmaSig->GetYaxis()->SetTitle("Events / GeV");
      sigmaSig->SetLineStyle(kDashed);
      sigmaSig->SetLineColor(abs(sigma->second->fSigma)+1);
      sigmaSig->SetLineWidth(2);
      padSpecSigTotal->cd();
      sigmaSig->DrawClone("hist same");
      TH1* sigmaRatioSig = (TH1*)sigmaSig->Clone("sigmaRatioSig");
      sigmaRatioSig->Divide(nomSig);
      padRatioSigTotal->cd();
      sigmaRatioSig->DrawClone("hist same");
      // background
      TH1* sigmaBkg  = (TH1*)sigma->second->fShiftedSpectrum.second->Clone();
      for (size_t bin = 1; bin <= size_t(sigmaBkg->GetNbinsX()); ++bin)
        sigmaBkg->SetBinContent(bin, sigmaBkg->GetBinContent(bin)/sigmaBkg->GetBinWidth(bin));
      sigmaBkg->GetYaxis()->SetTitle("Events / GeV");
      sigmaBkg->SetLineStyle(kDashed);
      sigmaBkg->SetLineColor(abs(sigma->second->fSigma)+1);
      sigmaBkg->SetLineWidth(2);
      padSpecBkgTotal->cd();
      sigmaBkg->DrawClone("hist same");
      TH1* sigmaRatioBkg = (TH1*)sigmaBkg->Clone("sigmaRatioBkg");
      sigmaRatioBkg->Divide(nomBkg);
      padRatioBkgTotal->cd();
      sigmaRatioBkg->DrawClone("hist same");
      legEntries[abs(sigma->second->fSigma)]
        = new TLegendEntry(sigmaSig->Clone(), Form("+/-%d#sigma", sigma->second->fSigma), "l");
      // Individual ------------
      // signal
      padSpecSig->cd();
      padSpecSig->Clear();
      nomSig->Draw("hist");
      sigmaSig->SetLineColor(kRed);
      sigmaSig->Draw("hist same");
      padRatioSig->cd();
      padRatioSig->Clear();
      nomRatio->Draw("hist");
      sigmaRatioSig->SetLineColor(kRed);
      sigmaRatioSig->Draw("hist same");
      canv->cd();
      canv->Clear();
      padSpecSig->DrawClone();
      padRatioSig->DrawClone();
      outDir->cd(dirSigShift_sigmas.c_str());
      canv->Write(Form("%s%dSigma", fName.c_str(), sigma->second->fSigma));
      // background
      padSpec_bg->cd();
      padSpec_bg->Clear();
      nomBkg->Draw("hist");
      sigmaBkg->SetLineColor(kRed);
      sigmaBkg->Draw("hist same");
      padRatio_bg->cd();
      padRatio_bg->Clear();
      nomRatio->Draw("hist");
      sigmaRatioBkg->SetLineColor(kRed);
      sigmaRatioBkg->Draw("hist same");
      canv->cd();
      canv->Clear();
      padSpec_bg->DrawClone();
      padRatio_bg->DrawClone();
      outDir->cd(dirBkgShiftSigmas.c_str());
      canv->Write(Form("%s%dSigma", fName.c_str(), sigma->second->fSigma));

    }

    // Draw legend ------------
    for (const auto& legEntry : legEntries)
      leg->AddEntry(legEntry.second->GetObject(),
        legEntry.second->GetLabel(),
        legEntry.second->GetOption());
    padSpecSigTotal->cd();
    leg->Draw();
    padSpecBkgTotal->cd();
    leg->Draw();

    // Draw all together ------------
    // signal
    canv->cd();
    canv->Clear();
    padSpecSigTotal->Draw();
    padRatioSigTotal->Draw();
    outDir->cd(dirSigShift.c_str());
    canv->Write(Form("%sSigAllSigmas", fName.c_str()));
    // background
    canv->cd();
    canv->Clear();
    padSpecBkgTotal->Draw();
    padRatioBkgTotal->Draw();
    outDir->cd(dirBkgShift.c_str());
    canv->Write(Form("%sBgAllSigmas", fName.c_str()));

    // clean up
    // need to figure out what needs to be deleted
    delete canv;

  } // function SystMakerShift::DrawShift

  // --------------------------------------------------------------------------
  string SystMakerShift::WriteTable(pair<TH1*, TH1*> nominal, string name) {

    ostringstream table;

    table << "\\newcommand{\\" << name << "}\n";
    table << "{\\centerline{\n";
    table << "\\begin{tabular}{ l | r r r | r r }\n";
    table << " & Nom.&(+)&(-)&\\%(+)&\\%(-) \\\\ \\hline\n";

    // We write +1 and -1 sigmas to a table. First check they both exist!

    vector<int> sigmas = this->GetSigmas();
    if (find(sigmas.begin(), sigmas.end(), 1) == sigmas.end())
      assert(false and "+1 shift required to call WriteTable()!");
    bool useDown = find(sigmas.begin(), sigmas.end(), -1) != sigmas.end();

    // Get nominal event counts for signal and background
    double sigNom = 0, bkgNom = 0;
    size_t nbins = nominal.first->GetNbinsX();
    for (size_t it = 1; it <= nbins; ++it) {
      sigNom += nominal.first->GetBinContent(it);
      bkgNom += nominal.second->GetBinContent(it);
    }

    // Now calculate the up and down totals
    double sigUp = sigNom, bkgUp = bkgNom, sigDown = 0, bkgDown = 0;
    if (useDown) {
      sigDown = sigNom;
      bkgDown = bkgNom;
    }

      for (size_t it = 1; it <= nbins; ++it) {
        sigUp += this->GetShiftedSpectrum(1).first->GetBinContent(it)
          - nominal.first->GetBinContent(it);
        sigDown += this->GetShiftedSpectrum(-1).first->GetBinContent(it)
          - nominal.first->GetBinContent(it);
        bkgUp += this->GetShiftedSpectrum(1).second->GetBinContent(it)
          - nominal.second->GetBinContent(it);
        if (useDown) {
          sigDown += this->GetShiftedSpectrum(-1).first->GetBinContent(it)
            - nominal.first->GetBinContent(it);
          bkgDown += this->GetShiftedSpectrum(-1).second->GetBinContent(it)
            - nominal.second->GetBinContent(it);
        }
      }

    // Write them to the table and return it
    table << "Signal & " << sigNom << " & " << sigUp << " & ";
    if (useDown) table << sigDown;
    else table << "N/A";
    table << " & " << 100*((sigUp/sigNom)-1) << " & ";
    if (useDown) table << 100*((sigDown/sigNom)-1);
    else table << "N/A";
    table << " \\\\\n";
    table << "Background & " << bkgNom << " & " << bkgUp << " & ";
    if (useDown) table << bkgDown;
    else table << "N/A";
    table << " & " << 100*((bkgUp/bkgNom)-1) << " & ";
    if (useDown) table << 100*((bkgDown/bkgNom)-1);
    else table << "N/A";
    table << " \\\\\n\\end{tabular}\n}}\n";

    return table.str();

  } // function SystMakerShift::WriteTable

  // -------------------------------------------------------------------------------------
  pair<TH1*, TH1*> SystMakerShift::GetShiftedSpectrum(int sigma) {
    if (!fSigmas.count(sigma)) {
      cout << "Error: attempting to retrieve " << sigma << " sigma shift from "
          << fName << "; does not exist.  Aborting..." << endl;
      abort();
    }
    return fSigmas[sigma]->fShiftedSpectrum;
  }

  // -------------------------------------------------------------------------------------
  void SystMakerShift::ProcessShift(pair<TH1*, TH1*>& nominal_spectrum,
               const Sample sample, IOscCalc* calc) {
    for (auto& sigma : fSigmas) {
      TH1* hShiftSig = GetSpectrum(sigma.second->fPrediction, true, sample, calc);
      TH1* hShiftBg = GetSpectrum(sigma.second->fPrediction, false, sample, calc);
      TH1* hShiftRatioSig = (TH1*)hShiftSig->Clone("shiftRatSig");
      hShiftRatioSig->Divide(nominal_spectrum.first);
      TH1* hShiftRatioBg = (TH1*)hShiftBg->Clone("shiftRatBkg");
      hShiftRatioBg->Divide(nominal_spectrum.second);
      double shiftSig = abs(hShiftSig->Integral() - nominal_spectrum.first->Integral()) /
        nominal_spectrum.first->Integral();
      double shiftBg = abs(hShiftBg->Integral() - nominal_spectrum.second->Integral()) /
        nominal_spectrum.second->Integral();
      sigma.second->fShiftedSpectrum = make_pair(hShiftSig, hShiftBg);
      sigma.second->fShiftedRatio = make_pair(hShiftRatioSig, hShiftRatioBg);
      sigma.second->fShift = make_pair(shiftSig, shiftBg);
    }
    return;
  }

  // -------------------------------------------------------------------------------------
  /// Base class implementation of MakePredictions, which just throws an error
  void SystMakerShift::MakePredictions(SystPredType predType,
    const HistAxis* axis, const HistAxis* numuAxis,
    const Cut* FDCut, const Cut* NDCut, const Cut* numuCut,
    const Var* weight, Loaders* loaders) {

    // This function has to be defined in the base class for <reasons> but the
    // user should never actually call it so it just throws an error

    assert(false and "Base class version of MakePredictions should never be called!");

  } // function SystMakerShift::MakePredictions

  // -------------------------------------------------------------------------------------
  /// SaveTo implementation for SystMakerShift
  void SystMakerShift::SaveTo(TDirectory* dir) const {

    TDirectory* tmp = gDirectory;

    dir->cd();
    TObjString("SystMakerShift").Write("type");

    TObjString(fName.c_str()).Write("name");
    TObjString(fLabel.c_str()).Write("label");

    const char* sigmaDir = "Sigmas";
    dir->mkdir(sigmaDir);
    for (const auto& sigma : fSigmas) {
      const char* sigmaSubdir = Form("%s/Sigma%d", sigmaDir, sigma.second->fSigma);
      dir->mkdir(sigmaSubdir);
      sigma.second->SaveTo(dir->GetDirectory(sigmaSubdir));
    }

    tmp->cd();

  } // function SystMakerShift::SaveTo

  // -------------------------------------------------------------------------------------
  /// LoadFrom implementation for SystMakerShift
  unique_ptr<SystMakerShift> SystMakerShift::LoadFrom(TDirectory* dir) {

    DontAddDirectory guard;

    TObjString* tag = (TObjString*)dir->Get("type");
    assert(tag);
    assert(tag->GetString() == "SystMakerShift");

    TObjString* name = (TObjString*)dir->Get("name");;
    TObjString* label = (TObjString*)dir->Get("label");;
    unique_ptr<SystMakerShift> systMakerShift
      = make_unique<SystMakerShift>(name->GetString().Data(), label->GetString().Data());

    map<int, SystMakerShiftSigma*> sigmas;
    TIter next(dir->GetDirectory("Sigmas")->GetListOfKeys());
    TKey* key;
    while ((key = (TKey*)next())) {
      SystMakerShiftSigma* sigma
        = SystMakerShiftSigma::LoadFrom(dir->GetDirectory(Form("Sigmas/%s", key->GetName())));
      sigmas[sigma->fSigma] = sigma;
    }
    systMakerShift->fSigmas = sigmas;

    return systMakerShift;

  } // function SystMakerShift::LoadFrom

  // -------------------------------------------------------------------------------------
  /// SystMakerLoaderShift constructor
  SystMakerLoaderShift::SystMakerLoaderShift(string name, string label) :
    SystMakerShift(name, label) {
  }

  // -------------------------------------------------------------------------------------
  /// Add sigma for loader shift with the same ND and FD loader
  void SystMakerLoaderShift::AddSigma(Loaders* loaders,
             int sigma) {

    this->AddSigma(loaders, loaders, sigma);

  } // function SystMakerLoaderShift::AddSigma

  // -------------------------------------------------------------------------------------
  /// Add sigma for loader shift with separate ND and FD loaders
  void SystMakerLoaderShift::AddSigma(Loaders* loadersND,
             Loaders* loaders_fd,
             int sigma) {

    SystMakerShiftSigma* shiftSigma = new SystMakerShiftSigma(fName, sigma);
    fSigmas[sigma] = shiftSigma;
    fLoaders[sigma] = make_pair(loadersND, loaders_fd);

  } // function SystMakerLoaderShift::AddSigma

  // -------------------------------------------------------------------------------------
  /// Add on/off for loader shift
  void SystMakerLoaderShift::AddOnOff(Loaders* loaders) {

    this->AddOnOff(loaders, loaders);

  } // function SystMakerLoaderShift::AddOnOff

  // -------------------------------------------------------------------------------------
  /// Add on/off for loader shift with separate ND and FD loaders
  void SystMakerLoaderShift::AddOnOff(Loaders* loadersND,
                                    Loaders* loaders_fd) {

    for (const int sigma : {-1,1}) {
      SystMakerShiftSigma* shiftSigma = new SystMakerShiftSigma(fName, sigma);
      fSigmas[sigma] = shiftSigma;
      fLoaders[sigma] = pair<Loaders*, Loaders*>(loadersND, loaders_fd);
    }

    // std::cout << "Calling CheckLoaders in AddOnOff" << std::endl;
    // CheckLoaders();

  } // function SystMakerLoaderShift::AddOnOff

  // -------------------------------------------------------------------------------------
  /// Make shifted prediction for systematic shift derived from alternate MC loader
  void SystMakerLoaderShift::MakePredictions(SystPredType predType,
    const HistAxis* axis, const HistAxis* numuAxis,
    const Cut* FDCut, const Cut* NDCut, const Cut* numuCut,
    const Var* weight, Loaders* loaders) {

    // CheckLoaders();

    for (auto& sigma : fSigmas) {

      sigma.second->fPrediction = GetPrediction(predType,
        axis, numuAxis,
        FDCut, NDCut, numuCut,
        &kNoShift, &kNoShift,
        weight,
        fLoaders[sigma.second->fSigma].second,
        fLoaders[sigma.second->fSigma].first);
    }

  } // function SystMakerLoaderShift::MakePredictions

  // -------------------------------------------------------------------------------------
  /// SystMakerWeightShift constructor
  SystMakerWeightShift::SystMakerWeightShift(string name, string label) :
    SystMakerShift(name, label) {
  }

  // -------------------------------------------------------------------------------------
  /// SystMakerWeightShift constructor
  SystMakerWeightShift::~SystMakerWeightShift() {
    for (auto& weight : fWeights)
      delete weight.second;
  }

  // -------------------------------------------------------------------------------------
  /// Add sigma for weight shift
  void SystMakerWeightShift::AddSigma(const Var* weight, int sigma) {

    SystMakerShiftSigma* shiftSigma = new SystMakerShiftSigma(fName, sigma);
    fSigmas[sigma] = shiftSigma;
    fWeights[sigma] = weight;
  } // function SystMakerWeightShift::AddSigma

  // -------------------------------------------------------------------------------------
  /// Add on/off for weight shift
  void SystMakerWeightShift::AddOnOff(const Var* weight) {

    for (const auto& sigma : {-1,1}) {
      SystMakerShiftSigma* shiftSigma = new SystMakerShiftSigma(fName, sigma);
      fSigmas[sigma] = shiftSigma;
      fWeights[sigma] = weight;
    }

  } // function SystMakerWeightShift::AddOnOff

  // -------------------------------------------------------------------------------------
  /// Make shifted prediction for systematic shift derived from CAFAna weight
  void SystMakerWeightShift::MakePredictions(SystPredType predType,
    const HistAxis* axis, const HistAxis* numuAxis,
    const Cut* FDCut, const Cut* NDCut, const Cut* numuCut,
    const Var* weight, Loaders* loaders) {

    for (auto& sigma : fSigmas)
      sigma.second->fPrediction = GetPrediction(predType, axis, numuAxis, FDCut,
      NDCut, numuCut, &kNoShift, &kNoShift, fWeights[sigma.second->fSigma], loaders);

  } // function SystMakerWeightShift::MakePrediction

  // -------------------------------------------------------------------------------------
  /// SystMakerSystShift constructor
  SystMakerSystShift::SystMakerSystShift(string name, string label) :
    SystMakerShift(name, label) {
  }

  // -------------------------------------------------------------------------------------
  /// Add sigma for syst shift
  void SystMakerSystShift::AddSigma(SystShifts* syst, int sigma) {
    SystMakerShiftSigma* shiftSigma = new SystMakerShiftSigma(fName, sigma);
    fSigmas[sigma] = shiftSigma;
    fSysts[sigma] = syst;
    return;
  }

  // -------------------------------------------------------------------------------------
  /// Make shifted prediction for systematic shift derived from SystShift
  void SystMakerSystShift::MakePredictions(SystPredType predType,
    const HistAxis* axis, const HistAxis* numuAxis,
    const Cut* FDCut, const Cut* NDCut, const Cut* numuCut,
    const Var* weight, Loaders* loaders) {

    for (auto& sigma : fSigmas) {
      sigma.second->fPrediction = GetPrediction(predType, axis, numuAxis, FDCut,
        NDCut, numuCut, &kNoShift, fSysts[sigma.second->fSigma], weight, loaders);
    }
  } // function SystMakerSystShift::MakePrediction

  // -------------------------------------------------------------------------------------
  /// SystMaker constructor
  SystMaker::SystMaker(string name, string label, string sampleName) :
    fName(name),
    fLabel(label),
    fSampleName(sampleName),
    fAxis(nullptr),
    fNuMuAxis(nullptr),
    fNDCut(nullptr),
    fFDCut(nullptr),
    fNuMuCut(nullptr),
    fLoaders(nullptr),
    fWeight(nullptr),
    fNominalPred(nullptr) {

  } // SystMaker constructor

  // -------------------------------------------------------------------------------------
  /// SystMaker destructor
  SystMaker::~SystMaker() {

  for (auto shift : fShifts) {
    if (shift.second) delete shift.second;
  }
  if (fNominalPred) delete fNominalPred;

  } // SystMaker destructor

  // -------------------------------------------------------------------------------------
  /// Get the shift for a given name
  SystMakerShift* SystMaker::GetShift(string shift_name) {
    if (!fShifts.count(shift_name)) {
      cout << "Error: requested shift " << shift_name << " does not exist.  Aborting..." << endl;
      abort();
    }
    return fShifts[shift_name];

  } // function SystMaker::GetShift

  // -------------------------------------------------------------------------------------
  /// Get the total shifts for a given sigma (sig+bkg)
  pair<double, double> SystMaker::GetTotalShift(int sigma) {
    if (!fTotalShifts.count(sigma)) {
      cout << "Error: sigma " << sigma << " does not exist.  Aborting..." << endl;
      abort();
    }
    return fTotalShifts[sigma];

  } // function SystMakerShift::GetTotalShift

  // -------------------------------------------------------------------------------------
  /// Get the total shift for a given sigma and channel selection (sig/bkg)
  double SystMaker::GetTotalShift(int sigma, bool signal) {

    pair<double, double> totalShift = this->GetTotalShift(sigma);
    return signal ? totalShift.first : totalShift.second;

  } // function SystMakerShift::GetTotalShift

  // -------------------------------------------------------------------------------------
  /// Get the sigmas for a given systematic shift
  vector<int> SystMaker::GetSigmas() {

    map<int, vector<string> > shifts_sigmas;
    for (const auto& shift : fShifts) {
      vector<int> shiftSigmas = shift.second->GetSigmas();
      for (const auto& shiftSigma : shiftSigmas)
        shifts_sigmas[shiftSigma].push_back(shift.second->GetName());
    }
    vector<int> sigmas;
    for (const auto& shiftSigmas : shifts_sigmas)
      if (shiftSigmas.second.size() == fShifts.size())
        sigmas.push_back(shiftSigmas.first);

    return sigmas;

  } // function SystMaker::GetSigmas

  // -------------------------------------------------------------------------------------
  /// Print the total shifts for each shift owned by this syst maker
  void SystMaker::PrintSyst(ostream& os) {

    os << endl << endl << showpos
       << "Systematic " << fLabel << " has total shift:" << endl
       << "  -1: " << fTotalShifts[-1].first << " (signal), " << fTotalShifts[-1].second << " (background)" << endl
       << "  +1: " << fTotalShifts[+1].first << " (signal), " << fTotalShifts[+1].second << " (background)" << endl;
    os << endl << "Individual shifts:" << endl;
    for (const auto& shift : fShifts)
      shift.second->PrintShift(os);

  } // function SystMaker::PrintSyst

  // -------------------------------------------------------------------------------------
  /// Add SystMakerShift to this syst maker
  void SystMaker::AddShift(SystMakerShift* shift) {

    fShifts[shift->GetName()] = shift;

  } // function SystMaker::AddShift

  // -------------------------------------------------------------------------------------
  /// Make plots of all the shifted spectra and ratios for this systmaker
  TCanvas* SystMaker::DrawSig() {

    // Nominal
    TH1* nomSig = (TH1*)fSpectra[0].first->Clone();
    for (size_t bin = 1; bin <= size_t(nomSig->GetNbinsX()); ++bin)
      nomSig->SetBinContent(bin, nomSig->GetBinContent(bin)/nomSig->GetBinWidth(bin));
    nomSig->GetYaxis()->SetTitle("Events / GeV");
    nomSig->SetLineColor(kBlack);
    nomSig->SetLineStyle(kSolid);
    nomSig->SetLineWidth(2);
    TH1* nomRatio = (TH1*)nomSig->Clone();
    nomRatio->GetYaxis()->SetTitle("Shifted / Nominal");
    for (size_t bin = 0; bin <= (size_t)nomRatio->GetNbinsX(); ++bin)
      nomRatio->SetBinContent(bin, 1.);

    // Draw total shifts -- all sigmas
    TCanvas* canv = new TCanvas("canvSig", "", 800, 800);
    TPad* padSpec  = new TPad("padSpecSig",  "", 0., 0.375, 1., 1.   );
    TPad* padRatio = new TPad("padRatioSig", "", 0., 0.,    1., 0.375);
    map<int, TLegendEntry*> legEntrySigma;
    padSpec->cd();
    nomSig->Draw("hist");
    padRatio->cd();
    nomRatio->Draw("hist");
    legEntrySigma[0] = new TLegendEntry(nomSig, "Nominal", "l");
    double ymaxRatio = 0;

    // Loop over each sigma
    for (const auto& sigma : this->GetSigmas()) {
      if (sigma == 0)
        continue;

      // Draw the spectrum
      TH1* sigmaSig = (TH1*)fSpectra[sigma].first->Clone();
      for (size_t bin = 1; bin <= size_t(sigmaSig->GetNbinsX()); ++bin)
        sigmaSig->SetBinContent(bin, sigmaSig->GetBinContent(bin)/sigmaSig->GetBinWidth(bin));
      sigmaSig->GetYaxis()->SetTitle("Events / GeV");
      double hmax = nomSig->GetMaximum() > sigmaSig->GetMaximum() ?
        nomSig->GetMaximum() : sigmaSig->GetMaximum();
      nomSig->SetMaximum(1.1*hmax);
      sigmaSig->SetLineColor(abs(sigma)+1);
      sigmaSig->SetLineStyle(kDashed);
      sigmaSig->SetLineWidth(2);
      padSpec->cd();
      sigmaSig->Draw("hist same");
      if (!legEntrySigma.count(abs(sigma)))
        legEntrySigma[abs(sigma)] = new TLegendEntry(sigmaSig, Form("+/-%d #sigma", abs(sigma)), "l");

      // Draw the ratio
      TH1* ratioSig = (TH1*)sigmaSig->Clone();
      ratioSig->Divide(nomSig);
      ratioSig->SetBinContent(0, 1);
      for (size_t i = 1; i <= (size_t)ratioSig->GetNbinsX(); ++i) {
        double val = ratioSig->GetBinContent(i);
        if (val == 0) {
          ratioSig->SetBinContent(i, 1);
          val = 1;
        }
        if (fabs(1-val) > ymaxRatio)
          ymaxRatio = fabs(1-val);
      }
      padRatio->cd();
      ratioSig->Draw("hist same");
    }
    TLegend* legSigma = new TLegend(0.5, 0.5, 0.88, 0.88);
    for (const auto& legEntry : legEntrySigma)
      legSigma->AddEntry(legEntry.second->GetObject(),
        legEntry.second->GetLabel(),
        legEntry.second->GetOption());
    padSpec->cd();
    legSigma->Draw("same");

    nomRatio->SetMinimum(1-(1.1*ymaxRatio));
    nomRatio->SetMaximum(1+(1.1*ymaxRatio));

    canv->cd();
    canv->Clear();
    padSpec->Draw();
    padRatio->Draw();

    return canv;

  } // function SystMaker::DrawSig

  // -------------------------------------------------------------------------------------
  /// Make plots of all the shifted spectra and ratios for this systmaker
  TCanvas* SystMaker::DrawBkg() {

    // Nominal
    TH1* nomBkg = (TH1*)fSpectra[0].second->Clone();
    for (size_t bin = 1; bin <= size_t(nomBkg->GetNbinsX()); ++bin)
      nomBkg->SetBinContent(bin, nomBkg->GetBinContent(bin)/nomBkg->GetBinWidth(bin));
    nomBkg->GetYaxis()->SetTitle("Events / GeV");
    nomBkg->SetLineColor(kBlack);
    nomBkg->SetLineStyle(kSolid);
    nomBkg->SetLineWidth(2);
    TH1* nomRatio = (TH1*)nomBkg->Clone();
    nomRatio->GetYaxis()->SetTitle("Shifted / Nominal");
    for (size_t bin = 0; bin <= (size_t)nomRatio->GetNbinsX(); ++bin)
      nomRatio->SetBinContent(bin, 1.);

    // Draw total shifts -- all sigmas
    TCanvas* canv = new TCanvas("canvBkg", "", 800, 800);
    TPad* padSpec  = new TPad("padSpecBkg",  "", 0., 0.375, 1., 1.   );
    TPad* padRatio = new TPad("padRatio_bkg", "", 0., 0.,    1., 0.375);
    map<int, TLegendEntry*> legEntrySigma;
    padSpec->cd();
    nomBkg->Draw("hist");
    padRatio->cd();
    nomRatio->Draw("hist");
    legEntrySigma[0] = new TLegendEntry(nomBkg, "Nominal", "l");
    double ymaxRatio = 0;

    // Loop over each sigma
    for (const auto& sigma : this->GetSigmas()) {
      if (sigma == 0)
        continue;

      // Draw the spectrum
      TH1* sigmaBkg = (TH1*)fSpectra[sigma].second->Clone();
      for (size_t bin = 1; bin <= size_t(sigmaBkg->GetNbinsX()); ++bin)
        sigmaBkg->SetBinContent(bin, sigmaBkg->GetBinContent(bin)/sigmaBkg->GetBinWidth(bin));
      sigmaBkg->GetYaxis()->SetTitle("Events / GeV");
      double hmax = nomBkg->GetMaximum() > sigmaBkg->GetMaximum() ?
        nomBkg->GetMaximum() : sigmaBkg->GetMaximum();
      nomBkg->SetMaximum(1.1*hmax);
      sigmaBkg->SetLineColor(abs(sigma)+1);
      sigmaBkg->SetLineStyle(kDashed);
      sigmaBkg->SetLineWidth(2);
      padSpec->cd();
      sigmaBkg->Draw("hist same");
      if (!legEntrySigma.count(abs(sigma)))
        legEntrySigma[abs(sigma)] = new TLegendEntry(sigmaBkg, Form("+/-%d #sigma", abs(sigma)), "l");

      // Draw the ratio
      TH1* ratioBkg = (TH1*)sigmaBkg->Clone();
      ratioBkg->Divide(nomBkg);
      ratioBkg->SetBinContent(0, 1);
      for (size_t i = 1; i <= (size_t)ratioBkg->GetNbinsX(); ++i) {
        double val = ratioBkg->GetBinContent(i);
        if (val == 0) {
          ratioBkg->SetBinContent(i, 1);
          val = 1;
        }
        if (fabs(1-val) > ymaxRatio)
          ymaxRatio = fabs(1-val);
      }
      padRatio->cd();
      ratioBkg->Draw("hist same");
    }
    TLegend* legSigma = new TLegend(0.5, 0.5, 0.88, 0.88);
    for (const auto& legEntry : legEntrySigma)
      legSigma->AddEntry(legEntry.second->GetObject(),
        legEntry.second->GetLabel(),
        legEntry.second->GetOption());
    padSpec->cd();
    legSigma->Draw("same");

    nomRatio->SetMinimum(1-(1.1*ymaxRatio));
    nomRatio->SetMaximum(1+(1.1*ymaxRatio));

    canv->cd();
    canv->Clear();
    padSpec->Draw();
    padRatio->Draw();

    return canv;

  } // function SystMaker::DrawBkg

  // -------------------------------------------------------------------------------------
  /// Make plots of all the shifted spectra and ratios for this systmaker
  void SystMaker::DrawSyst(TDirectory* outDir) {

    string dirTop = fName;
    string dirSignal = dirTop+"/Signal/";
    string dirBackground = dirTop+"/Background/";
    string dirSigShifts = dirSignal+"/Shifts/";
    string dirBkgShifts = dirBackground+"/Shifts/";
    outDir->mkdir(dirTop.c_str());
    outDir->mkdir(dirSignal.c_str());
    outDir->mkdir(dirBackground.c_str());
    outDir->mkdir(dirSigShifts.c_str());
    outDir->mkdir(dirBkgShifts.c_str());

    TCanvas* canvSig = DrawSig();
    TCanvas* canvBkg = DrawBkg();
    outDir->cd(dirSignal.c_str());
    canvSig->Write("TotalShiftSig");
    outDir->cd(dirBackground.c_str());
    canvBkg->Write("TotalShiftBg");

    // Draw total shifts -- +1/-1 sigma

    // Draw individual shifts
    for (const auto& shift : fShifts)
      shift.second->DrawShift(outDir->GetDirectory(dirTop.c_str()), fSpectra[0]);

    // clean up
    delete canvSig;
    delete canvBkg;

  } // function SystMaker::DrawSyst

  // --------------------------------------------------------------------------
  /// Write a table with up & down shifts for this systematic
  string SystMaker::WriteTable(string name) {

    ostringstream table;

    table << "\\newcommand{\\" << name << "}\n";
    table << "{\\centerline{\n";
    table << "\\begin{tabular}{ l | r r r | r r }\n";
    table << " & Nom.&(+)&(-)&\\%(+)&\\%(-) \\\\ \\hline\n";

    // We write +1 and -1 sigmas to a table. First check they both exist!
    vector<int> sigmas = this->GetSigmas();
    if (find(sigmas.begin(), sigmas.end(), 1) == sigmas.end())
      assert(false and "+1 sigma shift required to call WriteTable()!");
    bool useDown = find(sigmas.begin(), sigmas.end(), -1) != sigmas.end();

    // Get nominal event counts for signal and background
    double sigNom = 0, bkgNom = 0;
    size_t nbins = fSpectra[0].first->GetNbinsX();
    for (size_t it = 1; it <= nbins; ++it) {
      sigNom += fSpectra[0].first->GetBinContent(it);
      bkgNom += fSpectra[0].second->GetBinContent(it);
    }

    // Now calculate the up and down totals
    double sigUp = sigNom, bkgUp = bkgNom, sigDown = 0, bkgDown = 0;
    if (useDown) {
      sigDown = sigNom;
      bkgDown = bkgNom;
    }

    for (auto& shift : fShifts) {
      for (size_t it = 1; it <= nbins; ++it) {
        sigUp += shift.second->GetShiftedSpectrum(1).first->GetBinContent(it)
          - fSpectra[0].first->GetBinContent(it);
        bkgUp += shift.second->GetShiftedSpectrum(1).second->GetBinContent(it)
          - fSpectra[0].second->GetBinContent(it);
        if (useDown) {
          sigDown += shift.second->GetShiftedSpectrum(-1).first->GetBinContent(it)
            - fSpectra[0].first->GetBinContent(it);
          bkgDown += shift.second->GetShiftedSpectrum(-1).second->GetBinContent(it)
            - fSpectra[0].second->GetBinContent(it);
        }
      }
    }

    // Write them to the table and return it
    table << "Signal & " << sigNom << " & " << sigUp << " & ";
    if (useDown) table << sigDown;
    else table << "N/A";
    table << " & " << 100*((sigUp/sigNom)-1) << " & ";
    if (useDown) table << 100*((sigDown/sigNom)-1);
    else table << "N/A";
    table << " \\\\\n";
    table << "Background & " << bkgNom << " & " << bkgUp << " & ";
    if (useDown) table << bkgDown;
    else table << "N/A";
    table << " & " << 100*((bkgUp/bkgNom)-1) << " & ";
    if (useDown) table << 100*((bkgDown/bkgNom)-1);
    else table << "N/A";
    table << " \\\\\n\\end{tabular}\n}}\n";

    return table.str();

  } // function SystMaker::WriteTable

  // -------------------------------------------------------------------------------------
  /// Make predictions for all the SystMakerShift objects owned by this syst maker
  void SystMaker::MakePredictions(SystPredType predType) {

    fNominalPred = GetPrediction(predType,
               fAxis, fNuMuAxis,
               fFDCut, fNDCut, fNuMuCut,
               &kNoShift, &kNoShift,
               fWeight, fLoaders);

    for (const auto& shift : fShifts) {
      shift.second->MakePredictions(predType,
            fAxis, fNuMuAxis,
            fFDCut, fNDCut, fNuMuCut,
            fWeight, fLoaders);
    }
  } // function SystMaker::MakePredictions

  // -------------------------------------------------------------------------------------
  void SystMaker::ProcessSyst(const Sample sample, IOscCalc* calc) {

    // Get nominal spectrum
    fSpectra[0] = make_pair(GetSpectrum(fNominalPred, true, sample, calc),
      GetSpectrum(fNominalPred, false, sample, calc));

    // Process all shifts
    for (const auto& shift : fShifts) {
      shift.second->ProcessShift(fSpectra[0], sample, calc);
      SystMakerShift* systShift = static_cast<SystMakerShift*>(shift.second);
      fShifts[shift.first] = systShift;
    }

    // Get total +1/-1 sigma systematic shift -- required
    // Also get total sigma shift for all sigmas in this syst
    vector<int> sigmas = this->GetSigmas();
    Int_t nbins = fSpectra[0].first->GetNbinsX();
    for (const auto& sigma : sigmas) {
      vector<double> evtDiffSig(nbins, 0.), evtDiffBg(nbins, 0.);
      TH1* shiftSig = (TH1*)fSpectra[0].first->Clone("shiftSig");
      TH1* shiftBkg = (TH1*)fSpectra[0].second->Clone("shiftBkg");
      double totShiftSig = 0., totShiftBkg = 0.;
      for (Int_t bin = 1; bin <= nbins; ++bin) {
        double valSig = 0., valBkg = 0.;
        // If there's more than one shift, we should add them in quadrature
        if (fShifts.size() > 1) {
          for (auto& shift : fShifts) {
            double binShiftSig = shift.second->GetShiftedSpectrum(sigma).first->GetBinContent(bin) -
              fSpectra[0].first->GetBinContent(bin);
            double binShiftBkg = shift.second->GetShiftedSpectrum(sigma).second->GetBinContent(bin) -
              fSpectra[0].second->GetBinContent(bin);
            valSig += pow(binShiftSig,2);
            valBkg += pow(binShiftBkg,2);
          }
          double sign = sigma > 0 ? 1.0 : -1.0;
          shiftSig->SetBinContent(bin, fSpectra[0].first->GetBinContent(bin) + (sign * sqrt(valSig)));
          shiftBkg->SetBinContent(bin, fSpectra[0].second->GetBinContent(bin) + (sign * sqrt(valBkg)));
        }
        // If there's only one shift, just use the raw value
        else {
          valSig = fShifts.begin()->second->GetShiftedSpectrum(sigma).first->GetBinContent(bin)
            - fSpectra[0].first->GetBinContent(bin);
          valBkg = fShifts.begin()->second->GetShiftedSpectrum(sigma).second->GetBinContent(bin)
            - fSpectra[0].second->GetBinContent(bin);
          shiftSig->SetBinContent(bin, fShifts.begin()->second->GetShiftedSpectrum(sigma).first->GetBinContent(bin));
          shiftBkg->SetBinContent(bin, fShifts.begin()->second->GetShiftedSpectrum(sigma).second->GetBinContent(bin));
        }
        totShiftSig += valSig;
        totShiftBkg += valBkg;
      }
      fTotalShifts[sigma] = make_pair(totShiftSig, totShiftBkg);
      fSpectra[sigma] = make_pair(shiftSig, shiftBkg);
    }

  } // function SystMaker::ProcessSyst

  // -------------------------------------------------------------------------------------
  void SystMaker::SetAxes(const HistAxis* axis, const HistAxis* numuAxis) {
    fAxis = axis;
    fNuMuAxis = axis;
  }

  // -------------------------------------------------------------------------------------
  void SystMaker::SetCuts(const Cut* FDCut, const Cut* NDCut, const Cut* numuCut) {
    fNDCut = NDCut;
    fFDCut = FDCut;
    fNuMuCut = numuCut;
  }

  // -------------------------------------------------------------------------------------
  void SystMaker::SetNominalLoaders(Loaders* loaders) {
    fLoaders = loaders;
  }

  // -------------------------------------------------------------------------------------
  void SystMaker::SetNominalWeight(const Var* weight) {
    fWeight = weight;
  }

  // -------------------------------------------------------------------------------------
  /// Function to make ISyst for this SystMaker
  NuISyst* SystMaker::MakeISyst(covmx::Sample sample) {

    if (fSpectra.size() == 0) {
      cout << "There are no spectra owned by this SystMaker!"
                << " Did you forget to call Go()?" << endl;
      return nullptr;
    }

    // Make ratios from shifted spectra
    map<int, pair<TH1*, TH1*>> shifts;
    for (const auto& sigma : fSpectra) {
      TH1* sigRatio = (TH1*)sigma.second.first->Clone(Form("SigRatioSigma%d", sigma.first));
      TH1* bkgRatio = (TH1*)sigma.second.second->Clone(Form("BgRatioSigma%d", sigma.first));
      //sigRatio->Divide(fSpectra.at(0).first);
      //bkgRatio->Divide(fSpectra.at(0).second);
      // BEGIN DEBUG CODE
      for (int i = 1; i <= sigRatio->GetNbinsX(); ++i) {
        double eps = 1e-3;
        double num = sigRatio->GetBinContent(i);
        double den = fSpectra.at(0).first->GetBinContent(i);
        double val = (num+eps)/(den+eps);
        // if (val > 2) {
        //   cout << "Signal bin " << i << ", dividing thru "
        //        << num << " by " << den << " to get " << val
        //        << " - with epsilon it's " << valEps << endl;
        // }
        sigRatio->SetBinContent(i, val);
        num = bkgRatio->GetBinContent(i);
        den = fSpectra.at(0).second->GetBinContent(i);
        val = (num+eps)/(den+eps);
        bkgRatio->SetBinContent(i, val);
        // if (val > 2) {
        //   cout << "Background bin " << i << ", dividing thru "
        //        << num << " by " << den << " to get " << val
        //        << " - with epsilon it's " << valEps << endl;
        // }
      }
      // END DEBUG CODE
      shifts[sigma.first] = make_pair(sigRatio, bkgRatio);
    }

    // Return new ISyst
    string shortName = fSampleName.empty() ? fName : fSampleName + "_" + fName;
    return new NuISyst(shortName, fLabel, fSampleName, shifts);

  } // function SystMaker::MakeISyst

  // -------------------------------------------------------------------------------------
  /// SaveTo implementation for SystMaker
  void SystMaker::SaveTo(TDirectory* dir) const {

    TDirectory* tmp = gDirectory;

    dir->cd();
    TObjString("SystMaker").Write("type");

    TObjString(fName.c_str()).Write("name");
    TObjString(fLabel.c_str()).Write("label");
    TObjString(fSampleName.c_str()).Write("samplename");

    fNominalPred->SaveTo(dir, "nompred");

    const char* shifts_dir = "Shifts";
    dir->mkdir("Shifts");
    for (const auto& shift : fShifts) {
      const char* shiftSubdir = Form("%s/%s", shifts_dir, shift.first.c_str());
      dir->mkdir(shiftSubdir);
      shift.second->SaveTo(dir->GetDirectory(shiftSubdir));
    }

    tmp->cd();

  } // function SystMaker::SaveTo

  // -------------------------------------------------------------------------------------
  /// LoadFrom implementation for SystMaker
  unique_ptr<SystMaker> SystMaker::LoadFrom(TDirectory* dir) {

    DontAddDirectory guard;

    TObjString* tag = (TObjString*)dir->Get("type");
    assert(tag);
    assert(tag->GetString() == "SystMaker");

    TObjString* name = (TObjString*)dir->Get("name");
    assert(name);

    TObjString* label = (TObjString*)dir->Get("label");
    assert(label);

    TObjString* sampleName = (TObjString*)dir->Get("samplename");
    assert(sampleName);

    unique_ptr<SystMaker> systMaker
      = make_unique<SystMaker>(name->GetString().Data(),
        label->GetString().Data(), sampleName->GetString().Data());

    // Load prediction
    tag = (TObjString*)dir->Get("nompred/type");
    assert(tag);
    if (tag->GetString() == "NDPredictionSterile")
      systMaker->fNominalPred = NDPredictionSterile::LoadFrom(dir, "nompred").release();
    else if (tag->GetString() == "FDPredictionSterile")
      systMaker->fNominalPred = FDPredictionSterile::LoadFrom(dir, "nompred").release();
    else if (tag->GetString() == "PredictionSterile")
      systMaker->fNominalPred = PredictionSterile::LoadFrom(dir, "nompred").release();
    else assert(false && "Prediction type not recognised!");

    // Load shift makers
    map<string, SystMakerShift*> shifts;

    TIter shift_next(dir->GetDirectory("Shifts")->GetListOfKeys());
    TKey* shift_key;
    while ((shift_key = (TKey*)shift_next())) {
      SystMakerShift* shift = SystMakerShift::LoadFrom(dir->GetDirectory(Form("Shifts/%s",
        shift_key->GetName()))).release();
      shifts[shift->GetName()] = shift;
    }
    systMaker->fShifts = shifts;

    return move(systMaker);

  } // function SystMaker::LoadFrom

  // -------------------------------------------------------------------------------------
  /// SystematicsMaker constructor
  SystematicsMaker::SystematicsMaker(string name, SystPredType predType) :
    fName(name),
    fPredType(predType),
    fAxis(nullptr),
    fNuMuAxis(nullptr),
    fNDCut(nullptr),
    fFDCut(nullptr),
    fNuMuCut(nullptr),
    fLoaders(nullptr),
    fLightLoaders(nullptr),
    fWeight(nullptr),
    fSample(NullSample)
  {
    vector<SystPredType> supportedPredTypes = {SystPredType::kExtrap, SystPredType::kFD, SystPredType::kND};
    if (find(supportedPredTypes.begin(), supportedPredTypes.end(), predType) == supportedPredTypes.end()) {
      cout << "SystematicsMaker error: Unsupported prediction type.  Aborting..." << endl;
      abort();
    }
  } // SystematicsMaker constructor

  // -------------------------------------------------------------------------------------
  /// SystematicsMaker destructor
  SystematicsMaker::~SystematicsMaker() {
    for (auto systMaker : fSystMakers) {
      if (systMaker.second) delete systMaker.second;
    }
  } // function SystematicsMaker::~SystematicsMaker

  // -------------------------------------------------------------------------------------
  /// Get the names of all attached systematics
  vector<string> SystematicsMaker::GetSystematicsNames() {

    vector<string> systNames;
    for (const auto& syst : fSystMakers)
      systNames.push_back(syst.first);
    return systNames;

  } // function SystematicsMaker::GetSystematicsNames

  // -------------------------------------------------------------------------------------
  /// Make ISyst objects for each attached systematic and return them
  vector<NuISyst*> SystematicsMaker::MakeISysts() {
    
    vector<NuISyst*> systs;
    for (const auto& maker : fSystMakers) {
      NuISyst* syst = maker.second->MakeISyst(fSample);
      systs.push_back(syst);
    }
    return systs;

  } // function SystematicsMaker::MakeISysts

  // -------------------------------------------------------------------------------------
  void SystematicsMaker::AddSystematic(SystMaker* syst) {

    syst->SetAxes(fAxis, fNuMuAxis);
    syst->SetCuts(fFDCut, fNDCut, fNuMuCut);
    syst->SetNominalWeight(fWeight);
    if (syst->GetName().find("LightLevel") != string::npos)
      syst->SetNominalLoaders(fLightLoaders);
    else
      syst->SetNominalLoaders(fLoaders);
    syst->MakePredictions(fPredType);
    fSystMakers[syst->GetName()] = syst;
  }

  // -------------------------------------------------------------------------------------
  void SystematicsMaker::PrintSysts(ostream& os) {
    for (const auto& syst : fSystMakers)
      syst.second->PrintSyst(os);
    return;
  }

  // -------------------------------------------------------------------------------------
  void SystematicsMaker::DrawSysts(TDirectory* outDir) {
    for (const auto& syst : fSystMakers)
      syst.second->DrawSyst(outDir);
    return;
  }

  // -------------------------------------------------------------------------------------
  void SystematicsMaker::Go(Sample sample, IOscCalc* calc) {

    fSample = sample;
    for (const auto& syst : fSystMakers) {
      syst.second->ProcessSyst(fSample, calc);
      fSystMakers[syst.first] = syst.second;
    }
  }

  // -------------------------------------------------------------------------------------
  void SystematicsMaker::SetAxes(const HistAxis* axis, const HistAxis* numuAxis) {
    fAxis = axis;
    fNuMuAxis = numuAxis;
  }

  // -------------------------------------------------------------------------------------
  void SystematicsMaker::SetCuts(const Cut* FDCut, const Cut* NDCut,
    const Cut* numuCut) {
    fNDCut = NDCut;
    fFDCut = FDCut;
    fNuMuCut = numuCut;
  }

  // -------------------------------------------------------------------------------------
  void SystematicsMaker::SetNominalLoaders(Loaders* loaders, Loaders* loaders_light) {
    fLoaders = loaders;
    fLightLoaders = loaders_light;
  }

  // -------------------------------------------------------------------------------------
  void SystematicsMaker::SetNominalWeight(const Var* weight) {
    fWeight = weight;
  }

  // -------------------------------------------------------------------------------------
  /// function to add another SystematicsMaker to this one
  void SystematicsMaker::Add(SystematicsMaker* systMaker, bool overwrite) {

    map<string, SystMaker*> systMakers;
    for (const auto& systMaker : systMakers) {
      if (fSystMakers.count(systMaker.first) and !overwrite)
        continue;
      fSystMakers[systMaker.first] = systMaker.second;
    }
  } // function SystematicsMaker::Add

  // -------------------------------------------------------------------------------------
  void SystematicsMaker::SaveTo(TDirectory* dir, bool separate) const {

    TDirectory* tmp = gDirectory;

    dir->cd();
    TObjString("SystematicsMaker").Write("type");

    TObjString(fName.c_str()).Write("name");

    int predType = -1;
    switch (fPredType) {
    case SystPredType::kExtrap:
      predType = 0;
      break;
    case SystPredType::kFD:
      predType = 1;
      break;
    case SystPredType::kND:
      predType = 2;
      break;
    default:
      break;
    }
    TVectorD predtype(1);
    predtype[0] = predType;
    predtype.Write("predtype");

    TDirectory* systMakersDir = dir->mkdir("SystMakers");
    for (const auto& systMaker : fSystMakers) {
      systMaker.second->SaveTo(systMakersDir->mkdir(systMaker.first.c_str()));
    }

    tmp->cd();

  } // function SystematicsMaker::SaveTo

  //----------------------------------------------------------------------
  unique_ptr<SystematicsMaker> SystematicsMaker::LoadFrom(
    TDirectory* dir) {

    DontAddDirectory guard;

    TObjString* tag = (TObjString*)dir->Get("type");
    assert(tag);
    assert(tag->GetString() == "SystematicsMaker");

    TObjString* tsName = (TObjString*)dir->Get("name");
    assert(tsName);
    string name(tsName->GetString().Data());

    const TVectorD predtype = *(TVectorD*)dir->Get("predtype");
    SystPredType predType = SystPredType::kUnknown;
    switch (static_cast<int>(std::round(predtype[0]))) {
    case 0:
      predType = SystPredType::kExtrap;
      break;
    case 1:
      predType = SystPredType::kFD;
      break;
    case 2:
      predType = SystPredType::kND;
      break;
    }

    unique_ptr<SystematicsMaker> systsMaker = make_unique<SystematicsMaker>(name, predType);

    map<string, SystMaker*> systMakers;
    TIter next(dir->GetDirectory("SystMakers")->GetListOfKeys());
    TKey* key;
    while ((key = (TKey*)next())) {
      SystMaker* systMaker
        = SystMaker::LoadFrom((TDirectory*)key->ReadObj()).release();
      systMakers[systMaker->GetName()] = systMaker;
    }

    systsMaker->fSystMakers = systMakers;

    return move(systsMaker);

  } // function SystematicsMaker::LoadFrom

} // namespace ana