wrong_sign_uncertainty.C

Go to the documentation of this file.

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Vars/XsecTunes.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/XSec/Multiverse.h"
#include "CAFAna/XSec/SystematicDef.h"
#include "CAFAna/XSec/GenieMultiverseSyst.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Analysis/Exposures.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Systs/BeamSysts.h"


#include "TCanvas.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TLegendEntry.h"
#include "TPaletteAxis.h"
#include "THStack.h"

#include "NDAna/nuebarcc_inc/NuebarCCIncCuts.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncVars.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncBins.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncExtra.h"
// temporary measure to use prod4 trained muonid
#include "NDAna/muonid/NDXSecMuonPID.h"
namespace nuebarccinc {
  const ana::Cut kMuonIDProd4Cut = ana::kMuonIDProd4 < -0.2;
}

void load_libs_muonid()
{
  const int ret = gSystem->Load("libNDAnamuonid");
  if(ret != 0) exit(1);
}

using namespace ana;
TH2 * ToTH2(const TH1 * h, const HistAxis & xaxis, const HistAxis & yaxis);
TH1 * MostConservative(TH1 * up, TH1 * down, TH1* nominal);
void PlotMultiverse(Multiverse & mv, TH1 * nominal, std::string title, std::string name, TH1 * error = 0);
void Plot1D(TH1* h1, std::string h1_label, std::map<std::string, TH1*> systs, std::vector<std::string> to_plot, std::string name, double xmax=-1);
void Plot2D(TH1* total, std::map<std::string, TH1*> systs, std::vector<std::string> to_plot, std::string basename);
void abs(TH1* hist);
void MovePalette(TH1 * hist);

std::vector<std::string> all_systs  = {"ppfx", "genie", "calib_shape", "calib_shift", "ll", "cherenkov", "beam_focusing", "stat"};
std::vector<std::string> beam_systs = {"ppfx", "beam_focusing"};

std::map<std::string, nuebarccinc::Style> syst_styles = {{"genie"    , {"", kCyan    , 1}},
                                                         {"ppfx"     , {"", kMagenta , 1}},
                                                         {"ll"       , {"", kOrange  , 1}},
                                                         {"calib_shift"    , {"", kGreen   , 1}},
                                                         {"cherenkov", {"", kBlue+1  , 1}},
                                                         {"ll_up"    , {"", kOrange  , 9}},
                                                         {"ll_dw"    , {"", kOrange  , 2}},
                                                         {"calib_shift_up" , {"", kGreen   , 9}},
                                                         {"calib_shift_dw" , {"", kGreen   , 2}},
                                                         {"calib_shape" , {"", kRed-4   , 1}},
                                                         {"stat"     , {"", kBlack   , 2}},
                                                         
                                                         {"2kA"      , {"", kCyan, 1}},
                                                         {"02mmBeamSpotSize",         {"", kCyan + 1, 1}},
                                                         {"1mmBeamShiftX",            {"", kMagenta + 2, 1}},
                                                         {"1mmBeamShiftY",            {"", kMagenta+2,  2}},
                                                         {"3mmHorn1X",                {"", kViolet, 1}},
                                                         {"3mmHorn1Y",                {"", kViolet, 2}},
                                                         {"3mmHorn2X",                {"", kAzure, 1}},
                                                         {"3mmHorn2Y",                {"", kAzure, 2}},
                                                         {"7mmTargetZ",               {"", kBlue-6, 1}},
                                                         {"MagneticFieldinDecayPipe", {"", kGreen+2, 1}},
                                                         {"1mmHornWater",             {"", kYellow-1, 1}},
                                                         {"beam_focusing",            {"", kViolet+1}}
                                                         
};



std::map<std::string, double> xmax = {{"electron_kin", -1.},
                                      {"enu"         ,  6.},
                                      {"q2"          ,  2.}};

std::map<string, const ana::BeamSyst*> beamSystematics = {
  {"2kA",                      &ana::kBeamHornCurrent},
  {"02mmBeamSpotSize",         &ana::kBeamSpotSize},
  {"1mmBeamShiftX",            &ana::kBeamPosX},
  {"1mmBeamShiftY",            &ana::kBeamPosY},
  {"3mmHorn1X",                &ana::kBeamH1PosX},
  {"3mmHorn1Y",                &ana::kBeamH1PosY},
  {"3mmHorn2X",                &ana::kBeamH2PosX},
  {"3mmHorn2Y",                &ana::kBeamH2PosY},
  {"7mmTargetZ",               &ana::kBeamTarget},
  {"MagneticFieldinDecayPipe", &ana::kBeamMagField},
  {"1mmHornWater",             &ana::kBeamGeomWater}
};

void wrong_sign_uncertainty(bool make_plots = true,
                            std::string input_file_name = "/nova/ana/users/ddoyle/NuebarCCInc/WrongSign/wrong_sign_uncertainty.root",
                            std::string plot_dump = "/nova/ana/users/ddoyle/NuebarCCInc/WrongSign/plots")
{
  load_libs_muonid();
  std::map<std::string, const HistAxis*> axes;
  axes["electron_kin"] = new HistAxis(nuebarccinc::kTrueElectronEVsCosStandardAxis);
  axes["enu"         ] = new HistAxis(nuebarccinc::kTrueNeutrinoEnergyStandardAxis);
  axes["q2"          ] = new HistAxis(nuebarccinc::kTrueQ2StandardAxis            );

  std::map<std::string, const Cut*> cuts;
  cuts["nue"    ] = new Cut(nuebarccinc::kNueCC   );
  cuts["nuebar" ] = new Cut(nuebarccinc::kNuebarCC);
  cuts["nunubar"] = new Cut(nuebarccinc::kNuebarCC || nuebarccinc::kNueCC);

  std::map<std::string, SpectrumLoader*> loaders;
  loaders["nominal" ]     = 0;
  loaders["calib_shift_up"]     = 0;
  loaders["calib_shift_dw"]     = 0;
  loaders["ll_up"]        = 0;
  loaders["ll_dw"]        = 0;
  loaders["cherenkov"]    = 0;
  loaders["calib_shape"]  = 0;
  if(!make_plots) {
    loaders["nominal" ]    = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_NOMINAL);
    loaders["calib_shift_up"]    = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_UP);
    loaders["calib_shift_dw"]    = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_DOWN);
    loaders["calib_shape"] = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CALIB_SHAPE);
    loaders["ll_up"]       = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_LIGHT_UP);
    loaders["ll_dw"]       = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_LIGHT_DOWN);
    loaders["cherenkov"]   = new SpectrumLoader(nuebarccinc::PROD5_MC_RHC_CHERENKOV);
    const Var kCVWeight = kPPFXFluxCVWgt * kXSecCVWgt2020;

    const Cut kSelection = 
      nuebarccinc::kPreselectionLoose &&
      nuebarccinc::kMuonIDProd4Cut;

    std::map<std::string, std::map<std::string, std::map<std::string, Spectrum*> > > spectra;
    for(auto load_it = loaders.begin(); load_it != loaders.end(); load_it++) {
      load_it->second->SetSpillCut(kStandardSpillCuts);
      for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
        for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
          spectra[load_it->first][cut_it->first][ax_it->first] = new Spectrum(*load_it->second,
                                                                              *ax_it->second,
                                                                              kSelection && *cut_it->second,
                                                                              kNoShift,
                                                                              kCVWeight);
        }
      }
    }

    // beam systs
    for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
      for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
        for(auto beam_it = beamSystematics.begin(); beam_it != beamSystematics.end(); beam_it++) {
          auto name_up = beam_it->first + "_up";
          auto name_dw = beam_it->first + "_dw";

          spectra[name_up][cut_it->first][ax_it->first] = new Spectrum(*loaders["nominal"],
                                                                       *ax_it->second,
                                                                       kSelection && *cut_it->second,
                                                                       ana::SystShifts(beam_it->second, +1),
                                                                       kCVWeight);

          spectra[name_dw][cut_it->first][ax_it->first] = new Spectrum(*loaders["nominal"],
                                                                       *ax_it->second,
                                                                       kSelection && *cut_it->second,
                                                                       ana::SystShifts(beam_it->second, -1),
                                                                       kCVWeight);
        }
      }
    }

    const std::vector<Var> ppfx = GetkPPFXFluxUnivWgt();
    GenieMultiverseParameters genie_multiverse(100, getAllXsecSysts_2020());
    std::vector<SystShifts> genie_shifts = genie_multiverse.GetSystShifts();

    std::map<std::string, std::map<std::string, Multiverse *> > mv_ppfx;
    std::map<std::string, std::map<std::string, Multiverse *> > mv_genie;
    for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
      for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
        mv_ppfx[cut_it->first][ax_it->first] = new Multiverse(*loaders["nominal"],
                                                              *ax_it->second,
                                                              kSelection && *cut_it->second,
                                                              kNoShift,
                                                              ppfx,
                                                              kXSecCVWgt2020);
        mv_genie[cut_it->first][ax_it->first] = new Multiverse(*loaders["nominal"],
                                                               *ax_it->second,
                                                               kSelection && *cut_it->second,
                                                               genie_shifts,
                                                               kPPFXFluxCVWgt);
      }
    }

    for(auto load_it = loaders.begin(); load_it != loaders.end(); load_it++) 
      load_it->second->Go();


    TFile* output = new TFile("wrong_sign_uncertainty.root", "recreate");
    for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
      for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
        mv_ppfx .at(cut_it->first).at(ax_it->first)->SaveTo(output, ("ppfx_"  + cut_it->first + "_" + ax_it->first).c_str());
        mv_genie.at(cut_it->first).at(ax_it->first)->SaveTo(output, ("genie_" + cut_it->first + "_" + ax_it->first).c_str());
        for(auto load_it = loaders.begin(); load_it != loaders.end(); load_it++) {
          spectra.at(load_it->first).at(cut_it->first).at(ax_it->first)->SaveTo(output, TString::Format("%s_%s_%s",
                                                                                                        load_it->first.c_str(),
                                                                                                        cut_it ->first.c_str(),
                                                                                                        ax_it  ->first.c_str()).Data());
        }
        for(auto beam_it = beamSystematics.begin(); beam_it != beamSystematics.end(); beam_it++) {
          spectra.at(beam_it->first + "_up").at(cut_it->first).at(ax_it->first)->SaveTo(output, TString::Format("%s_%s_%s",
                                 cut_it  ->first.c_str(),
                                                                                                                ax_it   ->first.c_str()).Data());
          spectra.at(beam_it->first + "_dw").at(cut_it->first).at(ax_it->first)->SaveTo(output, TString::Format("%s_%s_%s",
                                                                                                                (beam_it->first + "_dw").c_str(),
                                                                                                                cut_it  ->first.c_str(),
                                                                                                                ax_it   ->first.c_str()).Data());
        }
      }
    }

    
  
    output->Close();
  }
  else {
    TFile * input = TFile::Open(input_file_name.c_str());   

    double MCPOT = ((TH1*) input->Get("nominal_nue_q2/pot"))->GetBinContent(1);
    std::map<std::string, double> pots;
    std::map<std::string, std::map<std::string, Multiverse *> > mv_ppfx;
    std::map<std::string, std::map<std::string, Multiverse *> > mv_genie;
    std::map<std::string, std::map<std::string, std::map<std::string, TH1*> > > hists;
    for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
      for(auto cut_it = cuts.begin(); cut_it != cuts.end(); cut_it++) {
        mv_ppfx [ax_it->first][cut_it->first] = Multiverse::LoadFrom(input, ("ppfx_"  + cut_it->first + "_" + ax_it->first).c_str())).release(;
        mv_genie[ax_it->first][cut_it->first] = Multiverse::LoadFrom(input, ("genie_" + cut_it->first + "_" + ax_it->first).c_str())).release(;

        mv_ppfx .at(ax_it->first).at(cut_it->first)->Scale(kAna2020RHCPOT / MCPOT);
        mv_genie.at(ax_it->first).at(cut_it->first)->Scale(kAna2020RHCPOT / MCPOT);

        for(auto load_it = loaders.begin(); load_it != loaders.end(); load_it++) {
          auto spec = Spectrum::LoadFrom(input, TString::Format("%s_%s_%s",
                                                                load_it->first.c_str(),
                                                                cut_it ->first.c_str(),
                                                                ax_it  ->first.c_str()).Data());
          pots[load_it->first] = spec->POT();
          auto pot = kAna2020RHCPOT;//spec->POT();
          hists[ax_it->first][cut_it->first][load_it->first] = spec->ToTH1(pot);
        }
        
        for(auto beam_it = beamSystematics.begin(); beam_it != beamSystematics.end(); beam_it++) {
          auto spec_up = Spectrum::LoadFrom(input, TString::Format("%s_%s_%s",
                                                                   (beam_it->first + "_up").c_str(),
                                                                   cut_it ->first         .c_str(),
                                                                   ax_it  ->first         .c_str()).Data());
          auto spec_dw = Spectrum::LoadFrom(input, TString::Format("%s_%s_%s",
                                                                   (beam_it->first + "_dw").c_str(),
                                                                   cut_it ->first         .c_str(),
                                                                   ax_it  ->first         .c_str()).Data());

          hists[ax_it->first][cut_it->first][beam_it->first + "_up"] = spec_up->ToTH1(kAna2020RHCPOT);
          hists[ax_it->first][cut_it->first][beam_it->first + "_dw"] = spec_dw->ToTH1(kAna2020RHCPOT);
        }
      }
    }
    
    // make fractions
    std::map<std::string, std::map<std::string, TH1*> > nuebar_fractions;
    std::map<std::string, Multiverse*> mv_ppfx_nuebar_fractions;
    std::map<std::string, Multiverse*> mv_genie_nuebar_fractions;
    std::map<std::string, Multiverse*> mv_ppfx_nue_nuebar_ratios;
    std::map<std::string, Multiverse*> mv_genie_nue_nuebar_ratios;
    std::map<std::string, std::map<std::string, TH1*> > nue_nuebar_ratios;
    for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
      mv_ppfx_nuebar_fractions [ax_it->first] = new Multiverse(*mv_ppfx .at(ax_it->first).at("nuebar"));
      mv_genie_nuebar_fractions[ax_it->first] = new Multiverse(*mv_genie.at(ax_it->first).at("nuebar"));
     
      *mv_ppfx_nuebar_fractions.at(ax_it->first)  /= *mv_ppfx .at(ax_it->first).at("nunubar");
      *mv_genie_nuebar_fractions.at(ax_it->first) /= *mv_genie.at(ax_it->first).at("nunubar");

      mv_ppfx_nue_nuebar_ratios[ax_it->first]  = new Multiverse(*mv_ppfx.at(ax_it->first).at("nue"));
      *mv_ppfx_nue_nuebar_ratios.at(ax_it->first) /= *mv_ppfx.at(ax_it->first).at("nuebar");

      mv_genie_nue_nuebar_ratios[ax_it->first]  = new Multiverse(*mv_genie.at(ax_it->first).at("nue"));
      *mv_genie_nue_nuebar_ratios.at(ax_it->first) /= *mv_genie.at(ax_it->first).at("nuebar");

      for(auto loader_it = loaders.begin(); loader_it != loaders.end(); loader_it++) {  
        nuebar_fractions[ax_it->first][loader_it->first] = (TH1*) hists.at(ax_it->first).at("nuebar").at(loader_it->first)->Clone(UniqueName().c_str());
        nuebar_fractions.at(ax_it->first).at(loader_it->first)->Divide(hists.at(ax_it->first).at("nunubar").at(loader_it->first));
        
        nue_nuebar_ratios[ax_it->first][loader_it->first] = (TH1*) hists.at(ax_it->first).at("nue").at(loader_it->first)->Clone(UniqueName().c_str());
        nue_nuebar_ratios.at(ax_it->first).at(loader_it->first)->Divide(hists.at(ax_it->first).at("nuebar").at(loader_it->first));
      }

      // two sided file systs
      for(std::string syst : {"calib_shift", "ll"}) {
        TH1 * tmp_up = (TH1*) hists.at(ax_it->first).at("nuebar").at(syst + "_up")->Clone(UniqueName().c_str());
        TH1 * tmp_dw = (TH1*) hists.at(ax_it->first).at("nuebar").at(syst + "_dw")->Clone(UniqueName().c_str());
        tmp_up->Divide(hists.at(ax_it->first).at("nunubar").at(syst + "_up"));
        tmp_dw->Divide(hists.at(ax_it->first).at("nunubar").at(syst + "_dw"));

        nuebar_fractions[ax_it->first][syst] = MostConservative(tmp_up, tmp_dw, nuebar_fractions.at(ax_it->first).at("nominal"));

        delete tmp_up;
        delete tmp_dw;

        tmp_up = (TH1*) hists.at(ax_it->first).at("nue").at(syst + "_up")->Clone(UniqueName().c_str());
        tmp_dw = (TH1*) hists.at(ax_it->first).at("nue").at(syst + "_dw")->Clone(UniqueName().c_str());
        tmp_up->Divide(hists.at(ax_it->first).at("nuebar").at(syst + "_up"));
        tmp_dw->Divide(hists.at(ax_it->first).at("nuebar").at(syst + "_dw"));

        nue_nuebar_ratios[ax_it->first][syst] = MostConservative(tmp_up, tmp_dw, nue_nuebar_ratios.at(ax_it->first).at("nominal"));
      }
      
      // treat all beam systs as two sided even though connor said there are some that are one sided. 
      for(auto beam_it = beamSystematics.begin(); beam_it != beamSystematics.end(); beam_it++) {
        TH1 * tmp_up = (TH1*) hists.at(ax_it->first).at("nuebar").at(beam_it->first + "_up")->Clone(UniqueName().c_str());
        TH1 * tmp_dw = (TH1*) hists.at(ax_it->first).at("nuebar").at(beam_it->first + "_dw")->Clone(UniqueName().c_str());
        tmp_up->Divide(hists.at(ax_it->first).at("nunubar").at(beam_it->first + "_up"));
        tmp_dw->Divide(hists.at(ax_it->first).at("nunubar").at(beam_it->first + "_dw"));

        nuebar_fractions[ax_it->first][beam_it->first] = MostConservative(tmp_up, tmp_dw, nuebar_fractions.at(ax_it->first).at("nominal"));

        delete tmp_up;
        delete tmp_dw;

        tmp_up = (TH1*) hists.at(ax_it->first).at("nue")  .at(beam_it->first + "_up")->Clone(UniqueName().c_str());
        tmp_dw = (TH1*) hists.at(ax_it->first).at("nue")  .at(beam_it->first + "_dw")->Clone(UniqueName().c_str());
        tmp_up->Divide(hists.at(ax_it->first).at("nuebar").at(beam_it->first + "_up"));
        tmp_dw->Divide(hists.at(ax_it->first).at("nuebar").at(beam_it->first + "_dw"));

        nue_nuebar_ratios[ax_it->first][beam_it->first] = MostConservative(tmp_up, tmp_dw, nue_nuebar_ratios.at(ax_it->first).at("nominal"));

      }

      // cherenkov
      nuebar_fractions.at(ax_it->first)["cherenkov"] = (TH1*) hists.at(ax_it->first).at("nuebar").at("cherenkov")->Clone(UniqueName().c_str());
      nuebar_fractions.at(ax_it->first).at("cherenkov")->Divide(hists.at(ax_it->first).at("nunubar").at("cherenkov"));

      nue_nuebar_ratios.at(ax_it->first)["cherenkov"] = (TH1*) hists.at(ax_it->first).at("nue").at("cherenkov")->Clone(UniqueName().c_str());
      nue_nuebar_ratios.at(ax_it->first).at("cherenkov")->Divide(hists.at(ax_it->first).at("nuebar").at("cherenkov"));
      
      // calib shape
      nuebar_fractions.at(ax_it->first)["calib_shape"] = (TH1*) hists.at(ax_it->first).at("nuebar").at("calib_shape")->Clone(UniqueName().c_str());
      nuebar_fractions.at(ax_it->first).at("calib_shape")->Divide(hists.at(ax_it->first).at("nunubar").at("calib_shape"));

      nue_nuebar_ratios.at(ax_it->first)["calib_shape"] = (TH1*) hists.at(ax_it->first).at("nue").at("calib_shape")->Clone(UniqueName().c_str());
      nue_nuebar_ratios.at(ax_it->first).at("calib_shape")->Divide(hists.at(ax_it->first).at("nuebar").at("calib_shape"));
      

      // multiverses
      TH1 * ppfx_up = mv_ppfx_nuebar_fractions.at(ax_it->first)->GetPlusOneSigmaShift (nuebar_fractions.at(ax_it->first).at("nominal"));
      TH1 * ppfx_dw = mv_ppfx_nuebar_fractions.at(ax_it->first)->GetMinusOneSigmaShift(nuebar_fractions.at(ax_it->first).at("nominal"));

      TH1 * genie_up = mv_genie_nuebar_fractions.at(ax_it->first)->GetPlusOneSigmaShift (nuebar_fractions.at(ax_it->first).at("nominal"));
      TH1 * genie_dw = mv_genie_nuebar_fractions.at(ax_it->first)->GetMinusOneSigmaShift(nuebar_fractions.at(ax_it->first).at("nominal"));

      nuebar_fractions[ax_it->first]["ppfx"]  = MostConservative(ppfx_up , ppfx_dw , nuebar_fractions.at(ax_it->first).at("nominal"));
      nuebar_fractions[ax_it->first]["genie"] = MostConservative(genie_up, genie_dw, nuebar_fractions.at(ax_it->first).at("nominal"));


      delete ppfx_up;
      delete ppfx_dw;
      delete genie_up;
      delete genie_dw;
      
      ppfx_up = mv_ppfx_nue_nuebar_ratios.at(ax_it->first)->GetPlusOneSigmaShift (nue_nuebar_ratios.at(ax_it->first).at("nominal"));
      ppfx_dw = mv_ppfx_nue_nuebar_ratios.at(ax_it->first)->GetMinusOneSigmaShift(nue_nuebar_ratios.at(ax_it->first).at("nominal"));
      nue_nuebar_ratios.at(ax_it->first)["ppfx"] = MostConservative(ppfx_up, ppfx_dw, nue_nuebar_ratios.at(ax_it->first).at("nominal"));

      genie_up = mv_genie_nue_nuebar_ratios.at(ax_it->first)->GetPlusOneSigmaShift (nue_nuebar_ratios.at(ax_it->first).at("nominal"));
      genie_dw = mv_genie_nue_nuebar_ratios.at(ax_it->first)->GetMinusOneSigmaShift(nue_nuebar_ratios.at(ax_it->first).at("nominal"));
      nue_nuebar_ratios.at(ax_it->first)["genie"] = MostConservative(genie_up, genie_dw, nue_nuebar_ratios.at(ax_it->first).at("nominal"));
    }

    for(auto it1 = nuebar_fractions.begin(); it1 != nuebar_fractions.end(); it1++)
      for(auto it2 = nuebar_fractions.at(it1->first).begin(); it2 != nuebar_fractions.at(it1->first).end(); it2++)
        it2->second->GetYaxis()->SetTitle("#bar{#nu_{e}} Fraction");
    
    // uncertainties
    std::map<std::string, std::map<std::string, TH1*> > nuebar_fractions_uncert_abs;
    std::map<std::string, std::map<std::string, TH1*> > nuebar_fractions_uncert_frac;
    std::map<std::string, std::map<std::string, TH1*> > nue_nuebar_ratios_uncert_abs;
    std::map<std::string, std::map<std::string, TH1*> > nue_nuebar_ratios_uncert_frac;
    for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
      for(std::string syst : {"ppfx", "genie", "calib_shift", "calib_shape", "ll", "cherenkov"} ) {
        nuebar_fractions_uncert_abs[ax_it->first][syst] = (TH1*) nuebar_fractions.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
        nuebar_fractions_uncert_abs[ax_it->first][syst]->Add(nuebar_fractions.at(ax_it->first).at("nominal"), -1);
        abs(nuebar_fractions_uncert_abs[ax_it->first][syst]);

        nuebar_fractions_uncert_frac[ax_it->first][syst] = (TH1*) nuebar_fractions_uncert_abs.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
        nuebar_fractions_uncert_frac.at(ax_it->first).at(syst)->Divide(nuebar_fractions.at(ax_it->first).at("nominal"));


        // ratios
        nue_nuebar_ratios_uncert_abs[ax_it->first][syst] = (TH1*) nue_nuebar_ratios.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
        nue_nuebar_ratios_uncert_abs[ax_it->first][syst]->Add(nue_nuebar_ratios.at(ax_it->first).at("nominal"), -1);
        abs(nue_nuebar_ratios_uncert_abs[ax_it->first][syst]);

        nue_nuebar_ratios_uncert_frac[ax_it->first][syst] = (TH1*) nue_nuebar_ratios_uncert_abs.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
        nue_nuebar_ratios_uncert_frac.at(ax_it->first).at(syst)->Divide(nue_nuebar_ratios.at(ax_it->first).at("nominal"));
      }

      // beam uncertainties
      for(auto beam_it = beamSystematics.begin(); beam_it != beamSystematics.end(); beam_it++) {
        std::string syst = beam_it->first;
        nuebar_fractions_uncert_abs[ax_it->first][syst] = (TH1*) nuebar_fractions.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
        nuebar_fractions_uncert_abs[ax_it->first][syst]->Add(nuebar_fractions.at(ax_it->first).at("nominal"), -1);
        abs(nuebar_fractions_uncert_abs[ax_it->first][syst]);

        nuebar_fractions_uncert_frac[ax_it->first][syst] = (TH1*) nuebar_fractions_uncert_abs.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
        nuebar_fractions_uncert_frac.at(ax_it->first).at(syst)->Divide(nuebar_fractions.at(ax_it->first).at("nominal"));

        // ratios
        nue_nuebar_ratios_uncert_abs[ax_it->first][syst] = (TH1*) nue_nuebar_ratios.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
        nue_nuebar_ratios_uncert_abs[ax_it->first][syst]->Add(nue_nuebar_ratios.at(ax_it->first).at("nominal"), -1);
        abs(nue_nuebar_ratios_uncert_abs[ax_it->first][syst]);

        nue_nuebar_ratios_uncert_frac[ax_it->first][syst] = (TH1*) nue_nuebar_ratios_uncert_abs.at(ax_it->first).at(syst)->Clone(UniqueName().c_str());
        nue_nuebar_ratios_uncert_frac.at(ax_it->first).at(syst)->Divide(nue_nuebar_ratios.at(ax_it->first).at("nominal"));
      }

      // combine beam uncertainties into one
      nuebar_fractions_uncert_frac.at(ax_it->first)["beam_focusing"] = (TH1*) nuebar_fractions_uncert_frac.at(ax_it->first).at("2kA")->Clone(UniqueName().c_str());
      nuebar_fractions_uncert_frac.at(ax_it->first).at("beam_focusing")->Reset();

      nue_nuebar_ratios_uncert_frac.at(ax_it->first)["beam_focusing"] = (TH1*) nue_nuebar_ratios_uncert_frac.at(ax_it->first).at("2kA")->Clone(UniqueName().c_str());
      nue_nuebar_ratios_uncert_frac.at(ax_it->first).at("beam_focusing")->Reset();

      auto NX = nuebar_fractions_uncert_frac.at(ax_it->first).at("beam_focusing")->GetNbinsX();
      auto NY = nuebar_fractions_uncert_frac.at(ax_it->first).at("beam_focusing")->GetNbinsY();
      auto NZ = nuebar_fractions_uncert_frac.at(ax_it->first).at("beam_focusing")->GetNbinsZ();
      for(auto x = 1; x <= NX; x++) {
        for(auto y = 1; y <= NY; y++) {
          for(auto z = 1; z <= NZ; z++) {
            double rat = 0;
            double val = 0;
            for(auto beam_it = beamSystematics.begin(); beam_it != beamSystematics.end(); beam_it++) {
              rat += std::pow(nue_nuebar_ratios_uncert_frac.at(ax_it->first).at(beam_it->first)->GetBinContent(x, y, z), 2);
              val += std::pow(nuebar_fractions_uncert_frac.at(ax_it->first).at(beam_it->first)->GetBinContent(x, y, z), 2);
            }
            nue_nuebar_ratios_uncert_frac.at(ax_it->first).at("beam_focusing")->SetBinContent(x, y, z, std::sqrt(rat));
            nuebar_fractions_uncert_frac.at(ax_it->first).at("beam_focusing")->SetBinContent(x, y, z, std::sqrt(val));
          }
        }
      }

      // MC statistics
      nuebar_fractions_uncert_abs .at(ax_it->first)["stat"] = (TH1*) nuebar_fractions.at(ax_it->first).at("nominal")->Clone(UniqueName().c_str());
      nuebar_fractions_uncert_frac.at(ax_it->first)["stat"] = (TH1*) nuebar_fractions.at(ax_it->first).at("nominal")->Clone(UniqueName().c_str());
      
      NX = nuebar_fractions.at(ax_it->first).at("nominal")->GetNbinsX();
      NY = nuebar_fractions.at(ax_it->first).at("nominal")->GetNbinsY();
      NZ = nuebar_fractions.at(ax_it->first).at("nominal")->GetNbinsZ();
      for(auto ix = 1; ix <= NX; ix++) {
        for(auto iy = 1; iy <= NY; iy++) {
          for(auto iz = 1; iz <= NZ; iz++) {
            // rescale back to full MC statistics

            auto num = hists.at(ax_it->first).at("nuebar") .at("nominal")->GetBinContent(ix, iy, iz) * MCPOT / kAna2020RHCPOT;
            auto den = hists.at(ax_it->first).at("nunubar").at("nominal")->GetBinContent(ix, iy, iz) * MCPOT / kAna2020RHCPOT;
            auto err = std::sqrt(1/num + 1/den);
            auto nom = num / den;
            nuebar_fractions_uncert_abs .at(ax_it->first).at("stat")->SetBinContent(ix, iy, iz, nom * err);
            nuebar_fractions_uncert_frac.at(ax_it->first).at("stat")->SetBinContent(ix, iy, iz, err      );

          }
        }
      }
    }

    // sum in quadrature
    std::map<std::string, TH1*> nsel_uncert;
    std::map<std::string, TH1*> mc_frac_uncert;
    for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
      nsel_uncert   [ax_it->first] = (TH1*) hists.at(ax_it->first).at("nuebar").at("nominal")->Clone(UniqueName().c_str());
      mc_frac_uncert[ax_it->first] = (TH1*) nuebar_fractions.at(ax_it->first).at("nominal")->Clone(UniqueName().c_str());
      nsel_uncert   .at(ax_it->first)->Reset();
      mc_frac_uncert.at(ax_it->first)->Reset();
      mc_frac_uncert.at(ax_it->first)->SetTitle("#bar{#nu}_{e} Fraction Uncertainty");
      nsel_uncert   .at(ax_it->first)->SetTitle("Selected #bar{#nu}_{e} Uncertainty");

      if(ax_it->first == "electron_kin")
        mc_frac_uncert.at(ax_it->first)->GetZaxis()->SetTitle("Fractional Uncertainty");
      else
        mc_frac_uncert.at(ax_it->first)->GetYaxis()->SetTitle("Fractional Uncertainty");
      
      auto NX = nuebar_fractions.at(ax_it->first).at("nominal")->GetNbinsX();
      auto NY = nuebar_fractions.at(ax_it->first).at("nominal")->GetNbinsY();
      auto NZ = nuebar_fractions.at(ax_it->first).at("nominal")->GetNbinsZ();      

      for(auto ix = 1; ix <= NX; ix++) {
        for(auto iy = 1; iy <= NY; iy++) {
          for(auto iz = 1; iz <= NZ; iz++) {
            double val = 0;
            for(std::string syst : {"ppfx", "genie", "calib_shape", "calib_shift", "ll", "cherenkov", "beam_focusing", "stat"}) {
              val += std::pow(nuebar_fractions_uncert_frac.at(ax_it->first).at(syst)->GetBinContent(ix, iy, iz), 2);
            }// syst
            mc_frac_uncert.at(ax_it->first)->SetBinContent(ix, iy, iz, std::sqrt(val));

            auto n  = hists.at(ax_it->first).at("nunubar").at("nominal")->GetBinContent(ix, iy, iz);
            nsel_uncert.at(ax_it->first)->SetBinContent(ix, iy, iz, std::sqrt(1/n + val));

          } // iz
        } // iy
      } // ix
    } // ax_it

    mc_frac_uncert.at("enu")->GetXaxis()->SetRangeUser(0, 6);
    mc_frac_uncert.at("q2" )->GetXaxis()->SetRangeUser(0, 2);
    nue_nuebar_ratios.at("enu").at("nominal")->GetXaxis()->SetRangeUser(0, 6);
    nue_nuebar_ratios.at("q2").at("nominal") ->GetXaxis()->SetRangeUser(0, 2);
    
    Plot1D(mc_frac_uncert.at("enu"), "Total", nuebar_fractions_uncert_frac.at("enu"), all_systs, plot_dump + "/enu_fractional_uncertainties_all.pdf");
    Plot1D(mc_frac_uncert.at("q2") , "Total", nuebar_fractions_uncert_frac.at("q2") , all_systs, plot_dump + "/q2_fractional_uncertainties_all.pdf" );
    Plot1D(mc_frac_uncert.at("enu"), "Total", nuebar_fractions_uncert_frac.at("enu"), beam_systs, plot_dump + "/enu_fractional_uncertainties_beam.pdf");
    Plot1D(mc_frac_uncert.at("q2") , "Total", nuebar_fractions_uncert_frac.at("q2") , beam_systs, plot_dump + "/q2_fractional_uncertainties_beam.pdf" );
    for(auto ax_it = axes.begin(); ax_it != axes.end(); ax_it++) {
      //      Plot1D(nuebar_fractions .at(ax_it->first).at("nominal"), "Nominal", nuebar_fractions.at(ax_it->first) , all_systs, plot_dump + "/" + ax_it->first + "_nuebar_fraction_all.pdf"  , xmax.at(ax_it->first));
      //      Plot1D(nue_nuebar_ratios.at(ax_it->first).at("nominal"), "Nominal", nue_nuebar_ratios.at(ax_it->first), all_systs, plot_dump + "/" + ax_it->first + "_nue_nuebar_ratios_all.pdf", xmax.at(ax_it->first));
      //
      //      Plot1D(nuebar_fractions .at(ax_it->first).at("nominal"), "Nominal", nuebar_fractions.at(ax_it->first) , beam_systs, plot_dump + "/" + ax_it->first + "_nuebar_fraction_beam.pdf"  , xmax.at(ax_it->first));
      //      Plot1D(nue_nuebar_ratios.at(ax_it->first).at("nominal"), "Nominal", nue_nuebar_ratios.at(ax_it->first), beam_systs, plot_dump + "/" + ax_it->first + "_nue_nuebar_ratios_beam.pdf", xmax.at(ax_it->first));

      //      Plot1D(hists.at(ax_it->first).at("nuebar" ).at("nominal"), "Nominal", hists.at(ax_it->first).at("nuebar" ), plot_dump + "/" + ax_it->first + "_nuebar_events.pdf" , xmax.at(ax_it->first));
      //      Plot1D(hists.at(ax_it->first).at("nue"    ).at("nominal"), "Nominal", hists.at(ax_it->first).at("nue"    ), plot_dump + "/" + ax_it->first + "_nue_events.pdf"    , xmax.at(ax_it->first));
      //      Plot1D(hists.at(ax_it->first).at("nunubar").at("nominal"), "Nominal", hists.at(ax_it->first).at("nunubar"), plot_dump + "/" + ax_it->first + "_nunubar_events.pdf", xmax.at(ax_it->first));

      PlotMultiverse(*mv_ppfx .at(ax_it->first).at("nuebar" ), hists.at(ax_it->first).at("nuebar" ).at("nominal"), "Selected #bar{#nu_{e}} PPFX" , plot_dump + "/" + ax_it->first + "_nuebar_events_ppfx.pdf" );
      PlotMultiverse(*mv_genie.at(ax_it->first).at("nuebar" ), hists.at(ax_it->first).at("nuebar" ).at("nominal"), "Selected #bar{#nu_{e}} GENIE", plot_dump + "/" + ax_it->first + "_nuebar_events_genie.pdf");
      PlotMultiverse(*mv_ppfx .at(ax_it->first).at("nue"    ), hists.at(ax_it->first).at("nue"    ).at("nominal"), "Selected #nu_{e} PPFX" , plot_dump + "/" + ax_it->first + "_nue_events_ppfx.pdf" );
      PlotMultiverse(*mv_genie.at(ax_it->first).at("nue"    ), hists.at(ax_it->first).at("nue"    ).at("nominal"), "Selected #nu_{e} GENIE", plot_dump + "/" + ax_it->first + "_nue_events_genie.pdf");
      PlotMultiverse(*mv_ppfx .at(ax_it->first).at("nunubar"), hists.at(ax_it->first).at("nunubar").at("nominal"), "Selected #bar{#nu_{e}} + #nu_{e} PPFX" , plot_dump + "/" + ax_it->first + "_nunubar_events_ppfx.pdf" );
      PlotMultiverse(*mv_genie.at(ax_it->first).at("nunubar"), hists.at(ax_it->first).at("nunubar").at("nominal"), "Selected #bar{#nu_{e}} + #nu_{e} GENIE", plot_dump + "/" + ax_it->first + "_nunubar_events_genie.pdf");

      PlotMultiverse(*mv_ppfx_nuebar_fractions .at(ax_it->first), nuebar_fractions.at(ax_it->first).at("nominal"), "#bar{#nu_{e}} Fraction PPFX" , plot_dump + "/" + ax_it->first + "_nuebar_fraction_ppfx.pdf",
                     nuebar_fractions.at(ax_it->first).at("ppfx"));
      
      PlotMultiverse(*mv_genie_nuebar_fractions.at(ax_it->first), nuebar_fractions.at(ax_it->first).at("nominal"), "#bar{#nu_{e}} Fraction GENIE", plot_dump + "/" + ax_it->first + "_nuebar_fraction_genie.pdf",
                     nuebar_fractions.at(ax_it->first).at("genie"));

      PlotMultiverse(*mv_ppfx_nue_nuebar_ratios.at(ax_it->first), nue_nuebar_ratios.at(ax_it->first).at("nominal"), "#nu_{e} / #bar{#nu}_{e} PPFX", plot_dump + "/" + ax_it->first + "_nue_nuebar_ratio_ppfx.pdf",
                     nue_nuebar_ratios.at(ax_it->first).at("ppfx"));

      PlotMultiverse(*mv_genie_nue_nuebar_ratios.at(ax_it->first), nue_nuebar_ratios.at(ax_it->first).at("nominal"), "#nu_{e} / #bar{#nu}_{e} GENIE", plot_dump + "/" + ax_it->first + "_nue_nuebar_ratio_genie.pdf",
                     nue_nuebar_ratios.at(ax_it->first).at("genie"));
                     

    }
    gStyle->SetPaintTextFormat(".2f");
    Plot2D(mc_frac_uncert.at("electron_kin"), nuebar_fractions_uncert_frac.at("electron_kin"), all_systs, plot_dump + "/electron_kin_fractional_uncertainty");

    TCanvas * c = new TCanvas();
    nsel_uncert.at("q2")->GetYaxis()->SetTitle("Fractional Uncertainty");
    nsel_uncert.at("q2")->GetXaxis()->SetRangeUser(0, xmax.at("q2"));
    nsel_uncert.at("q2")->Draw("hist");
    c->Print(TString::Format("%s/%s_nsel_uncert.pdf", plot_dump.c_str(), "q2").Data());

    nsel_uncert.at("enu")->GetYaxis()->SetTitle("Fractional Uncertainty");
    nsel_uncert.at("enu")->GetXaxis()->SetRangeUser(0, xmax.at("enu"));
    nsel_uncert.at("enu")->Draw("hist");
    c->Print(TString::Format("%s/%s_nsel_uncert.pdf", plot_dump.c_str(), "enu").Data());


    c->SetRightMargin(0.2);
    auto nsel_uncert_electron_kin_2d = ToTH2Helper(std::unique_ptr<TH1>((TH1*) nsel_uncert.at("electron_kin")->Clone(UniqueName().c_str())),
                                                   nuebarccinc::costhetabins,
                                                   nuebarccinc::eelecbins);
    nsel_uncert_electron_kin_2d->GetZaxis()->SetTitle("Fractional Uncertainty");

    nsel_uncert_electron_kin_2d->GetXaxis()->SetRangeUser(0.85, 1);
    nsel_uncert_electron_kin_2d->GetYaxis()->SetRangeUser(0, 6);
    nsel_uncert_electron_kin_2d->GetZaxis()->SetRangeUser(0, 1);
    nsel_uncert_electron_kin_2d->SetMarkerColor(kWhite);
    nsel_uncert_electron_kin_2d->GetXaxis()->SetTitle("True cos#theta_{e}");
    nsel_uncert_electron_kin_2d->GetYaxis()->SetTitle("True Electron Energy (GeV)");
    MovePalette(nsel_uncert_electron_kin_2d);
    nsel_uncert_electron_kin_2d->Draw("colz text");

    c->Print(TString::Format("%s/%s_nsel_uncert.pdf", plot_dump.c_str(), "electron_kin").Data());

    auto nnunubar_2d = ToTH2Helper(std::unique_ptr<TH1>((TH1*) hists.at("electron_kin").at("nunubar").at("nominal")->Clone(UniqueName().c_str())),
                                   nuebarccinc::costhetabins,
                                   nuebarccinc::eelecbins);

    nnunubar_2d->SetTitle("Selected N_{#bar{#nu}_{e}} + N_{#nu_{e}}");
    nnunubar_2d->GetXaxis()->SetRangeUser(0.85, 1);
    nnunubar_2d->GetYaxis()->SetRangeUser(0, 6);
    nnunubar_2d->SetMarkerColor(kWhite);
    nnunubar_2d->GetXaxis()->SetTitle("True cos#theta_{e}");
    nnunubar_2d->GetYaxis()->SetTitle("True Electron Energy (GeV)");
    MovePalette(nnunubar_2d);
    nnunubar_2d->Draw("colz text");

    c->Print(TString::Format("%s/%s_nsel_nunubar.pdf", plot_dump.c_str(), "electron_kin").Data());

    auto nom_nuebar_frac_2d = ToTH2Helper(std::unique_ptr<TH1>((TH1*) nuebar_fractions.at("electron_kin").at("nominal")->Clone(UniqueName().c_str())),
                                          nuebarccinc::costhetabins,
                                          nuebarccinc::eelecbins);

    nom_nuebar_frac_2d->SetTitle("s_{#bar{#nu}_{e}}");
    nom_nuebar_frac_2d->GetXaxis()->SetRangeUser(0.85, 1);
    nom_nuebar_frac_2d->GetYaxis()->SetRangeUser(0, 6);
    nom_nuebar_frac_2d->GetZaxis()->SetRangeUser(0, 1);
    nom_nuebar_frac_2d->SetMarkerColor(kWhite);
    nom_nuebar_frac_2d->GetXaxis()->SetTitle("True cos#theta_{e}");
    nom_nuebar_frac_2d->GetYaxis()->SetTitle("True Electron Energy (GeV)");
    MovePalette(nom_nuebar_frac_2d);
    nom_nuebar_frac_2d->Draw("colz text");

    c->Print(TString::Format("%s/%s_nuebar_fraction_2d.pdf", plot_dump.c_str(), "electron_kin").Data());



    // calculate required MC POT that would make MC stat. uncert
    // negligible -- ~1%

    std::map<std::string, TH1*> req_pot_1d;
    std::map<std::string, TH1*> req_pot_2d;
    std::map<std::string, TH1*> req_pot_wrt_prod5;
    auto target_uncert = 0.01;
    for(auto loader_it = loaders.begin(); loader_it != loaders.end(); loader_it++) {
      req_pot_1d[loader_it->first] = (TH1*) nuebar_fractions.at("electron_kin").at(loader_it->first)->Clone(UniqueName().c_str());
      req_pot_1d.at(loader_it->first)->Reset();

      for(auto i = 1; i <= req_pot_1d.at(loader_it->first)->GetNbinsX(); i++) {
        double s   = nuebar_fractions.at("electron_kin").at(loader_it->first)->GetBinContent(i);
        double n   = hists.at("electron_kin").at("nunubar").at(loader_it->first)->GetBinContent(i);
        double pot = pots.at(loader_it->first);

        if(s != 0. && s != 1)
          req_pot_1d.at(loader_it->first)->SetBinContent(i, (1+s)/s * std::pow(target_uncert, -2) * pot / n);
      }
      
      
      req_pot_2d[loader_it->first] = ToTH2Helper(std::unique_ptr<TH1>((TH1*) req_pot_1d.at(loader_it->first)->Clone(UniqueName().c_str())),
                                                 nuebarccinc::costhetabins,
                                                 nuebarccinc::eelecbins);
      req_pot_wrt_prod5[loader_it->first] = (TH1*) req_pot_2d.at(loader_it->first)->Clone(UniqueName().c_str());
      req_pot_wrt_prod5.at(loader_it->first)->Scale(1./pots.at(loader_it->first));


      gStyle->SetPaintTextFormat("1.1g");
      gStyle->SetPalette(kDeepSea);
      c->SetLogz();
      req_pot_2d.at(loader_it->first)->SetTitle(TString::Format("Required MC POT (%.1f%% Stat. Uncert)", target_uncert * 100));
      req_pot_2d.at(loader_it->first)->GetXaxis()->SetRangeUser(0.85, 1);
      req_pot_2d.at(loader_it->first)->GetYaxis()->SetRangeUser(0, 6);
      req_pot_2d.at(loader_it->first)->SetMarkerColor(kWhite);
      req_pot_2d.at(loader_it->first)->GetXaxis()->SetTitle("True cos#theta_{e}");
      req_pot_2d.at(loader_it->first)->GetYaxis()->SetTitle("True Electron Energy (GeV)");
      req_pot_2d.at(loader_it->first)->GetZaxis()->SetRangeUser(1e15, 1e25);
      MovePalette(req_pot_2d.at(loader_it->first));
      req_pot_2d.at(loader_it->first)->Draw("colz text");


      c->Print(TString::Format("%s/req_pot_%s.pdf", plot_dump.c_str(), loader_it->first.c_str()).Data());

      req_pot_wrt_prod5.at(loader_it->first)->SetTitle("Required MC POT / Prod5 POT");
      req_pot_wrt_prod5.at(loader_it->first)->GetXaxis()->SetRangeUser(0.85, 1);
      req_pot_wrt_prod5.at(loader_it->first)->GetYaxis()->SetRangeUser(0, 6);
      req_pot_wrt_prod5.at(loader_it->first)->SetMaximum(500);
      req_pot_wrt_prod5.at(loader_it->first)->SetMarkerColor(kWhite);
      req_pot_wrt_prod5.at(loader_it->first)->GetXaxis()->SetTitle("True cos#theta_{e}");
      req_pot_wrt_prod5.at(loader_it->first)->GetYaxis()->SetTitle("True Electron Energy (GeV)");
      MovePalette(req_pot_wrt_prod5.at(loader_it->first));
      req_pot_wrt_prod5.at(loader_it->first)->Draw("colz text");


      gStyle->SetPaintTextFormat(".1f");
      c->SetLogz(false);
      c->Print(TString::Format("%s/req_pot_wrt_prod5%s.pdf", plot_dump.c_str(), loader_it->first.c_str()).Data());
    }

    delete c;
    c   = new TCanvas();
    TLegend * leg = 0;
    auto kNuebarColor = kGreen;
    auto kNueColor = kMagenta;    
    // 1D stacked nnunubar plots
    for(auto loader_it = loaders.begin(); loader_it != loaders.end(); loader_it++) {
      for(std::string axis : {"enu", "q2"}) {
        if(leg) 
          delete leg;
        leg = new TLegend(0.6, 0.89, 0.8, 0.7);

        THStack * stack = new THStack(UniqueName().c_str(), "");
        hists.at(axis).at("nue")   .at(loader_it->first)->SetFillColor(kNueColor);
        hists.at(axis).at("nuebar").at(loader_it->first)->SetFillColor(kNuebarColor);
        hists.at(axis).at("nue")   .at(loader_it->first)->SetLineColor(kBlack);
        hists.at(axis).at("nuebar").at(loader_it->first)->SetLineColor(kBlack);
        hists.at(axis).at("nue")   .at(loader_it->first)->SetLineStyle(1);
        hists.at(axis).at("nuebar").at(loader_it->first)->SetLineStyle(1);

        stack->Add(hists.at(axis).at("nue"   ).at(loader_it->first));
        stack->Add(hists.at(axis).at("nuebar").at(loader_it->first));
        leg->AddEntry(hists.at(axis).at("nue"   ).at(loader_it->first), "#nu_{e}"      , "f");
        leg->AddEntry(hists.at(axis).at("nuebar").at(loader_it->first), "#bar{#nu}_{e}", "f");
        
        
        stack->Draw("hist");
        stack->GetXaxis()->SetRangeUser(0, xmax.at(axis));
        stack->GetXaxis()->SetTitle(hists.at(axis).at("nue").at(loader_it->first)->GetXaxis()->GetTitle());
        stack->GetYaxis()->SetTitle("Events");
        leg->Draw("same");
        c->Print(TString::Format("%s/%s_nunubar_stacked_%s.pdf", 
                                 plot_dump.c_str(),
                                 axis.c_str(),
                                 loader_it->first.c_str()).Data());
        
      }
    }
    

  }
}

void abs(TH1* hist)
{
  for(auto x = 1; x <= hist->GetNbinsX(); x++) {
    for(auto y = 1; y <= hist->GetNbinsY(); y++) {
      for(auto z = 1; z <= hist->GetNbinsZ(); z++) {
        auto val = hist->GetBinContent(x, y, z);
        hist->SetBinContent(x, y, z, std::abs(val));
      }
    }
  }
}

void PlotMultiverse(Multiverse & mv, TH1 * nominal, std::string title, std::string name, TH1 * error)
{
  auto hists = mv.GetUniversesHist(kCyan);

  double maxy = -1000;
  for(auto ihist = 0u; ihist < hists.size(); ihist++)
    maxy = std::max(maxy, hists[ihist]->GetMaximum());

  TCanvas * c = new TCanvas();
  nominal->GetYaxis()->SetRangeUser(0, 1.2 * maxy);
  nominal->SetLineColor(kBlue);
  nominal->SetTitle(title.c_str());

  nominal->Draw("hist");

  for(auto ihist = 0u; ihist < hists.size(); ihist++) 
    hists[ihist]->Draw("hist same");
  

  if(error) {
    for(auto x = 1; x <= nominal->GetNbinsX(); x++) {
      nominal->SetBinError(x, std::abs(error->GetBinContent(x) - nominal->GetBinContent(x)));
    }
    nominal->Draw("hist same e1");
  }
  else 
    nominal->Draw("hist same");


  c->Print(name.c_str());
}

void Plot2D(TH1* total, std::map<std::string, TH1*> systs, std::vector<std::string> to_plot, std::string basename)
{

  TCanvas * c = new TCanvas();
  c->SetRightMargin(0.2);
  std::map<std::string, TH1*> hists2d;
  hists2d["total"] = ToTH2Helper(std::unique_ptr<TH1>((TH1*) total->Clone(UniqueName().c_str())),
                                 nuebarccinc::costhetabins,
                                 nuebarccinc::eelecbins);
  for(auto syst : to_plot) {
    hists2d[syst] = ToTH2Helper(std::unique_ptr<TH1>((TH1*) systs.at(syst)->Clone(UniqueName().c_str())),
                                nuebarccinc::costhetabins,
                                nuebarccinc::eelecbins);
  }

  for(auto hist_it = hists2d.begin(); hist_it != hists2d.end(); hist_it++) {
    hist_it->second->SetTitle("");
    hist_it->second->GetXaxis()->SetTitle("True cos#theta_{e}");
    hist_it->second->GetYaxis()->SetTitle("True Electron Energy (GeV)");
    hist_it->second->GetXaxis()->SetRangeUser(0.85, 1);
    hist_it->second->GetYaxis()->SetRangeUser(0, 6);
    hist_it->second->SetMaximum(1);
    hist_it->second->SetMarkerColor(kWhite);
    MovePalette(hist_it->second);
    hist_it->second->Draw("colz text");
    c->Print((basename + "_" + hist_it->first + ".pdf").c_str());
  }
  
}

void Plot1D(TH1* h1, std::string h1_label, std::map<std::string, TH1*> systs, std::vector<std::string> to_plot, std::string name, double xmax) {
  TCanvas * c = new TCanvas();
  TLegend * leg = new TLegend();    

  // find max y
  double maxy = -999;
  for(auto syst : to_plot) 
    maxy = std::max(maxy, systs.at(syst)->GetMaximum());
  
  maxy = std::max(maxy, h1->GetMaximum());


  if(xmax > 0)
    h1->GetXaxis()->SetRangeUser(0, xmax);

  h1->GetYaxis()->SetRangeUser(0, maxy * 1.2);
  h1->SetLineColor(kBlack);
  h1->Draw("hist");
  leg->AddEntry(h1, h1_label.c_str(), "l");
  
  for(auto syst : to_plot ) {
    if(syst == "nominal") continue;
    systs.at(syst)->SetLineColor(syst_styles.at(syst).color);
    systs.at(syst)->SetLineStyle(syst_styles.at(syst).line_style);
    leg->AddEntry(systs.at(syst), syst.c_str(), "l");
    systs.at(syst)->Draw("hist same");
  }
  h1->Draw("hist same");
  leg->Draw();
  c->Print(name.c_str());
}
TH1 * MostConservative(TH1 * up, TH1 * down, TH1 * nominal)
{
  auto ret = (TH1*) up->Clone(UniqueName().c_str());
  auto NX = up->GetNbinsX();
  auto NY = up->GetNbinsY();
  auto NZ = up->GetNbinsZ();
  for(auto x = 1; x <= NX; x++) {
    for(auto y = 1; y <= NY; y++) {
      for(auto z = 1; z <= NZ; z++) {
        auto nom = nominal->GetBinContent(x, y, z);
        auto shift_up   = std::abs(nom - up->GetBinContent(x, y, z));
        auto shift_down = std::abs(nom - down->GetBinContent(x, y, z));

        ret->SetBinContent(x, y, z, std::max(shift_up, shift_down) + nom);
      }
    }
  }
  return ret;
}

void MovePalette(TH1 * hist)
{
  hist->Draw("colz");
  gPad->Update();
  auto palette = (TPaletteAxis*) hist->GetListOfFunctions()->FindObject("palette");
  palette->SetX1NDC(0.81);
  palette->SetX2NDC(0.85);
}