CutOptimization.cxx

Go to the documentation of this file.

#include "CAFAna/XSec/CutOptimization.h"

#include "TVectorD.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TFile.h"

#include <iostream>

namespace ana
{
  void 
  CutOptimization::PlotDebug(TH1 * hist, 
                             std::string draw_option, 
                             std::string title, 
                             std::string name) {
    TDirectory * tmp = gDirectory;
    TCanvas * c = new TCanvas();
    c->SetLeftMargin(0.15);
    c->SetRightMargin(0.01);


    c->cd();

    hist->Draw(draw_option.c_str());
    hist->SetTitle(title.c_str());
    c->Print(name.c_str());
      
    tmp->cd();
  }

  void 
  CutOptimization::PlotDebug(std::vector<TH1 *> hists,
                             std::vector<std::string> labels,
                             std::string draw_option, 
                             std::string title, 
                             std::string name, 
                             double ymax) {
    TDirectory * tmp = gDirectory;
    TCanvas * c = new TCanvas();
    c->SetLeftMargin(0.15);
    c->SetRightMargin(0.01);

    c->cd();
    int colors[] =  {kBlue, kViolet, kMagenta,
                     kPink, kGreen, kTeal,
                     kCyan, kAzure, kGray,
                     kGray+1, kGray+3};

    TLegend * leg = new TLegend(0.4,0.7,0.7,0.9);

    double max_val = -1e9;
    for(auto ihist = 0u; ihist < hists.size(); ihist++) {
      if(hists[ihist] == NULL) continue; 
      max_val = std::max(max_val, hists[ihist]->GetMaximum());
    }

    int drawn = 0;
    for(auto ihist = 0u; ihist < hists.size(); ihist++) {
      // pointer could be null
      if(hists[ihist] == NULL) continue;         

        // if valid pointer, style hist
      hists[ihist]->SetLineColor(colors[ihist]);
      hists[ihist]->SetTitle(title.c_str());

      leg->AddEntry(hists[ihist], labels[ihist].c_str(), "l");
      // use the axes of first valid histogram 
      if(drawn == 0) {      
        // adjust range so all hists will show up       
        if(ymax < 0)
          hists[ihist]->GetYaxis()->SetRangeUser(0, max_val * 1.2);
        else
          hists[ihist]->GetYaxis()->SetRangeUser(0, ymax);
        hists[ihist]->Draw(draw_option.c_str());        
      }
      else {
        hists[ihist]->Draw((draw_option + " same").c_str());    
      }
      drawn++;
    }
    
    leg->Draw("same");
    c->Print(name.c_str());
      
    tmp->cd();
  }

  void 
  CutOptimization::PlotDebug(const std::vector<TH1 *> universes,
                             Spectrum * nominal,
                             std::string title, 
                             std::string name)
  {
    TH1 * h_nominal = nominal->ToTH1(nominal->POT());
    PlotDebug(universes, h_nominal, title, name);
  }

  void 
  CutOptimization::PlotDebug(const std::vector<TH1 *> universes,
                             TH1 * h_nominal,
                             std::string title, 
                             std::string name)
  {
    TDirectory * tmp = gDirectory;
    TCanvas * c = new TCanvas();
    c->SetLeftMargin(0.15);
    c->SetRightMargin(0.01);
    c->cd();

    double max_val = -1e9;
    for(auto ihist = 0u; ihist < universes.size(); ihist++) {
      max_val = std::max(max_val, universes[ihist]->GetMaximum());
    }
    h_nominal->GetYaxis()->SetRangeUser(0, max_val * 1.2);
    h_nominal->SetTitle(title.c_str());
    h_nominal->Draw("hist");

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

    h_nominal->Draw("hist same");


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

    tmp->cd();
  }

  CutOptimization::CutOptimization(const HistAxis * opti,
                                   const Cut * signal_cut,
                                   const Cut * selection_cut)
    : fSignalCut(signal_cut), fSelectionCut(selection_cut), fOptiHistAxis(opti)
  {}

  void
  CutOptimization::DefineNominal(SystematicDef * nominal)
  {
    fNominalSystDef = nominal;

    fNominalSignal         = fNominalSystDef->BuildSpectrum(*fSignalCut);
    fNominalBkgd           = fNominalSystDef->BuildSpectrum(!*fSignalCut);
    fNominalSelectedSignal = fNominalSystDef->BuildSpectrum(*fSelectionCut && *fSignalCut);
    fNominalSelectedBkgd   = fNominalSystDef->BuildSpectrum(*fSelectionCut && !*fSignalCut);
    fNominalSelected       = fNominalSystDef->BuildSpectrum(*fSelectionCut);
  }
  

  void
  CutOptimization::DefineSystematic(SystematicDef * syst)
  {
    fSystDefs.push_back(syst);

    // for one sided systs, store in .up
    fSignal        [syst].up = syst->BuildSpectrumUp(*fSignalCut);
    fBkgd          [syst].up = syst->BuildSpectrumUp(!*fSignalCut);
    fSelectedSignal[syst].up = syst->BuildSpectrumUp(*fSelectionCut && *fSignalCut);
    fSelectedBkgd  [syst].up = syst->BuildSpectrumUp(*fSelectionCut && !*fSignalCut);
    fSelected      [syst].up = syst->BuildSpectrumUp(*fSelectionCut);    

    if(syst->IsTwoSided()) {
      fSignal        [syst].down = syst->BuildSpectrumDown(*fSignalCut);
      fBkgd          [syst].down = syst->BuildSpectrumDown(!*fSignalCut);
      fSelectedSignal[syst].down = syst->BuildSpectrumDown(*fSelectionCut && *fSignalCut);
      fSelectedBkgd  [syst].down = syst->BuildSpectrumDown(*fSelectionCut && !*fSignalCut);
      fSelected      [syst].down = syst->BuildSpectrumDown(*fSelectionCut);    
    }
  }

  void
  CutOptimization::SetSpillCuts(const SpillCut & spillcut)
  {
    for(auto iload = 0u; iload < fNominalSystDef->LoadersUp().size(); iload++)
      fNominalSystDef->Loaders()[iload]->SetSpillCut(spillcut);
    
    for(auto isyst = 0u; isyst < fSystDefs.size(); isyst++) {
      for(auto iload = 0u; iload < fSystDefs[isyst]->LoadersUp().size(); iload++)
        fSystDefs[isyst]->LoadersUp()[iload]->SetSpillCut(spillcut);
      if(fSystDefs[isyst]->IsTwoSided()) {
        for(auto iload = 0u; iload < fSystDefs[isyst]->LoadersDown().size(); iload++)
          fSystDefs[isyst]->LoadersDown()[iload]->SetSpillCut(spillcut);
      }
    }
  }


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

  void
  CutOptimization::DefineMultiverseSystematic(SystematicDef * syst,
                                              std::vector<Var> mv_weights)
  {
    fMVSystDefs.push_back(syst);
    fMVSignal[syst] = new Multiverse(syst->Loaders(),
                                      *fOptiHistAxis,
                                      *fSignalCut,
                                      *syst->Shifts(),
                                      mv_weights,
                                      *syst->Weight());

    fMVBkgd[syst] = new Multiverse(syst->Loaders(),
                                    *fOptiHistAxis,
                                    !*fSignalCut,
                                    *syst->Shifts(),
                                    mv_weights,
                                    *syst->Weight());
    fMVSelectedSignal[syst] =  new Multiverse(syst->Loaders(),
                                               *fOptiHistAxis,
                                               *fSelectionCut && *fSignalCut,
                                               *syst->Shifts(),
                                               mv_weights,
                                               *syst->Weight());
    fMVSelectedBkgd[syst] = new Multiverse(syst->Loaders(),
                                           *fOptiHistAxis,
                                           *fSelectionCut && !*fSignalCut,
                                           *syst->Shifts(),
                                           mv_weights,
                                           *syst->Weight());
    fMVSelected[syst] = new Multiverse(syst->Loaders(),
                                        *fOptiHistAxis,
                                        *fSelectionCut,
                                        *syst->Shifts(),
                                        mv_weights,
                                        *syst->Weight());                        
  }

  //////////////////////////////////////////////////////////////
  void
  CutOptimization::DefineMultiverseSystematic(SystematicDef * syst, 
                                              std::vector<SystShifts> mv_shifts) {
    fMVSystDefs.push_back(syst);
    fMVSignal[syst] = new Multiverse(syst->Loaders(),
                                     *fOptiHistAxis,
                                     *fSignalCut,
                                     mv_shifts,
                                     *syst->Weight());
    fMVBkgd[syst] = new Multiverse(syst->Loaders(),
                                   *fOptiHistAxis,
                                   !*fSignalCut,
                                   mv_shifts,
                                   *syst->Weight());
    fMVSelectedSignal[syst] = new Multiverse(syst->Loaders(),
                                             *fOptiHistAxis,
                                             *fSelectionCut && *fSignalCut,
                                             mv_shifts,
                                             *syst->Weight());
    fMVSelectedBkgd[syst] = new Multiverse(syst->Loaders(),
                                           *fOptiHistAxis,
                                           *fSelectionCut && !*fSignalCut,
                                           mv_shifts,
                                           *syst->Weight());
    fMVSelected[syst] = new Multiverse(syst->Loaders(),
                                       *fOptiHistAxis,
                                       *fSelectionCut,
                                       mv_shifts,
                                       *syst->Weight());                         
  }

  //////////////////////////////////////////////////////////////
  void
  CutOptimization::Go()
  {
    fNominalSystDef->Go();
    for(auto isyst = 0u; isyst < fSystDefs.size(); isyst++) 
      fSystDefs[isyst]->Go();
    for(auto isyst = 0u; isyst < fMVSystDefs.size(); isyst++)
      fMVSystDefs[isyst]->Go();
  }


  //////////////////////////////////////////////////////////////
  void 
  CutOptimization::SaveTo(TDirectory * dir, const std::string& name) const
  {
    TDirectory * tmp = gDirectory;

    dir = dir->mkdir(name.c_str()); // switch to subdir
    dir->cd();

    TObjString("CutOptimization").Write("type");

    TDirectory * nom_dir = dir->mkdir("Nominal");
    // save nominal
    fNominalSignal        ->SaveTo(nom_dir, "Signal");
    fNominalBkgd          ->SaveTo(nom_dir, "Bkgd");
    fNominalSelectedSignal->SaveTo(nom_dir, "SelectedSignal");
    fNominalSelectedBkgd  ->SaveTo(nom_dir, "SelectedBkgd");
    fNominalSelected      ->SaveTo(nom_dir, "Selected");

    dir->cd();

    TVectorD nsysts(1);
    nsysts[0] = fSystDefs.size();
    nsysts.Write("nsysts");

    // for each syst def, save up down pairs
    // Iterating over the vector insures systematics will be saved in the same
    // order everytime. Critical for hadd_cafana
    for(auto isyst = 0u; isyst < fSystDefs.size(); isyst++) {
      // save systdefs
      fSystDefs[isyst]->SaveTo(dir, TString::Format("syst_defs/%d", isyst).Data());
     
      TDirectory * syst_dir = dir->mkdir(TString::Format("syst_%s",fSystDefs[isyst]->GetName().c_str()).Data());
      SaveToUpDownSpectra( fSignal        .at(fSystDefs[isyst]), syst_dir->mkdir("fSignal"        ));
      SaveToUpDownSpectra( fBkgd          .at(fSystDefs[isyst]), syst_dir->mkdir("fBkgd"          ));
      SaveToUpDownSpectra( fSelectedSignal.at(fSystDefs[isyst]), syst_dir->mkdir("fSelectedSignal"));
      SaveToUpDownSpectra( fSelectedBkgd  .at(fSystDefs[isyst]), syst_dir->mkdir("fSelectedBkgd"  ));
      SaveToUpDownSpectra( fSelected      .at(fSystDefs[isyst]), syst_dir->mkdir("fSelected"      ));
    }
    

    TVectorD nmultiverse(1);
    nsysts[0] = fMVSystDefs.size();
    nsysts.Write("nmultiverse");

    // for each mv syst def, save multiverse
    // Iterating over the vector insures multiverses will be saved in the same
    // order everytime. Critical for hadd_cafana
    for(auto imv = 0u; imv < fMVSystDefs.size(); imv++) {
      // save mv syst_defs
      fMVSystDefs[imv]->SaveTo(dir, TString::Format("mv_defs/%d", imv).Data());

      TDirectory * mv_dir = dir->mkdir(TString::Format("mv_%s",fMVSystDefs[imv]->GetName().c_str()).Data());
      fMVSignal        .at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVSignal"        );
      fMVBkgd          .at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVBkgd"          );
      fMVSelectedSignal.at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVSelectedSignal");
      fMVSelectedBkgd  .at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVSelectedBkgd"  );
      fMVSelected      .at(fMVSystDefs[imv])->SaveTo(mv_dir, "fMVSelected"      );
    }

    dir->Write();
    delete dir;

    tmp->cd();

  }
  //////////////////////////////////////////////////////////////
  std::unique_ptr<CutOptimization> 
  CutOptimization::LoadFrom(TDirectory * dir, const std::string& name)
  {
    dir = dir->GetDirectory(name.c_str()); // switch to subdir
    assert(dir);

    Spectrum * nominal_signal;
    Spectrum * nominal_bkgd;
    Spectrum * nominal_selected_signal;
    Spectrum * nominal_selected_bkgd;
    Spectrum * nominal_selected;

    std::vector<SystematicDef *> syst_defs;
    std::map<SystematicDef *, UpDownPair<Spectrum> > syst_signal;
    std::map<SystematicDef *, UpDownPair<Spectrum> > syst_bkgd;
    std::map<SystematicDef *, UpDownPair<Spectrum> > syst_selected_signal;
    std::map<SystematicDef *, UpDownPair<Spectrum> > syst_selected_bkgd;
    std::map<SystematicDef *, UpDownPair<Spectrum> > syst_selected;

    std::vector<SystematicDef *> mv_defs;
    std::map<SystematicDef *, Multiverse * > mv_signal;
    std::map<SystematicDef *, Multiverse * > mv_bkgd;
    std::map<SystematicDef *, Multiverse * > mv_selected_signal;
    std::map<SystematicDef *, Multiverse * > mv_selected_bkgd;
    std::map<SystematicDef *, Multiverse * > mv_selected;

    TDirectory * nom_dir = dir->GetDirectory("Nominal");
    nominal_signal          = Spectrum::LoadFrom(nom_dir, "Signal"        ).release();
    nominal_bkgd            = Spectrum::LoadFrom(nom_dir, "Bkgd"          ).release();
    nominal_selected_signal = Spectrum::LoadFrom(nom_dir, "SelectedSignal").release();
    nominal_selected_bkgd   = Spectrum::LoadFrom(nom_dir, "SelectedBkgd"  ).release();
    nominal_selected        = Spectrum::LoadFrom(nom_dir, "Selected"      ).release();
    
    TVectorD * vsysts = (TVectorD*) dir->Get("nsysts");
    int nsysts = (*vsysts)[0];
    for(int isyst = 0; isyst < nsysts; isyst++) {
      syst_defs.push_back(SystematicDef::LoadFrom(dir, TString::Format("syst_defs/%d", isyst).Data()).release());

      TDirectory * syst_dir = dir->GetDirectory(TString::Format("syst_%s", syst_defs.back()->GetName().c_str()));      
      syst_signal         [syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fSignal"         ));
      syst_bkgd           [syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fBkgd"           ));
      syst_selected_signal[syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fSelectedSignal" ));
      syst_selected_bkgd  [syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fSelectedBkgd"   ));
      syst_selected       [syst_defs.back()] = LoadFromUpDownSpectra(syst_dir->GetDirectory("fSelected"       ));
    }
    
    TVectorD * vmultiverse = (TVectorD*) dir->Get("nmultiverse");
    int nmultiverse = (*vmultiverse)[0];
    for(int imv = 0; imv < nmultiverse; imv++) {
      mv_defs.push_back(SystematicDef::LoadFrom(dir, TString::Format("mv_defs/%d", imv).Data()).release());

      TDirectory * mv_dir = dir->GetDirectory(TString::Format("mv_%s", mv_defs.back()->GetName().c_str()));
      mv_signal         [mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVSignal"         ).release();
      mv_bkgd           [mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVBkgd"           ).release();
      mv_selected_signal[mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVSelectedSignal" ).release();
      mv_selected_bkgd  [mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVSelectedBkgd"   ).release();
      mv_selected       [mv_defs.back()] = Multiverse::LoadFrom(mv_dir, "fMVSelected"       ).release();
    }

    delete dir;

    return std::make_unique<CutOptimization> (nominal_signal,
                                              nominal_bkgd,
                                              nominal_selected_signal,
                                              nominal_selected_bkgd,
                                              nominal_selected,
                                              syst_defs,
                                              syst_signal,
                                              syst_bkgd,
                                              syst_selected_signal,
                                              syst_selected_bkgd,
                                              syst_selected,
                                              mv_defs,
                                              mv_signal,
                                              mv_bkgd,
                                              mv_selected_signal,
                                              mv_selected_bkgd,
                                              mv_selected);   
  }

  //////////////////////////////////////////////////////////////
  CutOptimization::CutOptimization(Spectrum * nominal_signal,
                                   Spectrum * nominal_bkgd,
                                   Spectrum * nominal_selected_signal,
                                   Spectrum * nominal_selected_bkgd,
                                   Spectrum * nominal_selected,
                                   std::vector<SystematicDef*> syst_defs,
                                   std::map<SystematicDef*, UpDownPair<Spectrum > > syst_signal,
                                   std::map<SystematicDef*, UpDownPair<Spectrum > > syst_bkgd,
                                   std::map<SystematicDef*, UpDownPair<Spectrum > > syst_selected_signal,
                                   std::map<SystematicDef*, UpDownPair<Spectrum > > syst_selected_bkgd,
                                   std::map<SystematicDef*, UpDownPair<Spectrum > > syst_selected,
                                   std::vector<SystematicDef*> mv_defs,
                                   std::map<SystematicDef*, Multiverse * > mv_signal,
                                   std::map<SystematicDef*, Multiverse * > mv_bkgd,
                                   std::map<SystematicDef*, Multiverse * > mv_selected_signal,
                                   std::map<SystematicDef*, Multiverse * > mv_selected_bkgd,
                                   std::map<SystematicDef*, Multiverse * > mv_selected)
    : fNominalSignal(nominal_signal),
      fNominalBkgd(nominal_bkgd),
      fNominalSelectedSignal(nominal_selected_signal),
      fNominalSelectedBkgd(nominal_selected_bkgd),
      fNominalSelected(nominal_selected),
      fSystDefs(syst_defs),
      fSignal(syst_signal),
      fBkgd(syst_bkgd),
      fSelectedSignal(syst_selected_signal),
      fSelectedBkgd(syst_selected_bkgd),
      fSelected(syst_selected),
      fMVSystDefs(mv_defs),
      fMVSignal(mv_signal),
      fMVBkgd(mv_bkgd),
      fMVSelectedSignal(mv_selected_signal),
      fMVSelectedBkgd(mv_selected_bkgd),
      fMVSelected(mv_selected)
  {}

  //////////////////////////////////////////////////////////////
  CutOptimization::CutOptimization(const CutOptimization & rhs)
  {
    this->fSignalCut = new Cut(*rhs.fSignalCut);
    this->fSelectionCut = new Cut(*rhs.fSignalCut);
    this->fOptiHistAxis = new HistAxis(*rhs.fOptiHistAxis);

    this->fNominalSystDef = new SystematicDef(*rhs.fNominalSystDef);
    this->fNominalSignal         = new Spectrum(*rhs.fNominalSignal        );
    this->fNominalBkgd           = new Spectrum(*rhs.fNominalBkgd          );
    this->fNominalSelectedSignal = new Spectrum(*rhs.fNominalSelectedSignal);
    this->fNominalSelectedBkgd   = new Spectrum(*rhs.fNominalSelectedBkgd  );
    this->fNominalSelected       = new Spectrum(*rhs.fNominalSelected      );
    for(auto isyst = 0u; isyst < rhs.fSystDefs.size(); isyst++) {
      this->fSystDefs.push_back(new SystematicDef(*rhs.fSystDefs[isyst]));
      this->fSignal        [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSignal        .at(rhs.fSystDefs[isyst]));
      this->fBkgd          [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fBkgd          .at(rhs.fSystDefs[isyst]));
      this->fSelectedSignal[this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelectedSignal.at(rhs.fSystDefs[isyst]));
      this->fSelectedBkgd  [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelectedBkgd  .at(rhs.fSystDefs[isyst]));
      this->fSelected      [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelected      .at(rhs.fSystDefs[isyst]));
    }
    for(auto imv = 0u; imv < rhs.fMVSystDefs.size(); imv++)  {
      this->fMVSystDefs.push_back(new SystematicDef(*rhs.fMVSystDefs[imv]));
      this->fMVSignal        [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSignal        .at(rhs.fMVSystDefs[imv]));
      this->fMVBkgd          [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVBkgd          .at(rhs.fMVSystDefs[imv]));
      this->fMVSelectedSignal[this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelectedSignal.at(rhs.fMVSystDefs[imv]));
      this->fMVSelectedBkgd  [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelectedBkgd  .at(rhs.fMVSystDefs[imv]));
      this->fMVSelected      [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelected      .at(rhs.fMVSystDefs[imv]));
    }
  }


  //////////////////////////////////////////////////////////////
  CutOptimization::CutOptimization(CutOptimization && rhs)
  {
    this->fSignalCut    = std::move(rhs.fSignalCut);
    this->fSelectionCut = std::move(rhs.fSignalCut);
    this->fOptiHistAxis = std::move(rhs.fOptiHistAxis);

    this->fNominalSystDef        = std::move(rhs.fNominalSystDef       );
    this->fNominalSignal         = std::move(rhs.fNominalSignal        );
    this->fNominalBkgd           = std::move(rhs.fNominalBkgd          );
    this->fNominalSelectedSignal = std::move(rhs.fNominalSelectedSignal);
    this->fNominalSelectedBkgd   = std::move(rhs.fNominalSelectedBkgd  );
    this->fNominalSelected       = std::move(rhs.fNominalSelected      );

    for(auto isyst = 0u; isyst < rhs.fSystDefs.size(); isyst++) {
      this->fSystDefs.push_back(std::move(rhs.fSystDefs[isyst]));
      this->fSignal        [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSignal        .at(rhs.fSystDefs[isyst])));
      this->fBkgd          [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fBkgd          .at(rhs.fSystDefs[isyst])));
      this->fSelectedSignal[this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelectedSignal.at(rhs.fSystDefs[isyst])));
      this->fSelectedBkgd  [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelectedBkgd  .at(rhs.fSystDefs[isyst])));
      this->fSelected      [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelected      .at(rhs.fSystDefs[isyst])));
    }
    for(auto imv = 0u; imv < rhs.fMVSystDefs.size(); imv++)  {
      this->fMVSystDefs.push_back(std::move(rhs.fMVSystDefs[imv]));
      this->fMVSignal        [this->fMVSystDefs.back()] = std::move(rhs.fMVSignal        .at(rhs.fMVSystDefs[imv]));
      this->fMVBkgd          [this->fMVSystDefs.back()] = std::move(rhs.fMVBkgd          .at(rhs.fMVSystDefs[imv]));
      this->fMVSelectedSignal[this->fMVSystDefs.back()] = std::move(rhs.fMVSelectedSignal.at(rhs.fMVSystDefs[imv]));
      this->fMVSelectedBkgd  [this->fMVSystDefs.back()] = std::move(rhs.fMVSelectedBkgd  .at(rhs.fMVSystDefs[imv]));
      this->fMVSelected      [this->fMVSystDefs.back()] = std::move(rhs.fMVSelected      .at(rhs.fMVSystDefs[imv]));
    }
    
  }

  //////////////////////////////////////////////////////////////
  CutOptimization &
  CutOptimization::operator=(const CutOptimization & rhs)
  {
    if(this == &rhs) return *this;

    this->fSignalCut = new Cut(*rhs.fSignalCut);
    this->fSelectionCut = new Cut(*rhs.fSignalCut);
    this->fOptiHistAxis = new HistAxis(*rhs.fOptiHistAxis);

    this->fNominalSystDef = new SystematicDef(*rhs.fNominalSystDef);
    this->fNominalSignal         = new Spectrum(*rhs.fNominalSignal        );
    this->fNominalBkgd           = new Spectrum(*rhs.fNominalBkgd          );
    this->fNominalSelectedSignal = new Spectrum(*rhs.fNominalSelectedSignal);
    this->fNominalSelectedBkgd   = new Spectrum(*rhs.fNominalSelectedBkgd  );
    this->fNominalSelected       = new Spectrum(*rhs.fNominalSelected      );
    for(auto isyst = 0u; isyst < rhs.fSystDefs.size(); isyst++) {
      this->fSystDefs.push_back(new SystematicDef(*rhs.fSystDefs[isyst]));
      this->fSignal        [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSignal        .at(rhs.fSystDefs[isyst]));
      this->fBkgd          [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fBkgd          .at(rhs.fSystDefs[isyst]));
      this->fSelectedSignal[this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelectedSignal.at(rhs.fSystDefs[isyst]));
      this->fSelectedBkgd  [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelectedBkgd  .at(rhs.fSystDefs[isyst]));
      this->fSelected      [this->fSystDefs.back()] = CopyUpDownSpectrum(rhs.fSelected      .at(rhs.fSystDefs[isyst]));
    }
    for(auto imv = 0u; imv < rhs.fMVSystDefs.size(); imv++)  {
      this->fMVSystDefs.push_back(new SystematicDef(*rhs.fMVSystDefs[imv]));
      this->fMVSignal        [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSignal        .at(rhs.fMVSystDefs[imv]));
      this->fMVBkgd          [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVBkgd          .at(rhs.fMVSystDefs[imv]));
      this->fMVSelectedSignal[this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelectedSignal.at(rhs.fMVSystDefs[imv]));
      this->fMVSelectedBkgd  [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelectedBkgd  .at(rhs.fMVSystDefs[imv]));
      this->fMVSelected      [this->fMVSystDefs.back()] = new Multiverse(*rhs.fMVSelected      .at(rhs.fMVSystDefs[imv]));
    }
    return *this;
  }

  //////////////////////////////////////////////////////////////
  CutOptimization &
  CutOptimization::operator=(CutOptimization && rhs)
  {
    this->fSignalCut    = std::move(rhs.fSignalCut);
    this->fSelectionCut = std::move(rhs.fSignalCut);
    this->fOptiHistAxis = std::move(rhs.fOptiHistAxis);

    this->fNominalSystDef        = std::move(rhs.fNominalSystDef       );
    this->fNominalSignal         = std::move(rhs.fNominalSignal        );
    this->fNominalBkgd           = std::move(rhs.fNominalBkgd          );
    this->fNominalSelectedSignal = std::move(rhs.fNominalSelectedSignal);
    this->fNominalSelectedBkgd   = std::move(rhs.fNominalSelectedBkgd  );
    this->fNominalSelected       = std::move(rhs.fNominalSelected      );

    for(auto isyst = 0u; isyst < rhs.fSystDefs.size(); isyst++) {
      this->fSystDefs.push_back(std::move(rhs.fSystDefs[isyst]));
      this->fSignal        [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSignal        .at(rhs.fSystDefs[isyst])));
      this->fBkgd          [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fBkgd          .at(rhs.fSystDefs[isyst])));
      this->fSelectedSignal[this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelectedSignal.at(rhs.fSystDefs[isyst])));
      this->fSelectedBkgd  [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelectedBkgd  .at(rhs.fSystDefs[isyst])));
      this->fSelected      [this->fSystDefs.back()] = CopyUpDownSpectrum(std::move(rhs.fSelected      .at(rhs.fSystDefs[isyst])));
    }
    for(auto imv = 0u; imv < rhs.fMVSystDefs.size(); imv++)  {
      this->fMVSystDefs.push_back(std::move(rhs.fMVSystDefs[imv]));
      this->fMVSignal        [this->fMVSystDefs.back()] = std::move(rhs.fMVSignal        .at(rhs.fMVSystDefs[imv]));
      this->fMVBkgd          [this->fMVSystDefs.back()] = std::move(rhs.fMVBkgd          .at(rhs.fMVSystDefs[imv]));
      this->fMVSelectedSignal[this->fMVSystDefs.back()] = std::move(rhs.fMVSelectedSignal.at(rhs.fMVSystDefs[imv]));
      this->fMVSelectedBkgd  [this->fMVSystDefs.back()] = std::move(rhs.fMVSelectedBkgd  .at(rhs.fMVSystDefs[imv]));
      this->fMVSelected      [this->fMVSystDefs.back()] = std::move(rhs.fMVSelected      .at(rhs.fMVSystDefs[imv]));
    }
    return *this;
  }

  //////////////////////////////////////////////////////////////
  void 
  CutOptimization::Integrate(TH1 *& result, const TH1 * tmp, OptimizeMethod_t method)
  {
    auto ndims = tmp->GetDimension();
    auto NX = tmp->GetNbinsX();
    auto NY = tmp->GetNbinsY();
    auto NZ = tmp->GetNbinsZ();

    // Useful for a 2-sided cut on a single variable like vertex position and fiducial 
    // constraints.
    // Fills a 2D histogram where 
    //  -- on the x axis is the value of the lower bound cut
    //  -- on the y axis is the value of the upper bound cut
    // bin contents are the integral of the 1D input histogram between the bounds corresponding to 
    // upper bin edges
    if(method == kIntegratedInsideBounds) {
      assert(ndims == 1 && "I don't know how to handle this histogram for kIntegratedInsideBounds");

      // don't leak clone of 1D
      delete result;

      // allow for variable bin width
      auto array =  tmp->GetXaxis()->GetXbins()->GetArray();
      if(array) {
        result = new TH2D(UniqueName().c_str(), "", 
                          NX, tmp->GetXaxis()->GetXbins()->GetArray(),
                          NX, tmp->GetXaxis()->GetXbins()->GetArray());
      }
      else {
        auto minx = tmp->GetXaxis()->GetBinLowEdge(1);
        auto maxx = tmp->GetXaxis()->GetBinLowEdge(NX+1);
        result = new TH2D(UniqueName().c_str(), "", 
                          NX, minx, maxx,
                          NX, minx, maxx);
                          
      }
      result->GetYaxis()->SetTitle(TString::Format("%s %s", tmp->GetXaxis()->GetTitle(), "(Upper)"));
      result->GetXaxis()->SetTitle(TString::Format("%s %s", tmp->GetXaxis()->GetTitle(), "(Lower)"));
      
      for(auto ibinlower = 1; ibinlower <= NX; ibinlower++) {
        for(auto ibinupper = ibinlower; ibinupper <=NX; ibinupper++) {
          double error = 0;
          double integral = tmp->IntegralAndError(ibinlower, ibinupper, error);
          result->SetBinContent(ibinlower, ibinupper, integral);
          result->SetBinContent(ibinlower, ibinupper, integral);
        }
      }
      return;
    }
    else {
      switch(ndims) {
      case 1: 
        {
          for(auto ibinx = 1; ibinx <= NX; ibinx++) {
            auto x_a = method == kIntegratedDown? 0     : ibinx;
            auto x_b = method == kIntegratedDown? ibinx : NX+1 ;
            double error    = 0;
            double integral = tmp->IntegralAndError(x_a, x_b, error);
            result->SetBinContent(ibinx, integral);
            result->SetBinError  (ibinx, error   );
        
          }
          break;
        }
      case 2:
        {
          // this way, the casts are only performed twice
          TH2 *       result_2d = static_cast<TH2*>      (result);
          const TH2 * tmp_2d    = static_cast<const TH2*>(tmp   );
          for(auto ibinx = 1; ibinx <= NX; ibinx++) {
            for(auto ibiny = 1; ibiny <= NY; ibiny++) {
              auto x_a = method == kIntegratedDown? 0     : ibinx;
              auto x_b = method == kIntegratedDown? ibinx : NX+1 ;
              auto y_a = method == kIntegratedDown? 0     : ibiny;
              auto y_b = method == kIntegratedDown? ibiny : NY+1 ;
            
              double error    = 0;
              double integral = tmp_2d->IntegralAndError(x_a, x_b,
                                                         y_a, y_b, error);
              result_2d->SetBinContent(ibinx, ibiny, integral);
              result_2d->SetBinError  (ibinx, ibiny, error   );
            }
          }
          result = static_cast<TH1*>(result_2d);
          break;
        }
      case 3: 
        {
          // this way, the casts are only performed twice
          TH3 *       result_3d = static_cast<TH3*>      (result);
          const TH3 * tmp_3d    = static_cast<const TH3*>(tmp   );
          for(auto ibinx = 1; ibinx <= NX; ibinx++) {
            for(auto ibiny = 1; ibiny <= NY; ibiny++) {
              for(auto ibinz = 1; ibinz <= NZ; ibinz++) {
                auto x_a = method == kIntegratedDown? 0     : ibinx;
                auto x_b = method == kIntegratedDown? ibinx : NX+1 ;
                auto y_a = method == kIntegratedDown? 0     : ibiny;
                auto y_b = method == kIntegratedDown? ibiny : NY+1 ;
                auto z_a = method == kIntegratedDown? 0     : ibinz;
                auto z_b = method == kIntegratedDown? ibinz : NZ+1 ;
            
                double error    = 0;
                double integral = tmp_3d->IntegralAndError(x_a, x_b,
                                                           y_a, y_b, 
                                                           z_a, z_b, error);
                result_3d->SetBinContent(ibinx, ibiny, ibinz, integral);
                result_3d->SetBinError  (ibinx, ibiny, ibinz, error   );
              }
            }
          }
          result = static_cast<TH1*>(result_3d);
          break;
        }
      }
    }
    
  }

  //////////////////////////////////////////////////////////////
  // Multiverse should already be scaled and integrated. 
  // This function only returns the absolute shifts that result from the multiverse
  UpDownPair<TH1> 
  CutOptimization::ToUpDownHist(Multiverse * mv, const TH1 * h_nominal)
  {
    TH1 * h_up   = mv->GetPlusOneSigmaShift (h_nominal);
    TH1 * h_down = mv->GetMinusOneSigmaShift(h_nominal);
    // transform into 1 sided shift as the most conservative shift
    TH1 * ret = (TH1*) h_up->Clone();
    for(auto ibinx = 1; ibinx <= ret->GetNbinsX(); ibinx++) {
      for(auto ibiny = 1; ibiny <= ret->GetNbinsY(); ibiny++) {
        for(auto ibinz = 1; ibinz <= ret->GetNbinsZ(); ibinz++) {
          double up      = h_up      ->GetBinContent(ibinx, ibiny, ibinz);
          double down    = h_down    ->GetBinContent(ibinx, ibiny, ibinz);
          double nominal = h_nominal ->GetBinContent(ibinx, ibiny, ibinz);
          double abs_shift = std::max(std::abs(up - nominal), std::abs(down - nominal));
          ret->SetBinContent(ibinx, ibiny, ibinz, abs_shift + nominal);
        }
      }
    }
    return UpDownPair<TH1>{ret};        
  }


  //////////////////////////////////////////////////////////////
  TH1 * 
  CutOptimization::ToHist(const Spectrum * spec, OptimizeMethod_t method, double POT)
  {
    TH1 * tmp;
    TH1 * ret;
    
    // this will happen for one sided systs
    if(spec == NULL) return NULL;

    // if no POT provided, use the spectrum's POT
    if(POT < 0) POT = spec->POT();

    // Call different function based on dimension
    auto ndims = spec->NDimensions();
    switch(ndims) {
    case 1:
      tmp = spec->ToTH1(POT);
      break;
    case 2:
      tmp = spec->ToTH2(POT);
      break;
    case 3:
      tmp = spec->ToTH3(POT);
      break;
    default:
      std::cout << "Spectrum appears to be empty. Something went terribly wrong here" << std::endl;
      exit(1);
    }

    if(method == kIntegratedUp || method == kIntegratedDown || method == kIntegratedInsideBounds) {
      ret = (TH1*) tmp->Clone();

      // integrated result gets saved back to ret
      Integrate(ret, tmp, method);      
    }
    else {
      // default is kBinByBin, just return the histogram then
      ret = tmp;
    }

    return ret;
  }

  // Efficiency is calculate either bin by bin for each numerator and denominator
  // or integrated to/from bin in numerator divided by the total number of signal events.
  // If using an integration method, get total number of signal events from the 
  // first or last bin of integrated denominator.
  //
  // Calculates binomial errors on the efficiency
  // https://home.fnal.gov/~paterno/images/effic.pdf
  TH1 * 
  CutOptimization::Efficiency(TH1* num, TH1 * denom, OptimizeMethod_t method)
  {
    // numerator should already be integrated appropriately 
    TH1 * eff = (TH1*) num->Clone(UniqueName().c_str());

    // if using an integration method, denom should also be integrated appropriately
    // and if so, scale numerator by total number of signal events, which is found
    // in either the first or last bin of denominator
    if(method == kIntegratedUp || method == kIntegratedDown || method == kIntegratedInsideBounds) {
      auto NX = denom->GetNbinsX();
      auto NY = denom->GetNbinsY();
      auto NZ = denom->GetNbinsZ();
      double integral = 0;
      if(method == kIntegratedUp || method == kIntegratedDown) {
        integral = method == kIntegratedUp ? denom->GetBinContent(1) : denom->GetBinContent(NX, NY, NZ);
      }
      else if(method == kIntegratedInsideBounds) {
        // in this case, the total number of signal events will be
        // where the upper bound is at a max and lower bound is at a minimum
        // y axis is upper bound,
        // x axis is lower bound
        integral = denom->GetBinContent(1, NY, 1);
      }

      eff->Scale(1/integral);
      // compute binomial errors by hand
      for(auto ibinx = 1; ibinx <= NX; ibinx++) {
        for(auto ibiny = 1; ibiny <= NY; ibiny++) {
          for(auto ibinz = 1; ibinz <= NZ; ibinz++) {
            double efficiency = eff->GetBinContent(ibinx, ibiny, ibinz);
            double binom_error = 1/integral * std::sqrt( efficiency * (1 - efficiency / integral) );
            eff->SetBinError(ibinx, ibiny, ibinz, binom_error);
          }     
        }       
      }
    }
    else {
      // else just divide numerator with denominator
      eff->Divide(eff, denom, 1, 1, "B"); // divide with binomial error propagation
    }

    eff->SetMaximum(1.2*eff->GetMaximum());

    return eff;
  }

  //////////////////////////////////////////////////////////////
  UpDownPair<TH1>
  CutOptimization::Efficiency(UpDownPair<TH1> num, UpDownPair<TH1> denom, OptimizeMethod_t method)
  {
    TH1 * up   = Efficiency(num.up, denom.up, method);
    TH1 * down = num.down == NULL? NULL : Efficiency(num.down, denom.down, method);
    return UpDownPair<TH1>{up, down};
  }
  
  Multiverse * 
  CutOptimization::Efficiency(Multiverse * num,
                              Multiverse * denom,
                              OptimizeMethod_t method)
  {
    // calculate uncertainty on efficiency of a multiverse as the 
    // +- 1 sigma band of the efficiency multiverse
    Multiverse * eff = new Multiverse(*num);
    if(method == kIntegratedUp || method == kIntegratedDown || method == kIntegratedInsideBounds) {
      eff->Transform(*denom,
                     [method](TH1 * univ_num, TH1 * univ_denom)
                     {
                       return  CutOptimization::Efficiency(univ_num, univ_denom, method);
                     });
    }
    else {
      eff->Divide(*denom, true);
    }
    return eff;
  }

  //////////////////////////////////////////////////////////////
  void 
  CutOptimization::Integrate(Multiverse * mv, OptimizeMethod_t method) 
  {
    mv->Transform([method](TH1 * univ)
                  {
                    TH1 * result = (TH1*) univ->Clone();
                    CutOptimization::Integrate(result, univ, method);
                    return result;
                  });
  }

  //////////////////////////////////////////////////////////////
  void 
  CutOptimization::OptimizedSigmaOverSigma(OptimizeMethod_t method, double data_pot,
                                           TDirectory * save_dir, std::string plot_dump,
                                           bool debug)
  {
    double mc_pot = fNominalSignal->POT();
    double s = data_pot / mc_pot;
    // transform the data into histograms
    // and integrate if method == kIntegratedUp || method == kIntegradedDown || method == kIntegratedInsideBounds
    TH1 * hNominalSignal         = ToHist(fNominalSignal        , method, mc_pot  ); // scaled to mc for eff   
    TH1 * hNominalSelectedSignal = ToHist(fNominalSelectedSignal, method, mc_pot  ); // scaled to mc for eff
    TH1 * hNominalBkgd           = ToHist(fNominalBkgd          , method, mc_pot  ); // scaled to for mc statistics
    TH1 * hNominalSelectedBkgd   = ToHist(fNominalSelectedBkgd  , method, mc_pot  ); // scaled to for mc statistics 
    TH1 * hNominalSelected       = ToHist(fNominalSelected      , method, data_pot); // scaled to data
    TH1 * hNominalEfficiency     = Efficiency(hNominalSelectedSignal, hNominalSignal, method); // scaled to mc
    
    
    // scale the multiverses and calculate efficiencies
    std::map<SystematicDef*, Multiverse *> fMVEfficiency;
    for(auto mv = fMVSystDefs.begin(); mv != fMVSystDefs.end(); mv++) {      
      fMVSignal        [*mv]->Scale(mc_pot   / mc_pot); // for posterity
      fMVSelectedSignal[*mv]->Scale(mc_pot   / mc_pot);
      fMVBkgd          [*mv]->Scale(mc_pot   / mc_pot);
      fMVSelectedBkgd  [*mv]->Scale(mc_pot   / mc_pot);
      fMVSelected      [*mv]->Scale(data_pot / mc_pot);

      if(method == kIntegratedUp || method == kIntegratedDown || method == kIntegratedInsideBounds) {
        Integrate(fMVSignal        [*mv], method);
        Integrate(fMVSelectedSignal[*mv], method);
        Integrate(fMVBkgd          [*mv], method);
        Integrate(fMVSelectedBkgd  [*mv], method);
        Integrate(fMVSelected      [*mv], method);
      }
      
      fMVEfficiency[*mv] = Efficiency(fMVSelectedSignal[*mv], fMVSignal[*mv], method);
    }
    

    std::vector<UpDownPair<TH1 > > hSystSignal;
    std::vector<UpDownPair<TH1 > > hSystBkgd;
    std::vector<UpDownPair<TH1 > > hSystSelectedSignal;
    std::vector<UpDownPair<TH1 > > hSystSelectedBkgd;
    std::vector<UpDownPair<TH1 > > hSystSelected;
    std::vector<UpDownPair<TH1 > > hSystEfficiency;

    std::vector<SystematicDef*> systs;
    for(auto isyst = 0u; isyst < fSystDefs.size(); isyst++) {
      SystematicDef * syst = fSystDefs[isyst];
      systs.push_back(syst);
      hSystSignal        .push_back({ToHist(fSignal         .at(syst).up, method, mc_pot  ), ToHist(fSignal         .at(syst).down, method, mc_pot  )});
      hSystSelectedSignal.push_back({ToHist(fSelectedSignal .at(syst).up, method, mc_pot  ), ToHist(fSelectedSignal .at(syst).down, method, mc_pot  )});

      hSystBkgd          .push_back({ToHist(fBkgd           .at(syst).up, method, mc_pot  ), ToHist(fBkgd           .at(syst).down, method, mc_pot  )});
      hSystSelectedBkgd  .push_back({ToHist(fSelectedBkgd   .at(syst).up, method, mc_pot  ), ToHist(fSelectedBkgd   .at(syst).down, method, mc_pot  )});
      hSystSelected      .push_back({ToHist(fSelected       .at(syst).up, method, data_pot), ToHist(fSelected       .at(syst).down, method, data_pot)});
      hSystEfficiency    .push_back(Efficiency(hSystSelectedSignal.back(), hSystSignal.back(), method));
    }
    for(auto imv = 0u; imv < fMVSystDefs.size(); imv++) {
      SystematicDef * mv = fMVSystDefs[imv];
      systs.push_back(mv);

      hSystSignal        .push_back(ToUpDownHist(fMVSignal         [mv], hNominalSignal        ));
      hSystSelectedSignal.push_back(ToUpDownHist(fMVSelectedSignal [mv], hNominalSelectedSignal));

      hSystBkgd          .push_back(ToUpDownHist(fMVBkgd           [mv], hNominalBkgd          ));
      hSystSelectedBkgd  .push_back(ToUpDownHist(fMVSelectedBkgd   [mv], hNominalSelectedBkgd  ));
      hSystSelected      .push_back(ToUpDownHist(fMVSelected       [mv], hNominalSelected      ));
      hSystEfficiency    .push_back(ToUpDownHist(fMVEfficiency     [mv], hNominalEfficiency    ));
    }

    TH1 * hFracUncertaintySelected     = (TH1*) hNominalSelected      ->Clone(UniqueName().c_str());
    TH1 * hFracUncertaintyBkgdStat     = (TH1*) hNominalSelectedBkgd  ->Clone(UniqueName().c_str());
    TH1 * hFracUncertaintyBkgdSyst     = (TH1*) hNominalBkgd          ->Clone(UniqueName().c_str());
    TH1 * hFracUncertaintyEfficiency   = (TH1*) hNominalSignal        ->Clone(UniqueName().c_str());
    TH1 * hFracUncertaintyXSec         = (TH1*) hNominalSignal        ->Clone(UniqueName().c_str());

    hFracUncertaintySelected     ->Scale(0);
    hFracUncertaintyBkgdStat     ->Scale(0);
    hFracUncertaintyBkgdSyst     ->Scale(0);
    hFracUncertaintyEfficiency   ->Scale(0);
    hFracUncertaintyXSec         ->Scale(0);

    std::vector<TH1*> vabs_uncert_eff_syst  (systs.size());
    std::vector<TH1*> vabs_uncert_bkgd_syst (systs.size());
    std::vector<TH1*> vfrac_uncert_eff_syst (systs.size());
    std::vector<TH1*> vfrac_uncert_bkgd_syst(systs.size());    
    TH1* sel_minus_bkgd = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
    sel_minus_bkgd->Scale(0);
    for(auto isyst = 0u; isyst < systs.size(); isyst++) {
      vabs_uncert_eff_syst  [isyst] = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
      vabs_uncert_bkgd_syst [isyst] = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
      vfrac_uncert_eff_syst [isyst] = (TH1*) hNominalSelected->Clone(UniqueName().c_str());
      vfrac_uncert_bkgd_syst[isyst] = (TH1*) hNominalSelected->Clone(UniqueName().c_str());    
      
      if(hNominalSelected->GetDimension() == 1) {
        vfrac_uncert_bkgd_syst[isyst] ->GetYaxis()->SetTitle("#frac{s#upoint#delta N^{syst}_{bkgd}}{N_{sel} - s#upoint N_{bkgd}}");
        vfrac_uncert_eff_syst [isyst] ->GetYaxis()->SetTitle("#frac{#delta#varepsilon}{#varepsilon}");
        vabs_uncert_bkgd_syst [isyst] ->GetYaxis()->SetTitle("s#upoint#delta N^{syst}_{bkgd}");
        vabs_uncert_eff_syst  [isyst] ->GetYaxis()->SetTitle("#delta#varepsilon");
      }

      else if(hNominalSelected->GetDimension() == 2) {
        vfrac_uncert_bkgd_syst[isyst] ->GetZaxis()->SetTitle("#frac{s#upoint#delta N^{syst}_{bkgd}}{N_{sel} - s#upoint N_{bkgd}}");
        vfrac_uncert_eff_syst [isyst] ->GetZaxis()->SetTitle("#frac{#delta#varepsilon}{#varepsilon}");
        vabs_uncert_bkgd_syst [isyst] ->GetZaxis()->SetTitle("s#upoint#delta N^{syst}_{bkgd}");
        vabs_uncert_eff_syst  [isyst] ->GetZaxis()->SetTitle("#delta#varepsilon");
      }
      vabs_uncert_eff_syst  [isyst] ->Scale(0);
      vabs_uncert_bkgd_syst [isyst] ->Scale(0);
      vfrac_uncert_eff_syst [isyst] ->Scale(0);
      vfrac_uncert_bkgd_syst[isyst] ->Scale(0);    
    }

    auto NX = hNominalSignal->GetNbinsX();
    auto NY = hNominalSignal->GetNbinsY();
    auto NZ = hNominalSignal->GetNbinsZ();
    for(auto ibinx = 1; ibinx <= NX; ibinx++) {
      for(auto ibiny = 1; ibiny <= NY; ibiny++) {
        for(auto ibinz = 1; ibinz <= NZ; ibinz++) {
          double abs_uncert_bkgd_syst = 0;
          double abs_uncert_eff       = 0;

          double nom_sel  = hNominalSelected     ->GetBinContent(ibinx, ibiny, ibinz); // scaled to data
          double nom_bkgd = hNominalSelectedBkgd ->GetBinContent(ibinx, ibiny, ibinz); // scaled to mc
          double nom_eff  = hNominalEfficiency   ->GetBinContent(ibinx, ibiny, ibinz);  
          double nom_signal = nom_sel - s * nom_bkgd;
          
          sel_minus_bkgd->SetBinContent(ibinx, ibiny, ibinz, nom_signal);
          for(auto isyst = 0u; isyst < systs.size(); isyst++) {
            double syst_bkgd = 0;
            double syst_eff  = 0;
            // if two sided, take largest shift
            if(systs[isyst]->IsTwoSided()) {
              
              double up_bkgd   = hSystSelectedBkgd[isyst].up  ->GetBinContent(ibinx, ibiny, ibinz);
              double down_bkgd = hSystSelectedBkgd[isyst].down->GetBinContent(ibinx, ibiny, ibinz);

              double up_eff    = hSystEfficiency  [isyst].up  ->GetBinContent(ibinx, ibiny, ibinz);
              double down_eff  = hSystEfficiency  [isyst].down->GetBinContent(ibinx, ibiny, ibinz);
              
              syst_bkgd = std::max(std::abs(up_bkgd - nom_bkgd), std::abs(down_bkgd - nom_bkgd));
              syst_eff  = std::max(std::abs(up_eff  - nom_eff ), std::abs(down_eff  - nom_eff ));
            }
            else {
              syst_bkgd = std::abs(hSystSelectedBkgd[isyst].up  ->GetBinContent(ibinx, ibiny, ibinz) - nom_bkgd);
              syst_eff  = std::abs(hSystEfficiency  [isyst].up  ->GetBinContent(ibinx, ibiny, ibinz) - nom_eff );
            }

            // scale for MC statistics
            syst_bkgd = s * syst_bkgd;
          
            abs_uncert_bkgd_syst += std::pow(syst_bkgd, 2);
            abs_uncert_eff       += std::pow(syst_eff , 2);         

            vabs_uncert_eff_syst  [isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_eff );
            vabs_uncert_bkgd_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_bkgd);
            if(nom_signal) 
              vfrac_uncert_bkgd_syst[isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_bkgd / nom_signal);    
            else
              vfrac_uncert_bkgd_syst[isyst]->SetBinContent(ibinx, ibiny, ibinz, 100);    
            if(nom_eff)
              vfrac_uncert_eff_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, syst_eff / nom_eff);
            else
              vfrac_uncert_eff_syst [isyst]->SetBinContent(ibinx, ibiny, ibinz, 100);
          } // end loop over systematics

          // calculate fractional uncertainty on efficiency
          double frac_uncert_sel        = 100;
          double frac_uncert_bkgd_stat  = 100;
          double frac_uncert_bkgd_syst  = 100;
          double frac_uncert_eff        = 100;
          if(nom_signal)  {
            frac_uncert_sel       =     std::sqrt(nom_sel)              / nom_signal;
            frac_uncert_bkgd_stat = s * std::sqrt(nom_bkgd)             / nom_signal;
            frac_uncert_bkgd_syst =     std::sqrt(abs_uncert_bkgd_syst) / nom_signal;
          }
          if(nom_eff)
            frac_uncert_eff       = std::sqrt(abs_uncert_eff)       / nom_eff;    

          // finally calculate FOM
          double frac_uncert_sigma = std::sqrt(std::pow(frac_uncert_sel      , 2) +
                                               std::pow(frac_uncert_bkgd_stat, 2) +
                                               std::pow(frac_uncert_bkgd_syst, 2) +
                                               std::pow(frac_uncert_eff      , 2)  );

          
          // fill histograms
          hFracUncertaintyXSec           ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_sigma    );
          hFracUncertaintySelected       ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_sel      );
          hFracUncertaintyBkgdStat       ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_bkgd_stat);
          hFracUncertaintyBkgdSyst       ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_bkgd_syst);
          hFracUncertaintyEfficiency     ->SetBinContent(ibinx, ibiny, ibinz, frac_uncert_eff      );             
        }// end loop over ibinz
      }// end loop over ibiny
    }// end loop over ibinx


    if(debug) {
      //////////////////////////////////// DEBUG //////////////////////////////////
      // save internal histograms to file for later plotting
      TFile * fdebug = new TFile((plot_dump + "/debug.root").c_str(), "recreate");
      
      hNominalSelectedBkgd->Write("hNominalSelectedBkgd");
      hNominalSelected->Write("hNominalSelected");
      hNominalSelectedSignal->Write("hNominalSelectedSignal");
      hNominalSignal->Write("hNominalSignal");
      hNominalBkgd->Write("hNominalBkgd");
      hNominalEfficiency->Write("hNominalEfficiency");

      hFracUncertaintyXSec       ->Write("hFracUncertaintyXSec"       );
      hFracUncertaintySelected   ->Write("hFracUncertaintySelected"   );
      hFracUncertaintyBkgdStat   ->Write("hFracUncertaintyBkgdStat"   );
      hFracUncertaintyBkgdSyst   ->Write("hFracUncertaintyBkgdSyst"   );
      hFracUncertaintyEfficiency ->Write("hFracUncertaintyEfficiency" );


      for(auto isyst = 0u; isyst < systs.size(); isyst++ ) {
        hSystSelectedBkgd  [isyst].up->Write((systs[isyst]->GetName() + "_hSystSelectedBkgd"   + "_UP").c_str());
        hSystSelected      [isyst].up->Write((systs[isyst]->GetName() + "_hSystSelected"       + "_UP").c_str());
        hSystSelectedSignal[isyst].up->Write((systs[isyst]->GetName() + "_hSystSelectedSignal" + "_UP").c_str());
        hSystBkgd          [isyst].up->Write((systs[isyst]->GetName() + "_hSystBkgd"           + "_UP").c_str());
        hSystSignal        [isyst].up->Write((systs[isyst]->GetName() + "_hSystSignal"         + "_UP").c_str());
        hSystEfficiency    [isyst].up->Write((systs[isyst]->GetName() + "_hSystEfficiency"     + "_UP").c_str());

        if(hSystSelectedBkgd[isyst].down != NULL) {
          hSystSelectedBkgd  [isyst].down->Write((systs[isyst]->GetName() + "_hSystSelectedBkgd"   + "_DOWN").c_str());
          hSystSelected      [isyst].down->Write((systs[isyst]->GetName() + "_hSystSelected"       + "_DOWN").c_str());
          hSystSelectedSignal[isyst].down->Write((systs[isyst]->GetName() + "_hSystSelectedSignal" + "_DOWN").c_str());
          hSystBkgd          [isyst].down->Write((systs[isyst]->GetName() + "_hSystBkgd"           + "_DOWN").c_str());
          hSystSignal        [isyst].down->Write((systs[isyst]->GetName() + "_hSystSignal"         + "_DOWN").c_str());
          hSystEfficiency    [isyst].down->Write((systs[isyst]->GetName() + "_hSystEfficiency"     + "_DOWN").c_str());
        }
        
        vabs_uncert_eff_syst  [isyst]->Write((systs[isyst]->GetName() + "_abs_uncert_eff_syst"  ).c_str());
        vabs_uncert_bkgd_syst [isyst]->Write((systs[isyst]->GetName() + "_abs_uncert_bkgd_syst" ).c_str());
        vfrac_uncert_eff_syst [isyst]->Write((systs[isyst]->GetName() + "_frac_uncert_eff_syst" ).c_str());
        vfrac_uncert_bkgd_syst[isyst]->Write((systs[isyst]->GetName() + "_frac_uncert_bkgd_syst").c_str());


//      debug_syst_sel_bkgd_up       .push_back(hSystSelectedBkgd  [isyst].up);
//      debug_syst_sel_bkgd_down     .push_back(hSystSelectedBkgd  [isyst].down);      
//      debug_syst_sel_up            .push_back(hSystSelected      [isyst].up);
//      debug_syst_sel_down          .push_back(hSystSelected      [isyst].down);      
//      debug_syst_sel_signal_up     .push_back(hSystSelectedSignal[isyst].up);
//      debug_syst_sel_signal_down   .push_back(hSystSelectedSignal[isyst].down);      

//      debug_syst_bkgd_up    .push_back(hSystBkgd[isyst].up);
//      debug_syst_bkgd_down  .push_back(hSystBkgd[isyst].down);      

        //      debug_syst_signal_up    .push_back(hSystSignal[isyst].up);
        //      debug_syst_signal_down  .push_back(hSystSignal[isyst].down);      

        //      debug_syst_eff_up    .push_back(hSystEfficiency[isyst].up);
        //      debug_syst_eff_down  .push_back(hSystEfficiency[isyst].down);      

        //      debug_syst_bkgd_labels.push_back(systs[isyst]->GetName());      
        //      debug_syst_bkgd_labels_no_nominal.push_back(systs[isyst]->GetName());      
      }
    

      /*
      PlotDebug(debug_syst_bkgd_up        , debug_syst_bkgd_labels, "hist", "Background +1#sigma"         , plot_dump + "/debug_syst_bkgd_up.pdf"        );
      PlotDebug(debug_syst_bkgd_down      , debug_syst_bkgd_labels, "hist", "Background -1#sigma"         , plot_dump + "/debug_syst_bkgd_down.pdf"      );
      PlotDebug(debug_syst_sel_up         , debug_syst_bkgd_labels, "hist", "Selected +1#sigma"           , plot_dump + "/debug_syst_sel_up.pdf"         );
      PlotDebug(debug_syst_sel_down       , debug_syst_bkgd_labels, "hist", "Selected -1#sigma"           , plot_dump + "/debug_syst_sel_down.pdf"       );
      PlotDebug(debug_syst_sel_bkgd_up    , debug_syst_bkgd_labels, "hist", "Selected Background +1#sigma", plot_dump + "/debug_syst_sel_bkgd_up.pdf"    );
      PlotDebug(debug_syst_sel_bkgd_down  , debug_syst_bkgd_labels, "hist", "Selected Background -1#sigma", plot_dump + "/debug_syst_sel_bkgd_down.pdf"  );
      PlotDebug(debug_syst_sel_signal_up  , debug_syst_bkgd_labels, "hist", "Selected Signal +1#sigma"    , plot_dump + "/debug_syst_sel_signal_up.pdf"  );
      PlotDebug(debug_syst_sel_signal_down, debug_syst_bkgd_labels, "hist", "Selected Signal -1#sigma"    , plot_dump + "/debug_syst_sel_signal_down.pdf");
      PlotDebug(debug_syst_signal_up      , debug_syst_bkgd_labels, "hist", "Signal +1#sigma"             , plot_dump + "/debug_syst_signal_up.pdf"      );
      PlotDebug(debug_syst_signal_down    , debug_syst_bkgd_labels, "hist", "Signal -1#sigma"             , plot_dump + "/debug_syst_signal_down.pdf"    );
      PlotDebug(debug_syst_eff_up         , debug_syst_bkgd_labels, "hist", "Efficiency +1#sigma"         , plot_dump + "/debug_syst_eff_up.pdf"         );
      PlotDebug(debug_syst_eff_down       , debug_syst_bkgd_labels, "hist", "Efficiency -1#sigma"         , plot_dump + "/debug_syst_eff_down.pdf"       );

      PlotDebug(vabs_uncert_eff_syst  , debug_syst_bkgd_labels_no_nominal, "hist", "Abs. Syst. Uncert Efficiency" , plot_dump + "/debug_syst_abs_uncert_eff.pdf"   );
      PlotDebug(vabs_uncert_bkgd_syst , debug_syst_bkgd_labels_no_nominal, "hist", "Abs. Syst. Uncert Background" , plot_dump + "/debug_syst_abs_uncert_bkgd.pdf"  );
      PlotDebug(vfrac_uncert_eff_syst , debug_syst_bkgd_labels_no_nominal, "hist", "Frac. Syst. Uncert Efficiency", plot_dump + "/debug_syst_frac_uncert_eff.pdf"  , 2);
      PlotDebug(vfrac_uncert_bkgd_syst, debug_syst_bkgd_labels_no_nominal, "hist", "Frac. Syst. Uncert Background", plot_dump + "/debug_syst_frac_uncert_bkgd.pdf" , 2);
      */
    
      for(auto syst : fMVSystDefs) {

        fMVSignal        [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_Signal"        );
        fMVBkgd          [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_Bkgd"          );
        fMVSelected      [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_Selected"      );
        fMVSelectedSignal[syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_SelectedSignal");
        fMVSelectedBkgd  [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_SelectedBkgd"  );
        fMVEfficiency    [syst]->SaveTo(fdebug, "MV_" + syst->GetName() + "_Efficiency"    );

//      const std::vector<TH1*> mv_signal_hists     = fMVSignal        [syst]->GetUniversesHist(kCyan);
//      const std::vector<TH1*> mv_bkgd_hists       = fMVBkgd          [syst]->GetUniversesHist(kCyan);
//      const std::vector<TH1*> mv_sel_hists        = fMVSelected      [syst]->GetUniversesHist(kCyan);
//      const std::vector<TH1*> mv_sel_signal_hists = fMVSelectedSignal[syst]->GetUniversesHist(kCyan);
//      const std::vector<TH1*> mv_sel_bkgd_hists   = fMVSelectedBkgd  [syst]->GetUniversesHist(kCyan);

//      std::string name = syst->GetName();

//      PlotDebug(mv_signal_hists    , hNominalSignal        , name + " Multiverse: Signal"        , plot_dump + "/debug_" + name + "_signal.pdf"    );
//      PlotDebug(mv_bkgd_hists      , hNominalBkgd          , name + " Multiverse: Bkgd"          , plot_dump + "/debug_" + name + "_bkgd.pdf"      );
//      PlotDebug(mv_sel_signal_hists, hNominalSelectedSignal, name + " Multiverse: SelectedSignal", plot_dump + "/debug_" + name + "_sel_signal.pdf");
//      PlotDebug(mv_sel_bkgd_hists  , hNominalSelectedBkgd  , name + " Multiverse: SelectedBkgd"  , plot_dump + "/debug_" + name + "_sel_bkgd.pdf"  );
//      PlotDebug(mv_sel_hists       , hNominalSelected      , name + " Multiverse: Selected"      , plot_dump + "/debug_" + name + "_sel.pdf"       );
      }
      
      /*
      PlotDebug(hNominalSignal         , "hist", "hNominalSignal"        , plot_dump + "/debug_nom_sig.pdf"); 
      PlotDebug(hNominalSelectedSignal , "hist", "hNominalSelectedSignal", plot_dump + "/debug_nom_sig_sel.pdf");
      PlotDebug(hNominalBkgd           , "hist", "hNominalBkgd"          , plot_dump + "/debug_nom_bkgd.pdf");
      PlotDebug(hNominalSelectedBkgd   , "hist", "hNominalSelectedBkgd"  , plot_dump + "/debug_nom_bkgd_sel.pdf");
      PlotDebug(hNominalSelected       , "hist", "hNominalSelected"      , plot_dump + "/debug_nom_sel.pdf");
      PlotDebug(hNominalEfficiency     , "hist", "hNominalEfficiency"    , plot_dump + "/debug_nom_eff.pdf");

      PlotDebug(sel_minus_bkgd         , "hist", "Nominal Selected - Background", plot_dump + "/debug_nom_sel_minus_bkgd.pdf");
      */
      fdebug->Close();
    //////////////////////////////////// END DEBUG //////////////////////////////////    
    }

    // plot histograms
    double maxy = 1;
    hFracUncertaintyXSec       ->SetTitle("Fractional Uncertainty on Cross-section"   );
    hFracUncertaintySelected   ->SetTitle("Fractional Stat. Uncertainty on Selection" );
    hFracUncertaintyBkgdStat   ->SetTitle("Fractional Stat. Uncertainty on Background");
    hFracUncertaintyBkgdSyst   ->SetTitle("Fractional Syst. Uncertainty on Background");
    hFracUncertaintyEfficiency ->SetTitle("Fractional Syst. Uncertainty on Efficiency");

    if(fNominalSignal->NDimensions() == 1 && method != kIntegratedInsideBounds) {

      hFracUncertaintyXSec       ->GetYaxis()->SetRangeUser(0, maxy);
      hFracUncertaintySelected   ->GetYaxis()->SetRangeUser(0, maxy);
      hFracUncertaintyBkgdStat   ->GetYaxis()->SetRangeUser(0, maxy);
      hFracUncertaintyBkgdSyst   ->GetYaxis()->SetRangeUser(0, maxy);
      hFracUncertaintyEfficiency ->GetYaxis()->SetRangeUser(0, maxy);

      hFracUncertaintyXSec       ->GetYaxis()->SetTitle("#frac{#delta#sigma}{#sigma}"                        );
      hFracUncertaintySelected   ->GetYaxis()->SetTitle("#frac{#sqrt{N_{sel}}}{N_{sel}-s#upoint N_{bkg}}"      );
      hFracUncertaintyBkgdStat   ->GetYaxis()->SetTitle("#frac{s#upoint#sqrt{N_{bkg}}}{N_{sel}-s#upoint N_{bkg}}");
      hFracUncertaintyBkgdSyst   ->GetYaxis()->SetTitle("#frac{#sqrt{s#upoint#deltaN^{syst}_{bkg}}}{N_{sel}-s#upoint N_{bkg}}");
      hFracUncertaintyEfficiency ->GetYaxis()->SetTitle("#frac{#delta#varepsilon}{#varepsilon}"              );

      hFracUncertaintyXSec       ->GetYaxis()->CenterTitle();
      hFracUncertaintySelected   ->GetYaxis()->CenterTitle();
      hFracUncertaintyBkgdStat   ->GetYaxis()->CenterTitle();
      hFracUncertaintyBkgdSyst   ->GetYaxis()->CenterTitle();
      hFracUncertaintyEfficiency ->GetYaxis()->CenterTitle();

      hFracUncertaintyXSec       ->SetLineColor(kRed);
      hFracUncertaintySelected   ->SetLineColor(kRed);
      hFracUncertaintyBkgdStat   ->SetLineColor(kRed);
      hFracUncertaintyBkgdSyst   ->SetLineColor(kRed);
      hFracUncertaintyEfficiency ->SetLineColor(kRed);

      // plot
      TCanvas * canvas = new TCanvas();
      canvas->SetLeftMargin(0.15);
      canvas->SetRightMargin(0.01);
      canvas->cd();
      std::string label = fNominalSignal->GetLabels()[0];
      hFracUncertaintyXSec       ->Draw("hist");
      canvas                     ->Print(TString::Format("%s/optimize_%s_frac_uncert_xsec.pdf",
                                                         plot_dump.c_str(),
                                                         label.c_str()));
      hFracUncertaintySelected   ->Draw("hist");
      canvas                     ->Print(TString::Format("%s/optimize_%s_frac_uncert_sel.pdf",
                                                         plot_dump.c_str(),
                                                         label.c_str()));
      hFracUncertaintyBkgdStat   ->Draw("hist");
      canvas                     ->Print(TString::Format("%s/optimize_%s_frac_uncert_bkgd_stat.pdf",
                                                         plot_dump.c_str(),
                                                         label.c_str()));
      hFracUncertaintyBkgdSyst   ->Draw("hist");
      canvas                     ->Print(TString::Format("%s/optimize_%s_frac_uncert_bkgd_syst.pdf",
                                                         plot_dump.c_str(),
                                                         label.c_str()));
      hFracUncertaintyEfficiency ->Draw("hist");
      canvas                     ->Print(TString::Format("%s/optimize_%s_frac_uncert_eff.pdf",
                                                         plot_dump.c_str(),
                                                         label.c_str()));


    }
    else if(fNominalSignal->NDimensions() == 2 || method == kIntegratedInsideBounds) {
      hFracUncertaintyXSec       ->GetZaxis()->SetRangeUser(0, maxy);
      hFracUncertaintySelected   ->GetZaxis()->SetRangeUser(0, maxy);
      hFracUncertaintyBkgdStat   ->GetZaxis()->SetRangeUser(0, maxy);
      hFracUncertaintyBkgdSyst   ->GetZaxis()->SetRangeUser(0, maxy);
      hFracUncertaintyEfficiency ->GetZaxis()->SetRangeUser(0, maxy);

      hFracUncertaintyXSec       ->GetZaxis()->SetTitle("#frac{#delta#sigma}{#sigma}"                 );
      hFracUncertaintySelected   ->GetZaxis()->SetTitle("#frac{#sqrt{N_{sel}}}{N_{sel}-N_{bkg}}"      );
      hFracUncertaintyBkgdStat   ->GetZaxis()->SetTitle("#frac{#sqrt{N_{bkg}}}{N_{sel}-N_{bkg}}"      );
      hFracUncertaintyBkgdSyst   ->GetZaxis()->SetTitle("#frac{#deltaN^{syst}_{bkg}}{N_{sel}-N_{bkg}}");
      hFracUncertaintyEfficiency ->GetZaxis()->SetTitle("#frac{#delta#varepsilon}{#varepsilon}"       );

      hFracUncertaintyXSec       ->GetYaxis()->CenterTitle();
      hFracUncertaintySelected   ->GetYaxis()->CenterTitle();
      hFracUncertaintyBkgdStat   ->GetYaxis()->CenterTitle();
      hFracUncertaintyBkgdSyst   ->GetYaxis()->CenterTitle();
      hFracUncertaintyEfficiency ->GetYaxis()->CenterTitle();

      // plot
      TCanvas * canvas = new TCanvas();
      canvas->cd();
      std::string labelx = fNominalSignal->GetLabels()[0];
      std::string labely = method == kIntegratedInsideBounds? labelx : fNominalSignal->GetLabels()[1];
      hFracUncertaintyXSec       ->Draw("colz");
      canvas                     ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_xsec.pdf",
                                                         plot_dump.c_str(),
                                                         labelx.c_str(),
                                                         labely.c_str()));
      hFracUncertaintySelected   ->Draw("colz");
      canvas                     ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_sel.pdf",
                                                         plot_dump.c_str(),
                                                         labelx.c_str(),
                                                         labely.c_str()));
      hFracUncertaintyBkgdStat   ->Draw("colz");
      canvas                     ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_bkgd_stat.pdf",
                                                         plot_dump.c_str(),
                                                         labelx.c_str(),
                                                         labely.c_str()));
      hFracUncertaintyBkgdSyst   ->Draw("colz");
      canvas                     ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_bkgd_syst.pdf",
                                                         plot_dump.c_str(),
                                                         labelx.c_str(),
                                                         labely.c_str()));
      hFracUncertaintyEfficiency ->Draw("colz");
      canvas                     ->Print(TString::Format("%s/optimize_%s_%s_frac_uncert_eff.pdf",
                                                         plot_dump.c_str(),
                                                         labelx.c_str(),
                                                         labely.c_str()));
    }
    else {
      std::cout << "Don't know how to draw a " << fNominalSignal->NDimensions() << " dimensional histogram.";
      exit(1);
    }

    // save      
    save_dir->cd();
    hFracUncertaintyXSec       ->Write("frac_uncert_xsec"     );
    hFracUncertaintySelected   ->Write("frac_uncert_sel"      );
    hFracUncertaintyBkgdStat   ->Write("frac_uncert_bkgd_stat");
    hFracUncertaintyBkgdSyst   ->Write("frac_uncert_bkgd_syst");
    hFracUncertaintyEfficiency ->Write("frac_uncert_eff"      );
  }

  //////////////////////////////////////////////////////////////
  TH1 * 
  CutOptimization::AbsUncertainty(const UpDownPair<TH1> syst, const TH1 * nominal) const
  {
    TH1 * ret = (TH1*) syst.up->Clone();
    int NX = syst.up->GetNbinsX();
    int NY = syst.up->GetNbinsY();
    int NZ = syst.up->GetNbinsZ();
    for(auto ibinx = 1; ibinx <= NX; ibinx++) {
      for(auto ibiny = 1; ibiny <= NY; ibiny++) {
        for(auto ibinz = 1; ibinz <= NZ; ibinz++) {
          double nom     = nominal->GetBinContent(ibinx, ibiny, ibinz);

          double up      = syst.up->GetBinContent(ibinx, ibiny, ibinz);
          double down    = nom;
          if(syst.down) down = syst.down->GetBinContent(ibinx, ibiny, ibinz);

          double error   = std::max(std::abs(up - nom), std::abs(down - nom));
          ret->SetBinContent(ibinx, ibiny, ibinz, error);                                   
        }
      }
    }
    
    return ret;
  }
  
  //////////////////////////////////////////////////////////////
  TH1 * 
  CutOptimization::AbsUncertaintySquared(const UpDownPair<TH1> syst, const TH1 * nominal) const
  {
    TH1 * ret = AbsUncertainty(syst, nominal);
    ret->Multiply(ret);
    return ret;
  }


  //////////////////////////////////////////////////////////////
  void 
  CutOptimization::Optimize(OptimizeMethod_t method, FOM_t fom, double data_pot,
                            TDirectory * save_dir, std::string plot_dump,
                            bool debug)
  {
    switch(fom) {
    case kdSigmaOverSigma: OptimizedSigmaOverSigma(method, data_pot, save_dir, plot_dump, debug);
      break;
    default:
      std::cout << "Your provided FOM is not yet implemented." << std::endl;
      exit(0);
    }
  }



}