PlotCovariance.C

Go to the documentation of this file.

/// \title  PlotCovariance.C
/// \author Connor Johnson
/// \date   Jan 16, 2020
/// \brief  A script which takes covariance matrices, and constructs our final pretty plots for NumuCC

#include <iostream>
#include <string>
#include <algorithm>

// ROOT includes
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TRatioPlot.h"
#include "THStack.h"
#include "TCanvas.h"
#include "TString.h"
#include "TSpline.h"
#include "TLegend.h"
#include "TLegendEntry.h"

#include "TMatrixD.h"
#include "TArray.h"
#include "TStyle.h"

#include "TMinuit.h"

// Needed ONLY for legend placement
#include "CAFAna/Analysis/Plots.h"

string cv_name = "cv";

// Binning logic

map<int,int> MuTperCosBins = 
{
  {1,  6 },
  {2,  7 },
  {3,  8 },
  {4,  9 },
  {5,  9 },
  {6,  9 },
  {7,  13},
  {8,  14},
  {9,  17},
  {10, 20},
  {11, 20},
  {12, 20},
  {13, 20}
};

const std::vector<double> angBinEdges {
  -1.0, 0.50, 0.56, 0.62, 0.68,
  0.74, 0.80, 0.85, 0.88, 0.91,
  0.94, 0.96, 0.98, 0.99, 1.0
};

const std::vector<double> angBinEdges_reduced {
  0.50, 0.56, 0.62, 0.68,
  0.74, 0.80, 0.85, 0.88, 0.91,
  0.94, 0.96, 0.98, 0.99, 1.0
};

const std::vector<double> mukeBinEdges {
  -10.0,
  0.5, 0.6, 0.7, 0.8, 0.9,
  1.0, 1.1, 1.2, 1.3, 1.4,
  1.5, 1.6, 1.7, 1.8, 1.9,
  2.0, 2.1, 2.2, 2.3, 2.4,
  2.5, 120
};

const std::vector<double> mukeBinEdges_reduced {
  0.5, 0.6, 0.7, 0.8, 0.9,
  1.0, 1.1, 1.2, 1.3, 1.4,
  1.5, 1.6, 1.7, 1.8, 1.9,
  2.0, 2.1, 2.2, 2.3, 2.4,
  2.5
};

const std::vector<double> enuBinEdges{
  // -10,
  //0.0, 0.25,
  0.50, 0.75,
  1.0, 1.25, 1.50, 1.75,
  2.0, 2.25, 2.50, 2.75,
  3.0, 3.25, 3.50, 3.75,
  4.0, 4.25, 4.50, 4.75,
  5.0//, 120
};

const std::vector<double> q2BinEdges {
  0.00,
  0.05, 0.10, 0.15, 0.20, 0.25,
  0.30, 0.40, 0.50, 0.60, 0.75,
  0.90, 1.10, 1.40, 1.80, 2.8//,
  //120.
};

// Axis titles
string mukeTitle  = "T_{#mu} (GeV)";
string muAngTitle = "cos(#theta_{#mu})";
string enuTitle   = "E_{#nu} (GeV)";
string q2Title    = "q^{2} (GeV^{2})";

// XSec titles
string xsec2DTitle  = "#frac{d^{2} #sigma}{d cos(#theta_{#mu}) d T_{#mu}} (cm^{2} / GeV / nucleon #times 10^{39})";
string xsecENuTitle = "#sigma / E_{#nu} (cm^{2} / GeV / nucleon #times 10^{39})";
string xsecQ2Title  = "#frac{d #sigma}{d q^{2}} (cm^{2} / GeV^{2} / nucleon #times 10^{39})";

TH2* Construct2D(TH1* histIn)
{
  TH2* histOut = new TH2D((histIn->GetName() + string("_2D")).c_str(), "", angBinEdges_reduced.size() - 1, angBinEdges_reduced.data(), mukeBinEdges_reduced.size() - 1, mukeBinEdges_reduced.data());

  int counter = 1;
  int NbinX = histOut->GetXaxis()->GetNbins();
  int NbinY = histOut->GetYaxis()->GetNbins();

  // Logic reverse-engineered from covariance making script
  for(int ibinx = 1; ibinx <= NbinX; ++ibinx){ // Throw out bin [-1, 0.5] in costheta
    for(int ibiny = 1; ibiny <= NbinY; ++ibiny){ // Throw out first _and_ last bin of MuKE (-10, 0.5), (2.5, 120)
      // Skip once we've reached how many MuKE bins to include in each costheta bin
      if(MuTperCosBins.find(ibinx) == MuTperCosBins.end() || ibiny > MuTperCosBins[ibinx])
        continue;
      histOut->SetBinContent(ibinx, ibiny, histIn->GetBinContent(counter));
      histOut->SetBinError(ibinx, ibiny, histIn->GetBinError(counter));
      counter++;
    }
  }
  histOut->GetXaxis()->SetRangeUser(0.5, 1.0);
  histOut->GetYaxis()->SetRangeUser(0.5, 2.5);
  return histOut;
}

TH1* ConstructENu(TH1* histIn)
{
  assert(histIn->GetNbinsX() == (int) enuBinEdges.size() - 1);
  TH1 * histOut = new TH1D((histIn->GetName() + string("_enuaxis")).c_str(), "", enuBinEdges.size() - 1, enuBinEdges.data());
  for (int ibinx = 1; ibinx <= histOut->GetNbinsX(); ++ibinx){
    histOut->SetBinContent(ibinx, histIn->GetBinContent(ibinx));
    histOut->SetBinError(ibinx, histIn->GetBinError(ibinx));
  }
  return histOut;
}

TH1* ConstructQ2(TH1* histIn)
{
  assert(histIn->GetNbinsX() == (int) q2BinEdges.size() - 1);
  TH1 * histOut = new TH1D((histIn->GetName() + string("_q2axis")).c_str(), "", q2BinEdges.size() - 1, q2BinEdges.data());
  for (int ibinx = 1; ibinx <= histOut->GetNbinsX(); ++ibinx){
    histOut->SetBinContent(ibinx, histIn->GetBinContent(ibinx));
    histOut->SetBinError(ibinx, histIn->GetBinError(ibinx));
  }
  return histOut;
}

// Function to create the map from covariance bin # to global bin number in 2D
TH1I * MakeCovTo2DBin(TH2 * hist, int nbins)
{
  TH1I * arr = new TH1I("covTo2DBin", "", nbins, 1, nbins + 1);
  int bincounter = 1;
  int NbinX = angBinEdges_reduced.size();
  int NbinY = mukeBinEdges_reduced.size();
  for(int ibinx = 1; ibinx <= NbinX; ++ibinx){ // Should this not be 2?
    for(int ibiny = 1; ibiny <= NbinY; ++ibiny){
      if(MuTperCosBins.find(ibinx) == MuTperCosBins.end() || ibiny > MuTperCosBins[ibinx])
        continue;
      arr->SetBinContent(bincounter, hist->GetBin(ibinx, ibiny));
      bincounter++;
    }
  }
  --bincounter; // Fix an off-by-one error for the very last valid bin incrementing.
  cout << "Counter=" << bincounter << " NBins=" << nbins << endl;
  assert(bincounter == nbins);
  return arr;
}

// Function to map x and y bin number to covariance bin number.
TH2I * Make2DToCovBin()
{
  TH2I * binmap = new TH2I("hist2DBinToCov", "", 
    angBinEdges_reduced.size() - 1, angBinEdges_reduced.data(),
    mukeBinEdges_reduced.size() - 1, mukeBinEdges_reduced.data());
  int bincounter = 1;
  int NbinX = angBinEdges_reduced.size();
  int NbinY = mukeBinEdges_reduced.size();
  for(int ibinx = 1; ibinx <= NbinX; ++ibinx){ // Should this not be 2?
    for(int ibiny = 1; ibiny <= NbinY; ++ibiny){
      if(MuTperCosBins.find(ibinx) == MuTperCosBins.end() || ibiny > MuTperCosBins[ibinx])
        continue;
      binmap->SetBinContent(ibinx, ibiny, bincounter);
      bincounter++;
    }
  }
  --bincounter; // Fix an off-by-one error for the very last valid bin incrementing.
  cout << "Counter=" << bincounter << endl;
  return binmap;
}

// Function to go from global bin number to bin in the covariance matrix
int GetCovBin(TH1I * covTo2DBin, int globalBin)
{
  assert(covTo2DBin && covTo2DBin->GetNbinsX() > 0);
  for(int i = 1; i <= covTo2DBin->GetNbinsX(); ++i)
    if (covTo2DBin->GetBinContent(i) == globalBin)
      return i;
  return -1;
}

// Set the errors of a histogram to the sqrt of diagonals of the corresponding covariance matrix
TH1 * SetErrorFromCovDiag(TH1* hist, TH2* mat)
{
  assert(hist->GetNbinsX() == mat->GetNbinsX());
  for (int ibin = 1; ibin <= mat->GetNbinsX(); ++ibin)
  {
    hist->SetBinError(ibin, sqrt(mat->GetBinContent(ibin, ibin))); // Fill with sqrt diagonals
  }
  return hist;
}

// Set the errors of a histogram to the sqrt of diagonals of the corresponding covariance matrix
TH1 * SetErrorFromCovDiag(TH1* hist, TMatrixD& mat)
{
  assert(hist->GetNbinsX() == mat.GetNrows());
  for (int irow = 0; irow < mat.GetNrows(); ++irow)
  {
    hist->SetBinError(irow+1, sqrt(mat[irow][irow])); // Fill with sqrt diagonals
  }
  return hist;
}

// Helpers for getting
enum EBinEdge{
  kBinLowEdge,
  kBinUpEdge
};

int GetBinByEdge(TAxis* axis, double binEdge, EBinEdge findBinBy = kBinLowEdge){
  double epsilon = 1E-5;
  int foundBin = -1;
  for(int ibin = 1; ibin <= axis->GetNbins(); ++ibin){
    double ibinEdge;
    if (findBinBy == kBinLowEdge)
      ibinEdge = axis->GetBinLowEdge(ibin);
    else if (findBinBy == kBinUpEdge)
      ibinEdge = axis->GetBinUpEdge(ibin);
    else
      throw runtime_error("EBinEdge has value " + to_string(findBinBy) + ", which is not defined.");
    if (abs(ibinEdge - binEdge) < epsilon){
      foundBin = ibin;
      break;
    }
  }
  if (foundBin < 1)
    throw runtime_error("Could not find bin matching " + to_string(binEdge));
  return foundBin;
}

// Overlay all generators over data
void OverlayGenerators(map<string, map<string, TH1* > > generators, string var, bool drawSpline = false, TLegend ** leg = NULL, double binLowEdge = -1.)
{
  int smoothFactor = 5;
  cout << "Starting bin checks for var: " << var << endl;
  for (pair<string, map<string, TH1* > > generator : generators){
    if (drawSpline)
    {
      TSpline3 * spline;
      if (var == "MuKin"){
        smoothFactor = 10;
        int generatorBin = GetBinByEdge(generator.second["MuKin"]->GetXaxis(), binLowEdge, kBinLowEdge);
        TH2 * hist = (TH2*) generator.second[var];
        TH1 * histSlice = hist->ProjectionY((hist->GetName() + string("_") + to_string(generatorBin)).c_str(), generatorBin, generatorBin);
        // cout << "Using bin low edge: " << hist->GetXaxis()->GetBinLowEdge(generatorBin) << " for generator." <<  endl;
        histSlice->Rebin(smoothFactor);
        histSlice->Scale(1.0/double(smoothFactor));
        // cout << "Generator: " << generator.first << " using generatorBin: " << generatorBin << " (" <<
          // hist->GetXaxis()->GetBinLowEdge(generatorBin) << " - " << hist->GetXaxis()->GetBinUpEdge(generatorBin) << ")" << endl;
        spline = new TSpline3(histSlice);
      }
      else{
        generator.second[var]->Rebin(smoothFactor);
        generator.second[var]->Scale(1.0/smoothFactor);
        spline = new TSpline3(generator.second[var]);
      }

      spline->SetLineColor(generator.second[var]->GetLineColor());
      spline->SetLineWidth(3);
      spline->Draw("same");
    }
    else
    {
      TH1 * hist = generator.second[var];
      if (var == "MuKin"){
        int generatorBin = GetBinByEdge(generator.second["MuKin"]->GetXaxis(), binLowEdge, kBinLowEdge);
        TH2 * hist2D = (TH2*) generator.second[var];
        hist = hist2D->ProjectionY((hist2D->GetName() + string("_") + to_string(generatorBin)).c_str(), generatorBin, generatorBin);
        // cout << "Generator: " << generator.first << " using generatorBin: " << generatorBin << " (" <<
          // hist2D->GetXaxis()->GetBinLowEdge(generatorBin) << " - " << hist2D->GetXaxis()->GetBinUpEdge(generatorBin) << ")" << endl;
      }
      hist->SetLineColor(generator.second[var]->GetLineColor());
      hist->SetLineWidth(3);
      hist->Draw("hist same");
    }

    if (leg)
      (*leg)->AddEntry(generator.second[var], generator.first.c_str());
  }
  if (leg){
    TLegend * oldLeg = (TLegend*) (*leg)->Clone();
    (*leg) = ana::AutoPlaceLegend(0.34, 0.2);
    (*leg)->SetFillStyle(0);
    (*leg)->SetTextSize(0.03);
    for (TObject * obj : *(oldLeg->GetListOfPrimitives())){
      TLegendEntry * entry = (TLegendEntry*) obj;
      (*leg)->AddEntry(entry->GetObject(), entry->GetLabel());
    }
  }
}

pair<TH1*, TH1*> GeneratorRatio(TH1 * gen, TH1 * data)
{
  TH1 * genClone = (TH1*) gen->Clone();
  TH1 * genRat   = (TH1*) gen->Clone();
  assert(gen->GetNbinsX() == data->GetNbinsX());
  genClone->SetTitle(data->GetTitle());
  genRat->Divide(data);
  genRat->GetYaxis()->SetTitle("Gen / Data");
  // genRat->SetLineWidth(2.);
  return {genClone, genRat};
}

void PlotGenRatios(map<string, map<string, TH1* > > generators, TH1* dataFull, TH1* dataStat, string var, TLegend ** leg = NULL, double binLowEdge = -1.)
{
  TRatioPlot * baseRatio = NULL;
  vector<pair<TH1*, TH1*> > ratios;
  for (pair<string, map<string, TH1* > > generator : generators){
    TH1 * hist = generator.second[var];
    if (var == "MuKin"){
      TH2 * hist2D = (TH2*) generator.second[var];
      int generatorBin = GetBinByEdge(hist2D->GetXaxis(), binLowEdge, kBinLowEdge);
      hist = hist2D->ProjectionY((hist->GetName() + string("_") + to_string(generatorBin)).c_str(), generatorBin, generatorBin);
      hist->SetLineColor(hist2D->GetLineColor());
    }
    pair<TH1*, TH1*> genRatio = GeneratorRatio(hist, dataFull);
    ratios.push_back(genRatio);
    if (!baseRatio){
      baseRatio = new TRatioPlot(hist, dataFull);
      // (*leg)->AddEntry(baseRatio->GetUpperRefObject(), generator.first.c_str());
    }
    if(leg)
      (*leg)->AddEntry(hist, generator.first.c_str());
  }

  double maxhist = (*std::max_element(ratios.begin(), ratios.end(), [](auto i, auto j){
      return (i.first->GetMaximum() < j.first->GetMaximum());
    })).first->GetMaximum();
  double maxrat  = (*std::max_element(ratios.begin(), ratios.end(), [](auto i, auto j){return i.second->GetMaximum() < j.second->GetMaximum();})).second->GetMaximum();
  double minrat  = (*std::min_element(ratios.begin(), ratios.end(), [](auto i, auto j){return i.second->GetMinimum() < j.second->GetMinimum();})).second->GetMinimum();

  baseRatio->Draw();
  baseRatio->SetSplitFraction(0.4);
  baseRatio->SetSeparationMargin(0.0);
  baseRatio->GetLowerRefYaxis()->SetTitle("Gen / Data");
  ((TH1*) baseRatio->GetUpperRefObject())->SetLineWidth(0.0);
  ((TH1*) baseRatio->GetUpperRefObject())->SetTitle("Generator Ratios");
  ((TH1*) baseRatio->GetUpperRefObject())->SetTitleSize(0.035, "Y");
  baseRatio->GetLowerRefGraph()->SetMarkerSize(0.0);
  baseRatio->GetLowerRefGraph()->GetYaxis()->SetTitleSize(0.035);
  baseRatio->SetGraphDrawOpt("X");
  baseRatio->GetUpperRefYaxis()->SetRangeUser(0.0, maxhist * 1.2);
  baseRatio->GetLowerRefYaxis()->SetRangeUser(minrat * 0.9, 1.1 * maxrat);
  baseRatio->GetUpperRefXaxis()->SetTitle(dataFull->GetXaxis()->GetTitle());
  baseRatio->GetLowYaxis()->SetTitle(dataFull->GetXaxis()->GetTitle());
  baseRatio->GetUpperRefYaxis()->SetTitle(dataFull->GetYaxis()->GetTitle());
  baseRatio->GetUpperRefYaxis()->CenterTitle(false);
  baseRatio->GetLowerRefYaxis()->CenterTitle(true);
  baseRatio->GetLowerRefYaxis()->SetTitleOffset(1.3);
  baseRatio->SetGridlines({});
  baseRatio->Draw("noconfint");
  for(pair<TH1*, TH1*> ratio : ratios)
  {
    TH1 * gen = ratio.first;
    TH1 * rat = ratio.second;
    baseRatio->GetUpperPad()->cd();
    gen->Draw("hist same");
    baseRatio->GetLowerPad()->cd();
    rat->Draw("hist same");
  }
  baseRatio->GetUpperPad()->cd();
  dataFull->Draw("e1 same");
  dataStat->Draw("e1 same");

  if (leg){
    TLegend * oldLeg = (TLegend*) (*leg)->Clone();
    (*leg) = ana::AutoPlaceLegend(0.34, var == "MuKin" ? 0.3 : 0.2);
    (*leg)->SetFillStyle(0);
    (*leg)->SetTextSize(0.03);
    for (TObject * obj : *(oldLeg->GetListOfPrimitives())){
      TLegendEntry * entry = (TLegendEntry*) obj;
      (*leg)->AddEntry(entry->GetObject(), entry->GetLabel());
    }
    (*leg)->Draw();
  }
}

// Take a TH2 and print all angle slices
void PrintAllSlices(TH2* histSystPlusStat, TH2* histStat, string histLabel, string filename, bool drawSpline, map<string, map<string, TH1* > > generators = {}, bool printRatios = false)
{
  TCanvas * canv = new TCanvas("sliceCanv", "sliceCanv", 1600, 1200);
  canv->SetLeftMargin(0.14);
  canv->SetRightMargin(0.05);
  canv->SetBottomMargin(0.12);
  TLegend * leg;
  for (int ibinx = 2; ibinx <= histSystPlusStat->GetNbinsX(); ++ibinx){
    TH1 * proj0 = histSystPlusStat->ProjectionY((histSystPlusStat->GetName() + string("_") + to_string(ibinx)).c_str(), ibinx, ibinx);
    proj0->SetTitle(TString::Format("%3.2f < cos(#theta_{#mu}) < %3.2f", histSystPlusStat->GetXaxis()->GetBinLowEdge(ibinx), histSystPlusStat->GetXaxis()->GetBinUpEdge(ibinx)));
    proj0->GetYaxis()->SetTitle(xsec2DTitle.c_str());
    proj0->GetXaxis()->SetRangeUser(0.5, 2.5);
    proj0->GetYaxis()->SetRangeUser(0.0, proj0->GetMaximum() * 1.25);
    proj0->SetMarkerStyle(histSystPlusStat->GetMarkerStyle());
    proj0->Draw("E1");
    leg = new TLegend();
    leg->SetFillStyle(0);
    leg->AddEntry(proj0, histLabel.c_str());
    proj0->Draw("E1 SAME");
    OverlayGenerators(generators, "MuKin", drawSpline, &leg, histSystPlusStat->GetXaxis()->GetBinLowEdge(ibinx));
    TH1 * projStat0 = histStat->ProjectionY((histStat->GetName() + string("_") + to_string(ibinx)).c_str(), ibinx, ibinx);
    projStat0->SetLineColor(kGray);
    projStat0->Draw("E1 SAME");
    leg->Draw();
    canv->SaveAs((filename + "_bin" + to_string(ibinx) + ".pdf").c_str());
    canv->Clear();
    if (printRatios)
    {
      leg->Clear();
      leg->AddEntry(proj0, histLabel.c_str());
      PlotGenRatios(generators, proj0, projStat0, "MuKin", &leg, histSystPlusStat->GetXaxis()->GetBinLowEdge(ibinx));
      leg->Draw();
      canv->SaveAs((filename + "_bin" + to_string(ibinx) + "_ratios.pdf").c_str());
      canv->Clear();
    }
    leg->Clear();
    delete leg;
  }
  canv->Clear();
  delete canv;
}

void SigmaDivideENu(TH1* hist)
{
  for (int ibinx = 1; ibinx <= hist->GetNbinsX(); ++ibinx){
    hist->SetBinContent(ibinx, hist->GetBinContent(ibinx) / hist->GetXaxis()->GetBinCenter(ibinx));
    hist->SetBinError(ibinx, hist->GetBinError(ibinx) / hist->GetXaxis()->GetBinCenter(ibinx));
  }
}

map<string, map<string, string> >  generator_hist_names
{
  {"genie3",
    {
      {"MuKin", "T_muCos_theta"},
      {"ENu",   "xsecVsEnu_2"},
      {"Q2",    "xsecVsq2_2"}
    }
  },
  {"genie3_rebin",
    {
      {"MuKin", "T_muCos_theta"},
      {"ENu",   "xsecVsEnu_2"},
      {"Q2",    "xsecVsq2VB_2"}
    }
  },
  {"default",
    {
      {"MuKin", "T_muCos_theta_Soup"},
      {"ENu",   "xsecVsEnu_2_Soup"},
      {"Q2",    "xsecVsq2_2_Soup"}
    }
  },
  {"rebin",
    {
      {"MuKin", "T_muCos_theta_Soup"},
      {"ENu",   "xsecVsEnu_2_Soup"},
      {"Q2",    "xsecVsq2VB_2_Soup"}
    }
  }
};

// Load the generator histograms from a file, scale if needed
map<string, TH1*> LoadGenerator(TFile * file, string label, Color_t col, double scaling = -1, string histLabel = "default")
{
  map<string, TH1*> generatorMap;

  assert(generator_hist_names.find(histLabel) != generator_hist_names.end());
  assert(generator_hist_names[histLabel].find("MuKin") != generator_hist_names[histLabel].end());
  assert(generator_hist_names[histLabel].find("ENu") != generator_hist_names[histLabel].end());
  assert(generator_hist_names[histLabel].find("Q2") != generator_hist_names[histLabel].end());

  generatorMap["MuKin"] = (TH2*) file->Get(generator_hist_names[histLabel]["MuKin"].c_str())->Clone();
  generatorMap["ENu"]   = (TH1*) file->Get(generator_hist_names[histLabel]["ENu"].c_str())->Clone();
  generatorMap["Q2"]    = (TH1*) file->Get(generator_hist_names[histLabel]["Q2"].c_str())->Clone();

  generatorMap["MuKin"]->SetLineColor(col);
  generatorMap["ENu"]  ->SetLineColor(col);
  generatorMap["Q2"]   ->SetLineColor(col);

  string appendLabel = histLabel == "rebin" || histLabel == "genie3_rebin" ? "rebin" : "";
  generatorMap["MuKin"]->SetName((label + "_MuKin" + appendLabel).c_str());
  generatorMap["ENu"]  ->SetName((label + "_ENu" +   appendLabel).c_str());
  generatorMap["Q2"]   ->SetName((label + "_Q2" +    appendLabel).c_str());

  if (scaling > 0)
  {
    generatorMap["MuKin"]->Scale(scaling);
    generatorMap["ENu"]  ->Scale(scaling);
    generatorMap["Q2"]   ->Scale(scaling);
  }

  generatorMap["MuKin"]->GetXaxis()->SetRangeUser(0.5, 1.0);
  generatorMap["MuKin"]->GetYaxis()->SetRangeUser(0.5, 2.5);
  generatorMap["ENu"]  ->GetXaxis()->SetRangeUser(0.5, 5.0);
  generatorMap["Q2"]   ->GetXaxis()->SetRangeUser(0.0, 2.8);

  // Set XSec of ENu to XSec / ENu of ENu
  SigmaDivideENu(generatorMap["ENu"]);

  return generatorMap;
}

map<string, TH1*> LoadGeniePrediction(TFile * file, string label, Color_t col, string weight)
{
  map<string, TH1*> geniePred;
  geniePred["MuKin"] = (TH2*) file->Get(("geniePrediction/" + weight + "/MuKin").c_str())->Clone();
  geniePred["ENu"]   = (TH1*) file->Get(("geniePrediction/" + weight + "/ENu").c_str())->Clone();
  geniePred["Q2"]    = (TH1*) file->Get(("geniePrediction/" + weight + "/Q2").c_str())->Clone();
  geniePred["MuKin"]->SetLineColor(col);
  geniePred["ENu"]  ->SetLineColor(col);
  geniePred["Q2"]   ->SetLineColor(col);
  geniePred["MuKin"]->SetName((label + "_MuKin").c_str());
  geniePred["ENu"]  ->SetName((label + "_ENu"  ).c_str());
  geniePred["Q2"]   ->SetName((label + "_Q2"   ).c_str());

  // Get rid of huge scaling factor
  geniePred["MuKin"]->Scale(1E39);
  geniePred["ENu"]  ->Scale(1E39);
  geniePred["Q2"]   ->Scale(1E39);

  geniePred["MuKin"]->GetXaxis()->SetRangeUser(0.5, 1.0);
  geniePred["MuKin"]->GetYaxis()->SetRangeUser(0.5, 2.5);
  geniePred["ENu"]  ->GetXaxis()->SetRangeUser(0.5, 5.0);
  geniePred["Q2"]   ->GetXaxis()->SetRangeUser(0.0, 2.8);

  // Set XSec of ENu to XSec / ENu of ENu
  SigmaDivideENu(geniePred["ENu"]);

  return geniePred;
}

void AreaNormalize(TH1* generator, TH1* data, const vector<double>& binning, const string label)
{
  double genBinLow  = generator->FindBin(binning[0]);
  double genBinUp   = generator->FindBin(binning[binning.size() - 1]) - 1;
  double dataBinLow = data->FindBin(binning[0]);
  double dataBinUp  = data->FindBin(binning[binning.size() - 1]) - 1;
  cout << "Scaling " << label << ": " << data->Integral(dataBinLow, dataBinUp, "width") / generator->Integral(genBinLow, genBinUp, "width") << endl;
  generator->Scale(data->Integral(dataBinLow, dataBinUp, "width") / generator->Integral(genBinLow, genBinUp, "width"));
}

// TH2 is slightly different. We need to normalize based on only the bins we present, not outside the phase space
void AreaNormalize2D(TH2* generator, TH2* data, const string label)
{
  double genIntegral  = 0.;
  double dataIntegral = 0.;
  for(int ibinx = 1; ibinx <= data->GetNbinsX(); ++ibinx)
    for(int ibiny = 1; ibiny <= data->GetNbinsY(); ++ibiny)
      if (data->GetBinContent(ibinx, ibiny) > 0.0){
        // cout << "\tFound (" << std::scientific << to_string(ibinx) << "," << to_string(ibiny) << ") in data.\n\tAdding: " 
          // << data->GetBinContent(ibinx, ibiny) << " * " << to_string(data->GetXaxis()->GetBinWidth(ibinx)) << " * " << to_string(data->GetYaxis()->GetBinWidth(ibiny)) << endl;
        dataIntegral += data->GetBinContent(ibinx, ibiny) * data->GetXaxis()->GetBinWidth(ibinx) * data->GetYaxis()->GetBinWidth(ibiny);
      }
  for(int ibinx = 1; ibinx <= generator->GetNbinsX(); ++ibinx)
    for(int ibiny = 1; ibiny <= generator->GetNbinsY(); ++ibiny){
        int dataBinX = data->GetXaxis()->FindBin(generator->GetXaxis()->GetBinCenter(ibinx));
        int dataBinY = data->GetYaxis()->FindBin(generator->GetYaxis()->GetBinCenter(ibiny));
        if (data->GetBinContent(dataBinX, dataBinY) > 0.0){
          // cout << "\tFound (" << std::scientific << to_string(ibinx) << "," << to_string(ibiny) << ") in generator.\n\tAdding: " 
            // << generator->GetBinContent(ibinx, ibiny) << " * " << to_string(generator->GetXaxis()->GetBinWidth(ibinx)) << " * " << to_string(generator->GetYaxis()->GetBinWidth(ibiny)) << endl;
          genIntegral += generator->GetBinContent(ibinx, ibiny) * generator->GetXaxis()->GetBinWidth(ibinx) * generator->GetYaxis()->GetBinWidth(ibiny);
        }
      }
  cout << "Scaling: " << label << ": " << dataIntegral / genIntegral << endl;
  generator->Scale(dataIntegral / genIntegral);
}

int Get2DBinArea(TH2* hist, int binxLow, int binxHigh, int binyLow, int binyHigh)
{
  double widthx = 0;
  double widthy = 0;
  for (int ibinx = binxLow; ibinx <= binxHigh; ++ibinx)
    widthx += hist->GetXaxis()->GetBinWidth(ibinx);
  for (int ibiny = binyLow; ibiny <= binyHigh; ++ibiny)
    widthy += hist->GetYaxis()->GetBinWidth(ibiny);
  return widthx * widthy;
}

map<string, map<string, TH1* > > RebinGenerators(map<string, map<string, TH1* > > generators)
{
  map<string, map<string, TH1* > > rebinnedGenerators;
  for (pair<string, map<string, TH1* > > generatorPair : generators){
    string genLabel = generatorPair.first;
    map<string, TH1* > & generator = generatorPair.second;

    // Initialize the empty histograms to the right binning
    TH2D * mukin_hist = new TH2D((generator["MuKin"]->GetName() + string("_databins")).c_str(), (genLabel + " Rebin").c_str(),
      angBinEdges_reduced.size() - 1, angBinEdges_reduced.data(), mukeBinEdges_reduced.size() - 1, mukeBinEdges_reduced.data());
    TH1D * enu_hist   = new TH1D((generator["ENu"]->GetName() + string("_databins")).c_str(), (genLabel + " Rebin").c_str(),
      enuBinEdges.size() - 1, enuBinEdges.data());
    TH1D * q2_hist    = new TH1D((generator["Q2"]->GetName() + string("_databins")).c_str(), (genLabel + " Rebin").c_str(),
      q2BinEdges.size() - 1, q2BinEdges.data());

    mukin_hist->SetLineColor(generator["MuKin"]->GetLineColor());
    enu_hist  ->SetLineColor(generator["ENu"]  ->GetLineColor());
    q2_hist   ->SetLineColor(generator["Q2"]->GetLineColor());

    // Loop through MuKin bin edges, fill with correct integral, scale by inverse bin width.
    for (unsigned int ibinx = 0; ibinx < angBinEdges_reduced.size() - 1; ++ibinx)
    {
      double binxEdgeLow = angBinEdges_reduced[ibinx];
      double binxEdgeUp  = angBinEdges_reduced[ibinx + 1];
      int binxLow        = GetBinByEdge(generator["MuKin"]->GetXaxis(), binxEdgeLow, kBinLowEdge);
      int binxHigh       = GetBinByEdge(generator["MuKin"]->GetXaxis(), binxEdgeUp, kBinUpEdge);
      for (unsigned int ibiny = 0; ibiny < mukeBinEdges_reduced.size() - 1; ++ibiny)
      {
        double binyEdgeLow = mukeBinEdges_reduced[ibiny];
        double binyEdgeUp  = mukeBinEdges_reduced[ibiny + 1];
        int binyLow        = GetBinByEdge(generator["MuKin"]->GetYaxis(), binyEdgeLow, kBinLowEdge);
        int binyHigh       = GetBinByEdge(generator["MuKin"]->GetYaxis(), binyEdgeUp, kBinUpEdge);
        double binIntegral = ((TH2*) generator["MuKin"])->Integral(binxLow, binxHigh, binyLow, binyHigh);
        int numBins = (binxHigh - binxLow + 1) * (binyHigh - binyLow + 1);
        mukin_hist->SetBinContent(ibinx + 1, ibiny + 1, binIntegral / numBins);
      }
    }

    // Simpler logic now for the ENu
    for (unsigned int ibinx = 0; ibinx < enuBinEdges.size() - 1; ++ibinx)
    {
      double binEdgeLow  = enuBinEdges[ibinx];
      double binEdgeUp   = enuBinEdges[ibinx + 1];
      int binLow         = GetBinByEdge(generator["ENu"]->GetXaxis(), binEdgeLow, kBinLowEdge);
      int binHigh        = GetBinByEdge(generator["ENu"]->GetXaxis(), binEdgeUp, kBinUpEdge);
      double binIntegral = generator["ENu"]->Integral(binLow, binHigh);
      enu_hist->SetBinContent(ibinx + 1, binIntegral / (binHigh - binLow + 1));
    }

    // Fill the Q2
    for (unsigned int ibinx = 0; ibinx < q2BinEdges.size() - 1; ++ibinx)
    {
      double binEdgeLow  = q2BinEdges[ibinx];
      double binEdgeUp   = q2BinEdges[ibinx + 1];
      int binLow         = GetBinByEdge(generator["Q2"]->GetXaxis(), binEdgeLow, kBinLowEdge);
      int binHigh        = GetBinByEdge(generator["Q2"]->GetXaxis(), binEdgeUp, kBinUpEdge);
      double binIntegral = generator["Q2"]->Integral(binLow, binHigh);
      q2_hist->SetBinContent(ibinx + 1, binIntegral / (binHigh - binLow + 1));
    }

    rebinnedGenerators[generatorPair.first]["MuKin"] = mukin_hist;
    rebinnedGenerators[generatorPair.first]["ENu"]   = enu_hist;
    rebinnedGenerators[generatorPair.first]["Q2"]    = q2_hist;
  }
  return rebinnedGenerators;
}

// Chi2 and minimization

double CalculateENuChi2(const TMatrixD& invCovMat, TH1 * data, TH1 * generator, string label)
{
  double chi2 = 0;
  double naiveChi2 = 0;
  assert(invCovMat.GetNrows() == invCovMat.GetNcols());
  assert(invCovMat.GetNrows() == (int) enuBinEdges.size() - 1);
  assert(invCovMat.GetNrows() == data->GetNbinsX());
  assert(invCovMat.GetNrows() == generator->GetNbinsX());
  TH2D * chi2Contribs = new TH2D((label + "_enu_chi2_contribs").c_str(), (label + " E_{#nu} #chi^{2} Contributions").c_str(), invCovMat.GetNrows(), 0, invCovMat.GetNrows(), invCovMat.GetNcols(), 0, invCovMat.GetNcols());
  for (int ibin = 0; ibin < invCovMat.GetNrows(); ++ibin)
    for (int jbin = 0; jbin < invCovMat.GetNcols(); ++jbin)
    {
      double binChi2 = (generator->GetBinContent(ibin + 1) - data->GetBinContent(ibin + 1)) * invCovMat[ibin][jbin] * (generator->GetBinContent(jbin + 1) - data->GetBinContent(jbin + 1));
      chi2Contribs->SetBinContent(ibin + 1, jbin + 1, binChi2);
      chi2 += binChi2;
      if (ibin == jbin) // Only diagonals for naive calculation
        naiveChi2 += pow(generator->GetBinContent(ibin + 1) - data->GetBinContent(ibin + 1), 2) / pow(data->GetBinError(ibin + 1), 2);
    }
  chi2Contribs->Draw("COLZ");
  cout << "\tNaive Chi2: " << naiveChi2 << endl;
  return chi2;
}

double Calculate1DChi2_simple(const TMatrixD& invCovMat, const TH1 * data, const TH1 * generator)
{
  double chi2 = 0;
  for (int ibin = 0; ibin < invCovMat.GetNrows(); ++ibin)
    for (int jbin = 0; jbin < invCovMat.GetNcols(); ++jbin)
      chi2 += (generator->GetBinContent(ibin + 1) - data->GetBinContent(ibin + 1)) * invCovMat[ibin][jbin] * (generator->GetBinContent(jbin + 1) - data->GetBinContent(jbin + 1));
  return chi2;
}

TMatrixD * chi2_min_enu_invCovMat;
TH1 *    chi2_min_enu_data;
TH1 *    chi2_min_enu_gen;
void fcn_enu(int &npar, double *gin, double &f, double *par, int iflag)
{
  cout << "Running fcn_enu() at scaling: " << to_string(par[0]) << endl;
  TH1 * scaled_gen = (TH1*) chi2_min_enu_gen->Clone("temp");
  scaled_gen->Scale(par[0]);
  f = Calculate1DChi2_simple(*chi2_min_enu_invCovMat, chi2_min_enu_data, scaled_gen);
  cout << "\tGives chi2: " << to_string(f) << endl;
  scaled_gen->Delete();
}

double CalculateQ2Chi2(const TMatrixD& invCovMat, TH1 * data, TH1 * generator, string label)
{
  double chi2 = 0;
  double naiveChi2 = 0;
  assert(invCovMat.GetNrows() == invCovMat.GetNcols());
  assert(invCovMat.GetNrows() == (int) q2BinEdges.size() - 1);
  assert(invCovMat.GetNrows() == data->GetNbinsX());
  assert(invCovMat.GetNrows() == generator->GetNbinsX());
  TH2D * chi2Contribs = new TH2D((label + "_q2_chi2_contribs").c_str(), (label + " q^{2} #chi^{2} Contributions").c_str(), invCovMat.GetNrows(), 0, invCovMat.GetNrows(), invCovMat.GetNcols(), 0, invCovMat.GetNcols());
  for (int ibin = 0; ibin < invCovMat.GetNrows(); ++ibin)
    for (int jbin = 0; jbin < invCovMat.GetNcols(); ++jbin)
    {
      double binChi2 = (generator->GetBinContent(ibin + 1) - data->GetBinContent(ibin + 1)) * invCovMat[ibin][jbin] * (generator->GetBinContent(jbin + 1) - data->GetBinContent(jbin + 1));
      chi2Contribs->SetBinContent(ibin + 1, jbin + 1, binChi2);
      chi2 += binChi2;
      if (ibin == jbin) // Only diagonals for naive calculation
        naiveChi2 += pow(generator->GetBinContent(ibin + 1) - data->GetBinContent(ibin + 1), 2) / pow(data->GetBinError(ibin + 1), 2);
    }
  chi2Contribs->Draw("COLZ");
  cout << "\tNaive Chi2: " << naiveChi2 << endl;
  return chi2;
}

TMatrixD * chi2_min_q2_invCovMat;
TH1 *    chi2_min_q2_data;
TH1 *    chi2_min_q2_gen;
void fcn_q2(int &npar, double *gin, double &f, double *par, int iflag)
{
  cout << "Running fcn_q2() at scaling: " << to_string(par[0]) << endl;
  TH1 * scaled_gen = (TH1*) chi2_min_q2_gen->Clone("temp");
  scaled_gen->Scale(par[0]);
  f = Calculate1DChi2_simple(*chi2_min_q2_invCovMat, chi2_min_q2_data, scaled_gen);
  cout << "\tGives chi2: " << to_string(f) << endl;
  scaled_gen->Delete();
}

double Calculate2DChi2(const TMatrixD& invCovMat, TH2 * data, TH2 * generator, TH1 * covToGlobalBinMap, string label)
{
  double chi2 = 0;
  double naiveChi2 = 0;
  assert(invCovMat.GetNrows() == invCovMat.GetNcols());
  assert(invCovMat.GetNrows() == covToGlobalBinMap->GetNbinsX());
  assert(data->GetNbinsX() == generator->GetNbinsX());
  assert(data->GetNbinsY() == generator->GetNbinsY());
  TH2D * chi2Contribs = new TH2D((label + "_mukin_chi2_contribs").c_str(), (label + "Muon Kinematic #chi^{2} Contributions").c_str(), invCovMat.GetNrows(), 0, invCovMat.GetNrows(), invCovMat.GetNcols(), 0, invCovMat.GetNcols());
  // Check the binnings are exact with edges
  for (int ibinx = 1; ibinx <= generator->GetNbinsX(); ++ibinx){
    for (int ibiny = 1; ibiny <= generator->GetNbinsY(); ++ibiny){
      double xval = generator->GetXaxis()->GetBinLowEdge(ibinx);
      double yval = generator->GetYaxis()->GetBinLowEdge(ibiny);
      if (find_if(angBinEdges.begin(), angBinEdges.end(), [xval](double b) {return abs(xval - b) < 1E-10;}) == angBinEdges.end() ||
          find_if(mukeBinEdges.begin(), mukeBinEdges.end(), [yval](double b) {return abs(yval - b) < 1E-10;}) == mukeBinEdges.end()){
        cerr << "\tDid not find bin: (" << generator->GetXaxis()->GetBinLowEdge(ibinx) << ", " << generator->GetYaxis()->GetBinLowEdge(ibiny) << ") in list of bins." << endl;
        cerr << "\tBinning (2D) does not match! Skipping generator." << endl;
        return -1;
      }
    }
  }
  for (int ibin = 0; ibin < invCovMat.GetNrows(); ++ibin)
    for (int jbin = 0; jbin < invCovMat.GetNcols(); ++jbin)
    {
      int iGlobalBin = covToGlobalBinMap->GetBinContent(ibin + 1);
      int jGlobalBin = covToGlobalBinMap->GetBinContent(jbin + 1);
      int ibinX = -1;
      int ibinY = -1;
      int jbinX = -1;
      int jbinY = -1;
      int temp  = -1;
      data->GetBinXYZ(iGlobalBin, ibinX, ibinY, temp); // Fill the binx and biny for global bin I
      data->GetBinXYZ(jGlobalBin, jbinX, jbinY, temp); // Fill the binx and biny for global bin J
      double binChi2 = (generator->GetBinContent(ibinX, ibinY) - data->GetBinContent(ibinX, ibinY)) * invCovMat[ibin][jbin] * (generator->GetBinContent(jbinX, jbinY) - data->GetBinContent(jbinX, jbinY));
      chi2Contribs->SetBinContent(ibin + 1, jbin + 1, binChi2);
      if (ibin == jbin) // Only diagonals for naive calculation
        naiveChi2 += pow(generator->GetBinContent(ibinX, ibinY) - data->GetBinContent(ibinX, ibinY), 2) / pow(data->GetBinError(ibinX, ibinY), 2);
      // cout << "\tUsing matrix bins (" << ibin << "," << jbin << ")" << endl;
      // cout << "\tCorresponding to MuKin:((" << ibinXCenter << ", " << ibinYCenter << "), (" << jbinXCenter << ", " << jbinYCenter << "))" << endl;
      // cout << "\tUsing generator bins ((" << ibinXGen << "," << ibinYGen << "), (" << jbinXGen << ", " << jbinYGen << "))" << endl;
      // cout << "\tUsing data bins ((" << ibinXData << "," << ibinYData << "), (" << jbinXData << ", " << jbinYData << "))" << endl;
      // cout << "\tAdding (" << generator->GetBinContent(ibinXGen, ibinYGen) << " - " << data->GetBinContent(ibinXData, ibinYData) << ") * " <<
        // invCovMat[ibin][jbin] << " * (" << generator->GetBinContent(jbinXGen, jbinYGen) << " - " << data->GetBinContent(jbinXData, jbinYData) << ")" << endl;
      // cout << "\t(=" << binChi2 << ")" << endl;
      chi2 += binChi2;
    }
  chi2Contribs->Draw("COLZ");
  cout << "\tNaive Chi2: " << naiveChi2 << endl;
  return chi2;
}

void PlotDataGenerators(TH1 * dataSystPlusStat, TH1 * dataStat, string dataTitle, map<string, map<string, TH1* > > &generatorMap, string varName, bool drawSpline = false)
{
  dataSystPlusStat->Draw("E1");
  TLegend * leg = new TLegend(0.6, 0.6, 0.9, 0.85);
  leg->SetFillStyle(0);
  leg->AddEntry(dataSystPlusStat, dataTitle.c_str());
  OverlayGenerators(generatorMap, varName, drawSpline, &leg);
  dataSystPlusStat->Draw("E1 SAME");
  dataStat->SetLineColor(kGray);
  dataStat->Draw("E1 SAME");
  leg->Draw();
}

TMatrixD StatMatFromHist(TH1* hist)
{
  int nrowcols = hist->GetNbinsX();
  TMatrixD statErrs(nrowcols, nrowcols);
  for(int irow = 0; irow < nrowcols; irow++)
    statErrs[irow][irow] = hist->GetBinError(irow+1);
  return statErrs;
}

void PlotCovMatrices(TMatrixD & cov, TMatrixD & inv, string outputFolder, string label, string title)
{
  TCanvas * tempCanv = new TCanvas("temp", "temp", 1600, 1200);
  tempCanv->SetRightMargin(0.15);

  TH2D * covHist = new TH2D(cov);
  covHist->SetName((label + "CovHist").c_str());
  covHist->SetTitle((title + " Covariance Matrix;;;#times 10^{-78}").c_str());
  covHist->Draw("COLZ");
  tempCanv->SaveAs((outputFolder + "/" + label + "CovMat.pdf").c_str());

  TH2D * invHist = new TH2D(inv);
  invHist->SetName((label + "InvCovHist").c_str());
  invHist->SetTitle((title + " Inv. Cov Matrix;;;#times 10^{78}").c_str());
  invHist->Draw("COLZ");
  tempCanv->SaveAs((outputFolder + "/" + label + "InvCovMat.pdf").c_str());

  TMatrixD identity = cov * inv;
  TH2D * identityHist = new TH2D(identity);
  identityHist->SetName((label + "Identity").c_str());
  identityHist->SetTitle((title + " Cov #times InvCov").c_str());
  identityHist->Draw("COLZ");
  tempCanv->SaveAs((outputFolder + "/" + label + "Identity.pdf").c_str());

  covHist->Write();
  invHist->Write();

  tempCanv->Clear();
  delete tempCanv;
}

void PlotCovariance(string outputFolder = ".")
{
  gStyle->SetTitleYSize(0.05);
  // gStyle->SetTitleZSize(0.15);

  TCanvas * c1 = new TCanvas("c1", "c1", 1600, 1200);
  string shapeOnlyFilename = "CovCurrent/Calculate_covariances_shape_all.root";
  string totalErrFilename = "CovCurrent/Calculate_covariances_total_all.root";

  TFile * outFile = TFile::Open((outputFolder + "/numuCCXSecResults.root").c_str(), "RECREATE");
  TDirectory * xsecDir = outFile->mkdir("results");
  TDirectory * generatorDir = outFile->mkdir("generators");

  // Get a 2D bin map of x-y to 
  TH2I * bin2DToCov = Make2DToCovBin();

  /////////////////////////////////
  // Load Statistical Errors
  /////////////////////////////////
  TFile * Stat2DFile = TFile::Open("/nova/ana/users/connorj/NumuCC/NERSC_Run_Dec2019/xsecUncertResults.root");
  TFile * Stat1DFile = TFile::Open("/nova/ana/users/connorj/NumuCC/Feb2020_MuonE_Angle_SystFix/xsecUncertResults_1DOnly.root");
  assert(Stat2DFile && Stat2DFile->IsOpen());
  assert(Stat1DFile && Stat1DFile->IsOpen());
  TH2 * statMuKin_originalbins    = (TH2*) Stat2DFile->Get("cv/xsec_cv_2DProj");
  TH1 * statENu_originalbins      = (TH1*) Stat1DFile->Get("cv/xsec_cv_ENu");
  TH1 * statQ2_originalbins       = (TH1*) Stat1DFile->Get("cv/xsec_cv_Q2");
  statMuKin_originalbins->Scale(1E39);
  statENu_originalbins->Scale(1E39);
  statQ2_originalbins->Scale(1E39);

  TH2 * statMuKin = new TH2D("statMuKin", ("Statistical Errors;" + muAngTitle + ";" + mukeTitle + ";" + xsec2DTitle).c_str(),
    angBinEdges_reduced.size() - 1, angBinEdges_reduced.data(),
    mukeBinEdges_reduced.size() - 1, mukeBinEdges_reduced.data());
  TH1 * statENu = statENu_originalbins->Rebin(enuBinEdges.size() - 1, "statENu", enuBinEdges.data());
  TH1 * statQ2  = statQ2_originalbins ->Rebin(q2BinEdges.size() - 1,  "statQ2",  q2BinEdges.data());
  // Loop through MuKin bin edges, fill with correct integral, scale by inverse bin width.
  for (unsigned int ibinx = 0; ibinx < angBinEdges_reduced.size() - 1; ++ibinx)
  {
    double binxEdgeLow = angBinEdges_reduced[ibinx];
    double binxEdgeUp  = angBinEdges_reduced[ibinx + 1];
    int binxLow        = GetBinByEdge(statMuKin_originalbins->GetXaxis(), binxEdgeLow, kBinLowEdge);
    int binxHigh       = GetBinByEdge(statMuKin_originalbins->GetXaxis(), binxEdgeUp, kBinUpEdge);
    for (unsigned int ibiny = 0; ibiny < mukeBinEdges_reduced.size() - 1; ++ibiny)
    {
      double binyEdgeLow = mukeBinEdges_reduced[ibiny];
      double binyEdgeUp  = mukeBinEdges_reduced[ibiny + 1];
      int binyLow        = GetBinByEdge(statMuKin_originalbins->GetYaxis(), binyEdgeLow, kBinLowEdge);
      int binyHigh       = GetBinByEdge(statMuKin_originalbins->GetYaxis(), binyEdgeUp, kBinUpEdge);
      statMuKin->SetBinContent(ibinx + 1, ibiny + 1, statMuKin_originalbins->GetBinContent(binxLow, binyLow));
      statMuKin->SetBinError(ibinx + 1, ibiny + 1, statMuKin_originalbins->GetBinError(binxLow, binyLow));
    }
  }

  // Add the statistical errors to the diagonals of covariance matrix
  TDirectory * statErrDir = xsecDir->mkdir("statErrs");
  statErrDir->cd();
  statMuKin->Write();
  statENu->Write();
  statQ2->Write();

  TMatrixD statMatMuKin(bin2DToCov->GetMaximum(), bin2DToCov->GetMaximum());
  for(int ibin = 0; ibin < statMatMuKin.GetNrows(); ++ibin){
    unsigned int binx = 0;
    unsigned int biny = 0;
    for(int ibinx = 1; ibinx <= bin2DToCov->GetNbinsX(); ++ibinx)
      for(int ibiny = 1; ibiny <= bin2DToCov->GetNbinsY(); ++ibiny)
        if (ibin == bin2DToCov->GetBinContent(ibinx, ibiny))
        {
          binx = ibinx;
          biny = ibiny;
        }
    assert(binx > 0);
    assert(biny > 0);
    statMatMuKin[ibin][ibin] = statMuKin->GetBinError(binx, biny);
  }
  TMatrixD statMatENu(StatMatFromHist(statENu));
  TMatrixD statMatQ2 (StatMatFromHist(statQ2));

  /////////////////////////////////
  // Load Total Error
  /////////////////////////////////

  // Load the nominal hists
  TFile * totalErrFile = TFile::Open(totalErrFilename.c_str());
  TH1 * nomMuKinTotal1D       = (TH1*) totalErrFile->Get("hcv_mukin");
  TH1 * nomENuTotal_nobininfo = (TH1*) totalErrFile->Get("hcv_enu");
  TH1 * nomQ2Total_nobininfo  = (TH1*) totalErrFile->Get("hcv_q2");
  nomMuKinTotal1D      ->Scale(1E39);
  nomENuTotal_nobininfo->Scale(1E39);
  nomQ2Total_nobininfo ->Scale(1E39);

  // Load the covariance matrices
  TH2 * covTotalMuKin = (TH2*) totalErrFile->Get("cov_mukin");
  TH2 * covTotalENu   = (TH2*) totalErrFile->Get("cov_enu");
  TH2 * covTotalQ2    = (TH2*) totalErrFile->Get("cov_q2");
  covTotalMuKin->Scale(1E78);
  covTotalENu  ->Scale(1E78);
  covTotalQ2   ->Scale(1E78);

  // Make matricies and inverses
  TMatrixD covMatTotalMuKin(covTotalMuKin->GetNbinsX(), covTotalMuKin->GetNbinsY());
  TMatrixD covMatTotalENu  (covTotalENu  ->GetNbinsX(), covTotalENu  ->GetNbinsY());
  TMatrixD covMatTotalQ2   (covTotalQ2   ->GetNbinsX(), covTotalQ2   ->GetNbinsY());
  for (int ibinx = 1; ibinx <= covTotalMuKin->GetNbinsX(); ++ibinx)
    for (int ibiny = 1; ibiny <= covTotalMuKin->GetNbinsY(); ++ibiny)
      covMatTotalMuKin[ibinx - 1][ibiny - 1] = covTotalMuKin->GetBinContent(ibinx, ibiny);
  for (int ibinx = 1; ibinx <= covTotalENu->GetNbinsX(); ++ibinx)
    for (int ibiny = 1; ibiny <= covTotalENu->GetNbinsY(); ++ibiny)
      covMatTotalENu[ibinx - 1][ibiny - 1] = covTotalENu->GetBinContent(ibinx, ibiny);
  for (int ibinx = 1; ibinx <= covTotalQ2->GetNbinsX(); ++ibinx)
    for (int ibiny = 1; ibiny <= covTotalQ2->GetNbinsY(); ++ibiny)
      covMatTotalQ2[ibinx - 1][ibiny - 1] = covTotalQ2->GetBinContent(ibinx, ibiny);
  TMatrixD covMatTotalMuKinFull = covMatTotalMuKin + statMatMuKin;
  TMatrixD covMatTotalENuFull   = covMatTotalENu + statMatENu;
  TMatrixD covMatTotalQ2Full    = covMatTotalQ2 + statMatQ2;
  TMatrixD invCovMatTotalMuKin = covMatTotalMuKinFull;
  TMatrixD invCovMatTotalENu   = covMatTotalENuFull;
  TMatrixD invCovMatTotalQ2    = covMatTotalQ2Full;
  invCovMatTotalMuKin.Invert();
  invCovMatTotalENu  .Invert();
  invCovMatTotalQ2   .Invert();

  // Set nominal error bars to the sqrt of covariance diagonals
  SetErrorFromCovDiag(nomMuKinTotal1D      , covMatTotalMuKinFull);
  SetErrorFromCovDiag(nomENuTotal_nobininfo, covMatTotalENuFull);
  SetErrorFromCovDiag(nomQ2Total_nobininfo , covMatTotalQ2Full);

  // Load the right binnings from the nominal
  TH2 * nomMuKinTotal = Construct2D(nomMuKinTotal1D);
  TH1 * nomENuTotal   = ConstructENu(nomENuTotal_nobininfo);
  TH1 * nomQ2Total    = ConstructQ2(nomQ2Total_nobininfo);
  SigmaDivideENu(nomENuTotal); // ENu then gets scaled.

  // Load a second copy with the stats errors included
  TH1 * nomMuKinTotalStats1D       = (TH1*) nomMuKinTotal1D      ->Clone("hmukin_cv_2D_stats");
  TH1 * nomENuTotalStats_nobininfo = (TH1*) nomENuTotal_nobininfo->Clone("henu_cv_enuaxis_stats");
  TH1 * nomQ2TotalStats_nobininfo  = (TH1*) nomQ2Total_nobininfo ->Clone("hq2_cv_q2axis_stats");
  SetErrorFromCovDiag(nomMuKinTotalStats1D      , statMatMuKin);
  SetErrorFromCovDiag(nomENuTotalStats_nobininfo, statMatENu);
  SetErrorFromCovDiag(nomQ2TotalStats_nobininfo , statMatQ2);
  TH2 * nomMuKinStats = Construct2D(nomMuKinTotalStats1D);
  TH1 * nomENuStats   = ConstructENu(nomENuTotalStats_nobininfo);
  TH1 * nomQ2Stats    = ConstructQ2(nomQ2TotalStats_nobininfo);
  SigmaDivideENu(nomENuStats); // ENu then gets scaled.

  // Set titles, axis, markers, etc.
  nomMuKinTotal->SetTitle((";" + muAngTitle + ";" + mukeTitle + ";" + xsec2DTitle).c_str());
  nomENuTotal->  SetTitle((";" + enuTitle + ";" + xsecENuTitle).c_str());
  nomQ2Total->   SetTitle((";" + q2Title + ";" + xsecQ2Title).c_str());
  nomMuKinTotal->SetMarkerStyle(20);
  nomENuTotal->  SetMarkerStyle(20);
  nomQ2Total->   SetMarkerStyle(20);

  // Save to the output
  TDirectory * totalUncertResults = xsecDir->mkdir("totalUncertResults");
  totalUncertResults->cd();
  nomMuKinTotal->Write("MuKin");
  nomENuTotal->Write("ENu");
  nomQ2Total->Write("Q2");

  /////////////////////////////////
  // Load Shape-only Error
  /////////////////////////////////

  // Load the nominal hists
  TFile * shapeOnlyFile = TFile::Open(shapeOnlyFilename.c_str());
  TH1 * nomMuKinShape1D       = (TH1*) shapeOnlyFile->Get("hcv_mukin");
  TH1 * nomENuShape_nobininfo = (TH1*) shapeOnlyFile->Get("hcv_enu");
  TH1 * nomQ2Shape_nobininfo  = (TH1*) shapeOnlyFile->Get("hcv_q2");
  nomMuKinShape1D      ->Scale(1E39);
  nomENuShape_nobininfo->Scale(1E39);
  nomQ2Shape_nobininfo ->Scale(1E39);

  // Load the covariance matrices
  TH2 * covShapeMuKin = (TH2*) shapeOnlyFile->Get("cov_mukin");
  TH2 * covShapeENu   = (TH2*) shapeOnlyFile->Get("cov_enu");
  TH2 * covShapeQ2    = (TH2*) shapeOnlyFile->Get("cov_q2");
  covShapeMuKin->Scale(1E78);
  covShapeENu  ->Scale(1E78);
  covShapeQ2   ->Scale(1E78);

  // Make matricies and inverseShape
  TMatrixD covMatShapeMuKin(covShapeMuKin->GetNbinsX(), covShapeMuKin->GetNbinsY());
  TMatrixD covMatShapeENu  (covShapeENu  ->GetNbinsX(), covShapeENu  ->GetNbinsY());
  TMatrixD covMatShapeQ2   (covShapeQ2   ->GetNbinsX(), covShapeQ2   ->GetNbinsY());
  for (int ibinx = 1; ibinx <= covShapeMuKin->GetNbinsX(); ++ibinx)
    for (int ibiny = 1; ibiny <= covShapeMuKin->GetNbinsY(); ++ibiny)
      covMatShapeMuKin[ibinx - 1][ibiny - 1] = covShapeMuKin->GetBinContent(ibinx, ibiny);
  for (int ibinx = 1; ibinx <= covShapeENu->GetNbinsX(); ++ibinx)
    for (int ibiny = 1; ibiny <= covShapeENu->GetNbinsY(); ++ibiny)
      covMatShapeENu[ibinx - 1][ibiny - 1] = covShapeENu->GetBinContent(ibinx, ibiny);
  for (int ibinx = 1; ibinx <= covShapeQ2->GetNbinsX(); ++ibinx)
    for (int ibiny = 1; ibiny <= covShapeQ2->GetNbinsY(); ++ibiny)
      covMatShapeQ2[ibinx - 1][ibiny - 1] = covShapeQ2->GetBinContent(ibinx, ibiny);
  TMatrixD covMatShapeMuKinFull  = covMatShapeMuKin + statMatMuKin;
  TMatrixD covMatShapeENuFull    = covMatShapeENu + statMatENu;
  TMatrixD covMatShapeQ2Full     = covMatShapeQ2 + statMatQ2;
  TMatrixD invCovMatShapeMuKin = covMatShapeMuKinFull;
  TMatrixD invCovMatShapeENu   = covMatShapeENuFull;
  TMatrixD invCovMatShapeQ2    = covMatShapeQ2Full;
  invCovMatShapeMuKin.Invert();
  invCovMatShapeENu  .Invert();
  invCovMatShapeQ2   .Invert();

  // Plot and save covariances
  TDirectory * covarianceDir = outFile->mkdir("covariance");
  covarianceDir->cd();
  PlotCovMatrices(covMatTotalENuFull,   invCovMatTotalENu,   outputFolder, "enuTotal",   "E_{#nu} Total");
  PlotCovMatrices(covMatTotalQ2Full,    invCovMatTotalQ2,    outputFolder, "q2Total",    "q^{2}  Total");
  PlotCovMatrices(covMatTotalMuKinFull, invCovMatTotalMuKin, outputFolder, "mukinTotal", "Muon Kinematic Total");
  PlotCovMatrices(covMatShapeENuFull,   invCovMatShapeENu,   outputFolder, "enuShape",   "E_{#nu} Shape-only");
  PlotCovMatrices(covMatShapeQ2Full,    invCovMatShapeQ2,    outputFolder, "q2Shape",    "q^{2}  Shape-only");
  PlotCovMatrices(covMatShapeMuKinFull, invCovMatShapeMuKin, outputFolder, "mukinShape", "Muon Kinematic Shape-only");

  // Set nominal error bars to the sqrt of covariance diagonals
  SetErrorFromCovDiag(nomMuKinShape1D      , covMatShapeMuKinFull);
  SetErrorFromCovDiag(nomENuShape_nobininfo, covMatShapeENuFull);
  SetErrorFromCovDiag(nomQ2Shape_nobininfo , covMatShapeQ2Full);

  // Load the right binnings from the nominal
  TH2 * nomMuKinShape = Construct2D(nomMuKinShape1D);
  TH1 * nomENuShape   = ConstructENu(nomENuShape_nobininfo);
  TH1 * nomQ2Shape    = ConstructQ2(nomQ2Shape_nobininfo);
  SigmaDivideENu(nomENuShape); // ENu then gets scaled.

  // // Load a second copy with the stats errors included
  // TH1 * nomMuKinShapeStats1D       = (TH1*) nomMuKinShape1D      ->Clone("hmukin_cv_2D_shape_stats");
  // TH1 * nomENuShapeStats_nobininfo = (TH1*) nomENuShape_nobininfo->Clone("henu_cv_enuaxis_shape_stats");
  // TH1 * nomQ2ShapeStats_nobininfo  = (TH1*) nomQ2Shape_nobininfo ->Clone("hq2_cv_q2axis_shape_stats");
  // SetErrorFromCovDiag(nomMuKinShapeStats1D      , covMatShapeMuKinStats);
  // SetErrorFromCovDiag(nomENuShapeStats_nobininfo, covMatShapeENuStats);
  // SetErrorFromCovDiag(nomQ2ShapeStats_nobininfo , covMatShapeQ2Stats);
  // TH2 * nomMuKinShapeStats = Construct2D(nomMuKinShapeStats1D);
  // TH1 * nomENuShapeStats   = ConstructENu(nomENuShapeStats_nobininfo);
  // TH1 * nomQ2ShapeStats    = ConstructQ2(nomQ2ShapeStats_nobininfo);
  // SigmaDivideENu(nomENuShapeStats); // ENu then gets scaled.

  // Set titles, axis, markers, etc.
  nomMuKinShape->SetTitle((";" + muAngTitle + ";" + mukeTitle + ";" + xsec2DTitle).c_str());
  nomENuShape->  SetTitle((";" + enuTitle + ";" + xsecENuTitle).c_str());
  nomQ2Shape->   SetTitle((";" + q2Title + ";" + xsecQ2Title).c_str());
  nomMuKinShape->SetMarkerStyle(20);
  nomENuShape->  SetMarkerStyle(20);
  nomQ2Shape->   SetMarkerStyle(20);

  // Save to the output
  TDirectory * shapeUncertResults = xsecDir->mkdir("shapeUncertResults");
  shapeUncertResults->cd();
  nomMuKinShape->Write("MuKin");
  nomENuShape->Write("ENu");
  nomQ2Shape->Write("Q2");

  /////////////////////////////////
  // Make map of covariance bin to global bin
  /////////////////////////////////
  TH1I * covTo2DBin = MakeCovTo2DBin(nomMuKinTotal, covMatTotalMuKin.GetNrows());
  outFile->cd();
  covTo2DBin->Write("CovTo2DGlobalBin");
  bin2DToCov->Write("Bin2DToCov");

  /////////////////////////////////
  // Load the Generators
  /////////////////////////////////

  cout << "Started loading generators" << endl;

  // Load Generator Predictions
  TFile * gibuuFile = TFile::Open("GiBUU_numuCC_predictions.root");
  TFile * neutFile  = TFile::Open("NEUT54_numuCC_predictions.root");
  TFile * nuwroFile = TFile::Open("NuWro_numuCC_predictions.root");
  TFile * genieFile = TFile::Open("GeniePrediction.root");
  TFile * genie3File = TFile::Open("Genie3Predictions.root");
  assert(gibuuFile && gibuuFile->IsOpen());
  assert(neutFile  && neutFile->IsOpen());
  assert(nuwroFile && nuwroFile->IsOpen());
  assert(genieFile && genieFile->IsOpen());
  assert(genie3File && genie3File->IsOpen());
  map<string, map<string, TH1* > > generators;
  generators["GiBUU 2019"] = LoadGenerator(gibuuFile, "gibuu", kBlue, 1E39);
  generators["NEUT 5.4.0"]  = LoadGenerator(neutFile, "neut",  kGreen, 10);
  generators["NuWro 2019"] = LoadGenerator(nuwroFile, "nuwro", kViolet, 1E39);
  generators["GENIE 3.00.06"] = LoadGeniePrediction(genie3File, "genie3", kOrange, "ppfx");
  generators["GENIE 2.12.2 - Untuned"] = LoadGeniePrediction(genieFile,   "novatune", kCyan, "ppfx");
  generators["GENIE 2.12.2 - NOvA Tune"] = LoadGeniePrediction(genieFile, "untuned",  kRed, "both");

  // Copy of each generator that has been normalized to area.
  map<string, map<string, TH1* > > generatorsAreaNorm;
  for(pair<string, map<string, TH1* > > generator : generators)
    for(pair<string, TH1*> var : generator.second)
      generatorsAreaNorm[generator.first + " (area norm.)"][var.first] = (TH1*) var.second->Clone((var.second->GetName() + string("_norm")).c_str());
  for (pair<string, map<string, TH1* > > generator : generatorsAreaNorm)
  {
    AreaNormalize2D((TH2*) generator.second["MuKin"], nomMuKinShape, generator.first + " 2D");
    AreaNormalize(generator.second["ENu"],   nomENuShape, enuBinEdges, generator.first + " ENu");
    AreaNormalize(generator.second["Q2"],    nomQ2Shape, q2BinEdges, generator.first + " Q2");
  }

  // Rebinned generators
  map<string, map<string, TH1* > > generatorsVarBin;
  generatorsVarBin["GiBUU 2019"] = LoadGenerator(gibuuFile, "gibuu", kBlue, 1E39, "rebin");
  generatorsVarBin["NEUT 5.4.0"]  = LoadGenerator(neutFile, "neut",  kGreen, 10, "rebin");
  generatorsVarBin["NuWro 2019"] = LoadGenerator(nuwroFile, "nuwro", kViolet, 1E39, "rebin");
  generatorsVarBin["GENIE 3.00.06"] = LoadGeniePrediction(genie3File, "genie3_rebin", kOrange, "ppfx");
  generatorsVarBin["GENIE 2.12.2 - Untuned"] = LoadGeniePrediction(genieFile,   "novatune_rebin", kCyan, "ppfx"); // Genie are identical
  generatorsVarBin["GENIE 2.12.2 - NOvA Tune"] = LoadGeniePrediction(genieFile, "untuned_rebin",  kRed, "both");
  // Need to hardcode fix a few things
  generatorsVarBin["GiBUU 2019"]["ENu"]->Rebin(5);
  generatorsVarBin["GiBUU 2019"]["ENu"]->Scale(1./5);
  generatorsVarBin["GiBUU 2019"]["ENu"]->GetXaxis()->SetRangeUser(0.0, 5.0);
  generatorsVarBin["NuWro 2019"]["ENu"]->Rebin(5);
  generatorsVarBin["NuWro 2019"]["ENu"]->Scale(1./5);
  // generatorsVarBin["NuWro 2019"]["Q2"]->Scale(1.0/22);
  // generatorsVarBin["NEUT 5.4.0"]["Q2"]->Scale(1.0/30);
  generatorsVarBin["NEUT 5.4.0"]["ENu"]->Rebin(10);
  generatorsVarBin["NEUT 5.4.0"]["ENu"]->Scale(1./10);

  // Copy of each generator that has been normalized to area.
  map<string, map<string, TH1* > > generatorsAreaNormVarBin;
  for(pair<string, map<string, TH1* > > generator : generatorsVarBin)
    for(pair<string, TH1*> var : generator.second)
      generatorsAreaNormVarBin[generator.first + " (area norm.)"][var.first] = (TH1*) var.second->Clone((var.second->GetName() + string("_norm")).c_str());
  for (pair<string, map<string, TH1* > > generator : generatorsAreaNormVarBin)
  {
    AreaNormalize2D((TH2*) generator.second["MuKin"], nomMuKinShape, generator.first + " 2D Rebin");
    AreaNormalize(generator.second["ENu"],   nomENuShape, enuBinEdges, generator.first + " ENu Rebin");
    AreaNormalize(generator.second["Q2"],    nomQ2Shape, q2BinEdges, generator.first + " Q2 Rebin");
  }
  c1->SetRightMargin(0.15);

  // Rebinned generators
  generatorDir->cd();
  map<string, map<string, TH1* > > generatorsRebin = RebinGenerators(generatorsVarBin);
  map<string, map<string, TH1* > > generatorsAreaNormRebin = RebinGenerators(generatorsAreaNormVarBin);
  generatorDir->mkdir("standard")->cd();
  for (auto generator : generators)
    for (auto var : generator.second)
      var.second->Write((generator.first + "_" + var.first).c_str());
  generatorDir->mkdir("normalized")->cd();
  for (auto generator : generatorsAreaNorm)
    for (auto var : generator.second)
      var.second->Write((generator.first + "_" + var.first + "_norm").c_str());
  generatorDir->mkdir("databinned")->cd();
  for (auto generator : generatorsRebin)
    for (auto var : generator.second){
      var.second->Write((generator.first + "_" + var.first + "_databins").c_str());
    }
  generatorDir->mkdir("databinned_normalized")->cd();
  for (auto generator : generatorsAreaNormRebin)
    for (auto var : generator.second){
      var.second->Write((generator.first + "_" + var.first + "_databins_norm").c_str());
    }
  // Delete the ENu from Chi2
  // generatorsRebin["NuWro 2019"]->Delete("ENu");

  /////////////////////////////////
  // Calculate Chi2s
  /////////////////////////////////
  map<string, map<string, double > > generatorChi2s;
  map<string, map<string, TH1* > >   chi2_min_scaled_generators;
  map<string, map<string, double > > generatorScaledChi2s;
  map<string, string> gen_short_forms
  {
    {"GiBUU 2019", "GiBUU"},
    {"NEUT 5.4.0", "NEUT"},
    {"NuWro 2019", "NuWro"},
    {"GENIE 2.12.2 - Untuned", "Genie_Untuned"},
    {"GENIE 2.12.2 - NOvA Tune", "Genie_Tuned"},
    {"GENIE 3.00.06", "Genie_3"}
  };
  // for (auto generatorPair : generatorsRebin){
  //   double enuChi2 = CalculateENuChi2(invCovMatTotalENu, nomENuTotal, generatorPair.second["ENu"], gen_short_forms[generatorPair.first]);
  //   generatorChi2s[generatorPair.first]["ENu"] = enuChi2;
  //   c1->SaveAs((outputFolder + "/enu_chi2_contribs_" + gen_short_forms[generatorPair.first] + ".pdf").c_str());

  //   // TMinuit * enu_minimize = new TMinuit(1);
  //   // chi2_min_enu_invCovMat = new TMatrixD(invCovMatTotalENu);
  //   // chi2_min_enu_data      = nomENuTotal;
  //   // chi2_min_enu_gen       = generatorPair.second["ENu"];
  //   // enu_minimize->SetFCN(fcn_enu);
  //   // enu_minimize->DefineParameter(0, "scale", 1.0, 10.0, 0., 0.);
  //   // enu_minimize->Migrad();
  //   // double enu_scale, enu_scale_err;
  //   // enu_minimize->GetParameter(0, enu_scale, enu_scale_err);
  //   // cout << gen_short_forms[generatorPair.first] + " ENu scaling: " << enu_scale << " (error = " << enu_scale_err << ")" << endl;
  //   // TH1 * scaled_enu_gen = (TH1*)generatorPair.second["ENu"]->Clone((generatorPair.second["ENu"]->GetName() + string("_minchi2")).c_str());
  //   // scaled_enu_gen->Scale(enu_scale);
  //   // chi2_min_scaled_generators[generatorPair.first]["ENu"] = scaled_enu_gen;
  //   // generatorScaledChi2s[generatorPair.first]["ENu"] = CalculateENuChi2(invCovMatTotalENu, nomENuTotal, scaled_enu_gen, gen_short_forms[generatorPair.first] + "_chi2min");
  //   // c1->SaveAs((outputFolder + "/enu_chi2_contribs_" + gen_short_forms[generatorPair.first] + "_minchi2.pdf").c_str());

  //   double q2Chi2 = CalculateQ2Chi2(invCovMatTotalQ2, nomQ2Total, generatorPair.second["Q2"], gen_short_forms[generatorPair.first]);
  //   generatorChi2s[generatorPair.first]["Q2"] = q2Chi2;
  //   c1->SaveAs((outputFolder + "/q2_chi2_contribs_" + gen_short_forms[generatorPair.first] + ".pdf").c_str());

  //   // TMinuit * q2_minimize = new TMinuit(1);
  //   // chi2_min_q2_invCovMat = new TMatrixD(invCovMatTotalQ2);
  //   // chi2_min_q2_data      = nomQ2Total;
  //   // chi2_min_q2_gen       = generatorPair.second["Q2"];
  //   // q2_minimize->SetFCN(fcn_q2);
  //   // q2_minimize->DefineParameter(0, "scale", 1.0, 10.0, 0., 0.);
  //   // q2_minimize->Migrad();
  //   // double q2_scale, q2_scale_err;
  //   // q2_minimize->GetParameter(0, q2_scale, q2_scale_err);
  //   // cout << gen_short_forms[generatorPair.first] + " Q2 scaling: " << q2_scale << " (error = " << q2_scale_err << ")" << endl;
  //   // TH1 * scaled_q2_gen = (TH1*)generatorPair.second["Q2"]->Clone((generatorPair.second["Q2"]->GetName() + string("_minchi2")).c_str());
  //   // scaled_q2_gen->Scale(q2_scale);
  //   // chi2_min_scaled_generators[generatorPair.first]["Q2"] = scaled_q2_gen;
  //   // generatorScaledChi2s[generatorPair.first]["Q2"] = CalculateQ2Chi2(invCovMatTotalQ2, nomQ2Total, scaled_q2_gen, gen_short_forms[generatorPair.first] + "_chi2min");
  //   // c1->SaveAs((outputFolder + "/q2_chi2_contribs_" + gen_short_forms[generatorPair.first] + "_minchi2.pdf").c_str());

  //   cout << "2D Chi2: " << generatorPair.first << endl;
  //   double mukinChi2 = Calculate2DChi2(invCovMatTotalMuKin, nomMuKinTotal, (TH2*) generatorPair.second["MuKin"], covTo2DBin, gen_short_forms[generatorPair.first]);
  //   cout << "Found 2D chi2 = " << mukinChi2 << " for generator: " << generatorPair.first << endl;
  //   generatorChi2s[generatorPair.first]["MuKin"] = mukinChi2;
  //   c1->SaveAs((outputFolder + "/mukin_chi2_contribs_" + gen_short_forms[generatorPair.first] + ".pdf").c_str());
  // }
  // map<string, map<string, double > > generatorNormChi2s;
  // for (auto generatorPair : generatorsAreaNormRebin){
  //   string gen_short_form = gen_short_forms[generatorPair.first.substr(0, generatorPair.first.size() - 13)] + "_areanorm";

  //   cout << "ENu Chi2: " << generatorPair.first << endl;
  //   double enuChi2 = CalculateENuChi2(invCovMatShapeENu, nomENuShape, generatorPair.second["ENu"], gen_short_form);
  //   cout << "Found ENu chi2 = " << enuChi2 << " for generator: " << generatorPair.first << endl;
  //   generatorNormChi2s[generatorPair.first]["ENu"] = enuChi2;
  //   c1->SaveAs((outputFolder + "/enu_chi2_contribs_" + gen_short_form + ".pdf").c_str());

  //   cout << "Q2 Chi2: " << generatorPair.first << endl;
  //   double q2chi2 = CalculateQ2Chi2(invCovMatShapeQ2, nomQ2Shape, generatorPair.second["Q2"], gen_short_form);
  //   cout << "Found Q2 chi2 = " << q2chi2 << " for generator: " << generatorPair.first << endl;
  //   generatorNormChi2s[generatorPair.first]["Q2"] = q2chi2;
  //   c1->SaveAs((outputFolder + "/q2_chi2_contribs_" + gen_short_form + ".pdf").c_str());

  //   cout << "2D Chi2: " << generatorPair.first << endl;
  //   double mukinChi2 = Calculate2DChi2(invCovMatShapeMuKin, nomMuKinShape, (TH2*) generatorPair.second["MuKin"], covTo2DBin, gen_short_form);
  //   cout << "Found 2D chi2 = " << mukinChi2 << " for generator: " << generatorPair.first << endl;
  //   generatorNormChi2s[generatorPair.first]["MuKin"] = mukinChi2;
  //   c1->SaveAs((outputFolder + "/mukin_chi2_contribs_" + gen_short_form + ".pdf").c_str());
  // }

  /////////////////////////////////
  // Make plots
  /////////////////////////////////

  string dataTotalTitle = "Data - Total error";
  string dataShapeTitle = "Data - Shape only errs";  

  // Canvas
  TLegend * leg;

  // Draw the full 2D XSec to check
  c1->SetRightMargin(0.15);
  nomMuKinTotal->Draw("COLZ");
  c1->SaveAs((outputFolder + "/xsec_2dProj.pdf").c_str());

  // Tweak the margins again
  c1->SetLeftMargin(0.13);
  c1->SetRightMargin(0.05);
  c1->SetBottomMargin(0.10);

  // 2D Errors
  PrintAllSlices(nomMuKinTotal, nomMuKinStats, dataTotalTitle, outputFolder + "/muKinSlices",            true,  generators);
  PrintAllSlices(nomMuKinShape, nomMuKinStats, dataShapeTitle, outputFolder + "/muKinSlices_norm",       true,  generatorsAreaNorm);
  PrintAllSlices(nomMuKinTotal, nomMuKinStats, dataTotalTitle, outputFolder + "/muKinSlices_rebin",      false, generatorsRebin, true);
  PrintAllSlices(nomMuKinShape, nomMuKinStats, dataShapeTitle, outputFolder + "/muKinSlices_norm_rebin", false, generatorsAreaNormRebin, true);

  // ENu Plots
  nomENuTotal->GetYaxis()->SetRangeUser(0.0, nomENuTotal->GetMaximum() * 1.35);
  PlotDataGenerators(nomENuTotal, nomENuStats, dataTotalTitle, generators, "ENu", true);
  c1->SaveAs((outputFolder + "/results_enu.pdf").c_str());
  PlotDataGenerators(nomENuShape, nomENuStats, dataShapeTitle, generatorsAreaNorm, "ENu", true);
  c1->SaveAs((outputFolder + "/results_enu_norm.pdf").c_str());
  PlotDataGenerators(nomENuTotal, nomENuStats, dataTotalTitle, generatorsRebin, "ENu", false);
  c1->SaveAs((outputFolder + "/results_enu_rebin.pdf").c_str());
  PlotDataGenerators(nomENuShape, nomENuStats, dataShapeTitle, generatorsAreaNormRebin, "ENu", false);
  c1->SaveAs((outputFolder + "/results_enu_norm_rebin.pdf").c_str());

  // Q2 Plots
  nomQ2Total->GetYaxis()->SetRangeUser(0.0, nomQ2Total->GetMaximum() * 1.25);
  nomQ2Shape->GetYaxis()->SetRangeUser(0.0, nomQ2Shape->GetMaximum() * 1.25);
  PlotDataGenerators(nomQ2Total, nomQ2Stats, dataTotalTitle, generators, "Q2", true);
  c1->SaveAs((outputFolder + "/results_q2.pdf").c_str());
  PlotDataGenerators(nomQ2Shape, nomQ2Stats, dataShapeTitle, generatorsAreaNorm, "Q2", true);
  c1->SaveAs((outputFolder + "/results_q2_norm.pdf").c_str());
  PlotDataGenerators(nomQ2Total, nomQ2Stats, dataTotalTitle, generatorsRebin, "Q2", false);
  c1->SaveAs((outputFolder + "/results_q2_rebin.pdf").c_str());
  PlotDataGenerators(nomQ2Shape, nomQ2Stats, dataShapeTitle, generatorsAreaNormRebin, "Q2", false);
  c1->SaveAs((outputFolder + "/results_q2_norm_rebin.pdf").c_str());

  // Ratios
  c1->Clear();
  delete c1;
  c1 = new TCanvas("c1", "c1", 1600, 1320);
  leg = new TLegend();
  leg->AddEntry(nomENuTotal, dataTotalTitle.c_str());
  PlotGenRatios(generatorsRebin, nomENuTotal, nomENuStats, "ENu", &leg);
  c1->SaveAs((outputFolder + "/results_enu_ratios.pdf").c_str());

  c1->Clear();
  leg->Clear();
  delete leg;
  leg = new TLegend();
  leg->AddEntry(nomENuShape, dataShapeTitle.c_str());
  PlotGenRatios(generatorsAreaNormRebin, nomENuShape, nomENuStats, "ENu", &leg);
  c1->SaveAs((outputFolder + "/results_enu_shape_ratios.pdf").c_str());

  c1->Clear();
  leg->Clear();
  delete leg;
  leg = new TLegend();
  leg->AddEntry(nomQ2Total, dataTotalTitle.c_str());
  PlotGenRatios(generatorsRebin, nomQ2Total, nomQ2Stats, "Q2", &leg);
  c1->SaveAs((outputFolder + "/results_q2_ratios.pdf").c_str());

  c1->Clear();
  leg->Clear();
  delete leg;
  leg = new TLegend();
  leg->AddEntry(nomQ2Shape, dataShapeTitle.c_str());
  PlotGenRatios(generatorsAreaNormRebin, nomQ2Shape, nomQ2Stats, "Q2", &leg);
  c1->SaveAs((outputFolder + "/results_q2_shape_ratios.pdf").c_str());

  outFile->Close();
}