PPFXHelper.h

Go to the documentation of this file.

/*
  This file contains function systematic studies
  In this file:


  std::string StringFromInt(int i, bool withE20);
  std::string StringFromDouble(double pot);

  void SqrtError(TH1* nom);
  void ZeroError(TH1* nom);

  void AddErrorInQuadrature(TH1* nom, TH1* up, bool use_max);
  void AddErrorInQuadrature(TH1* nom, TH1* up, TH1* down, bool use_max);

  TH1* GetNC(IDecomp* specs, double POT);
  TH1* GetNC(IPrediction* specs, double POT);

  TH1* GetBG(IDecomp* specs, double POT);
  TH1* GetBG(IPrediction* specs, double POT);

  void InitializeSystText(FILE* text, strings strs);

  void LoadMaps(TDirectory* dir,
                std::map<std::string,                                     IDecomp*>    * nominal,
                std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > >* shifted);
  void LoadMaps(TDirectory* dir,
                std::map<std::string,                                     PredictionNoExtrap*>    * nominal,
                std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > >* shifted);
  void LoadMaps(TDirectory* dir,
                std::map<std::string,                                     PredictionSterile*>    * nominal,
                std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >* shifted);

  void PlotSyst(TH1* h0, TH1* hp1,
                TDirectory* out, FILE* text, strings strs);
  void PlotSyst(TH1* h0, TH1* hp1, TH1* hm1, TH1* hp2, TH1* hm2,
                TDirectory* out, FILE* text, strings strs);

  void PlotSyst(TH1* nom, IDecomp* shifts,
                TDirectory* out, FILE* text, strings strs,
                double POT, bool NC = true, bool use_max);
  void PlotSyst(TH1* nom, PredictionNoExtrap* shifts,
                TDirectory* out, FILE* text, strings strs,
                double POT, bool NC = true, bool use_max);
  void PlotSyst(TH1* nom, PredictionSterile* shifts,
                TDirectory* out, FILE* text, strings strs,
                double POT, bool NC = true, bool use_max);

  void PlotSyst(TH1* nom, std::map<int, IDecomp*> shifts,
                TDirectory* out, FILE* text, strings strs,
                double POT, bool NC = true, bool use_max);
  void PlotSyst(TH1* nom, std::map<int, PredictionNoExtrap*> shifts,
                TDirectory* out, FILE* text, strings strs,
                double POT, bool NC = true, bool use_max);
  void PlotSyst(TH1* nom, std::map<int, PredictionSterile*> shifts,
                TDirectory* out, FILE* text, strings strs,
                double POT, bool NC = true, bool use_max);

  void PlotSystBand(TH1* h, TDirectory* out, FILE* text, strings strs);

  void SaveMaps(
    TDirectory* out,
    std::map<std::string, IDecomp*> decomps_nominal,
    std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted,
    std::map<std::string, PredictionNoExtrap*> predsNE_nominal,
    std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted,
    std::map<std::string, PredictionSterile*> predsSt_nominal,
    std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted
  );
*/

#pragma once

#include <map>
#include <iostream>
#include <string>
#include <vector>

#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Decomp/IDecomp.h"
#include "CAFAna/Prediction/IPrediction.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "NuXAna/Prediction/PredictionSterile.h"
#include "NuXAna/macros/NuSLoadMacroDefs.h"
#include "NuXAna/macros/NuSPlotFunctions.h"

#include "OscLib/OscCalcPMNSOpt.h"
#include "OscLib/OscCalcSterile.h"

#include "TCanvas.h"
#include "TDirectory.h"
#include "TH1.h"
#include "TKey.h"
#include "TLatex.h"
#include "TLegend.h"

namespace ana
{
  /// A helper structure to contain a group of string for plotting
  struct strings
  {
    std::string fComponent;
    std::string fDet;
    std::string fPOT;
    std::string fSample;
    std::string fSystType;
    std::string fSystS;
    std::string fSystL;
    std::string fXLabel;
  };

  //----------------------------------------------------------------------------
  /// Helper function that takes an integer and turns it into a string
  std::string StringFromInt(int i, bool withE20)
  {
    char pot[3];
    sprintf(pot, "%d", i);
    std::string potStr = std::string(pot);
    if(withE20) { potStr += " x 10^{20} POT"; }

    return potStr;
  }

  //----------------------------------------------------------------------------
  /// Helper function that takes a double POT value and turns it into a string
  /// The output has the format "# x 10^{20} POT"
  std::string StringFromDouble(double pot)
  {
    double dropE20 = pot/1e20;
    int i = (int)dropE20;
    if(dropE20 > (double)i + 0.5) { ++i; }

    return StringFromInt(i, true);
  }

  //----------------------------------------------------------------------------
  /// Set the errors of the input histogram to be their square roots
  void SqrtError(TH1* nom)
  {
    for(int i = 0, n = nom->GetNbinsX(); i <= n; ++i) {
      nom->SetBinError(i, sqrt(nom->GetBinError(i)));
    }
    return;
  }

  //----------------------------------------------------------------------------
  /// Set the errors of the input histogram to 0
  void ZeroError(TH1* nom)
  {
    for(int i = 0, n = nom->GetNbinsX(); i <= n; ++i) {
      nom->SetBinError(i, 0.);
    }
    return;
  }

  //----------------------------------------------------------------------------
  /// Combine errors in quadrature, setting error bins in nominal spectra
  /// Leaves errors squared -- square root must be applied later
  /// \param nom The nominal histogram
  /// \param up The shifted spectra
  /// \param use_max Only keep the max error, don't combine in quadrature
  void AddErrorInQuadrature(TH1* nom, TH1* up, bool use_max)
  {
    double error;
    for(int i = 0, n = nom->GetNbinsX(); i <= n; ++i) {
      error = std::abs(up->GetBinContent(i) - nom->GetBinContent(i));
      double error_sq = error*error;
      double err_prev = std::abs(nom->GetBinError(i));
      if(use_max) {
        nom->SetBinError(i, std::max(err_prev, error_sq));
      }
      else {
        nom->SetBinError(i, err_prev + error_sq);
      }
    }

    return;
  }

  //----------------------------------------------------------------------------
  /// Combine errors in quadrature, setting error bins in nominal spectra
  /// Takes the larger difference from nominal from either up or down shift
  /// Leaves errors squared -- square root must be applied later
  /// \param nom The nominal histogram
  /// \param up The +1 sigma shifted spectra
  /// \param down The -1 sigma shifted spectra
  /// \param use_max Only keep the max error, don't combine in quadrature
  void AddErrorInQuadrature(TH1* nom, TH1* up, TH1* down, bool use_max)
  {
    double error;
    for(int i = 0, n = nom->GetNbinsX(); i <= n; ++i) {
      error = std::max(std::abs(up ->GetBinContent(i) - nom ->GetBinContent(i)),
                       std::abs(nom->GetBinContent(i) - down->GetBinContent(i)));
      double error_sq = error*error;
      double err_prev = std::abs(nom->GetBinError(i));
      if(use_max) {
        nom->SetBinError(i, std::max(err_prev, error_sq));
      }
      else {
        nom->SetBinError(i, err_prev + error_sq);
      }
    }

    return;
  }

  //----------------------------------------------------------------------------
  /// Get the NC spectrum from a decomposition at a specific POT value
  /// Sets errors to 0
  TH1* GetNC(IDecomp* specs, double POT)
  {
    Spectrum s = specs->NCTotalComponent();
    TH1* ret = s.ToTH1(POT);
    ZeroError(ret);
    return ret;
  }

  //----------------------------------------------------------------------------
  /// Get the NC spectrum from a prediction object at a specific POT value
  /// Sets errors to 0
  TH1* GetNC(IPrediction* specs, double POT)
  {
    osc::OscCalcPMNSOpt* calc = (osc::OscCalcPMNSOpt*)DefaultOscCalc();
    Spectrum s = specs->PredictComponent(
      calc, Flavors::kAll, Current::kNC, Sign::kBoth
    );
    TH1* ret = s.ToTH1(POT);
    ZeroError(ret);
    return ret;
  }

  //----------------------------------------------------------------------------
  /// Get the background spectrum from a decomposition at a specific POT value
  /// Sets errors to 0
  TH1* GetBG(IDecomp* specs, double POT)
  {
    Spectrum s = specs->NueComponent();
    s += specs->AntiNueComponent();
    s += specs->NumuComponent();
    s += specs->AntiNumuComponent();
    TH1* ret = s.ToTH1(POT);
    ZeroError(ret);
    return ret;
  }

  //----------------------------------------------------------------------------
  /// Get the background spectrum from a prediction object at a specific POT value
  /// Sets errors to 0
  TH1* GetBG(IPrediction* specs, double POT)
  {
    osc::OscCalcPMNSOpt* calc = (osc::OscCalcPMNSOpt*)DefaultOscCalc();
    Spectrum s = specs->PredictComponent(
      calc, Flavors::kAll, Current::kCC, Sign::kBoth
    );
    TH1* ret = s.ToTH1(POT);
    ZeroError(ret);
    return ret;
  }

  //----------------------------------------------------------------------------
  /// Print some initial text about a systematic -- the systematic type and sample
  void InitializeSystText(FILE* text, strings strs)
  {
    fprintf(text, "%s %s\n\n", strs.fSystType.c_str(), strs.fSample.c_str());

    return;
  }

  //----------------------------------------------------------------------------
  /// Load nominal and shifted decompositions from a TDirectory
  /// The maps to contain these objects must already exist
  void LoadMaps(TDirectory* dir,
                std::map<std::string,                                     IDecomp*>    * nominal,
                std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > >* shifted)
  {
    TDirectory* tmp = gDirectory;
    TDirectory* keyDir = gDirectory;

    // Loop over all keys in this directory
    dir->cd();
    TIter nextkey(dir->GetListOfKeys());
    TKey* key;
    TKey* oldkey = 0;

    std::string sep = "__"; // This is the token separator
    std::string cafanaType = "decomp";
    std::string sampleLabel = "";
    std::string systLabel = "";
    int sigma = 0;

    while((key = (TKey*)nextkey())) {
      // Keep only the highest cycle number for each key
      if(oldkey && !strcmp(oldkey->GetName(), key->GetName())) { continue; }
      if(key->ReadObj()->IsFolder()) {
        keyDir = (TDirectory*)key->ReadObj();
      }
      else {
        continue;
      }

      std::string name(key->GetName());
      std::vector<std::string> tokens;
      std::size_t pos = 0, prev = 0;
      while(pos != std::string::npos) { // Loop over the TObject name
        pos = name.find(sep, prev); // Find the next separator

        if(pos != std::string::npos) {
          tokens.push_back(name.substr(prev, pos-prev));
          prev = pos+2; // Move past the current separator
        }
        else {
          tokens.push_back(name.substr(prev, name.length()-prev));
        }
      }

      // Make sure the object we found is the type we want
      if(cafanaType.compare(tokens[0]) != 0) { continue; }

      sampleLabel = tokens[1];
      if(tokens.size() > 2) {
        systLabel = tokens[2];
        sigma = std::stoi(tokens[3]);

        (*shifted)[sampleLabel][systLabel][sigma] = LoadFrom<IDecomp>(keyDir).release();
      }
      else {
        (*nominal)[sampleLabel] = LoadFrom<IDecomp>(keyDir).release();
      }
    }

    tmp->cd();
    return;
  }

  //----------------------------------------------------------------------------
  /// Load nominal and shifted PredictionNoExtrap objects from a TDirectory
  /// The maps to contain these objects must already exist
  void LoadMaps(TDirectory* dir,
                std::map<std::string,                                     PredictionNoExtrap*>    * nominal,
                std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > >* shifted)
  {
    TDirectory* tmp = gDirectory;
    TDirectory* keyDir = gDirectory;

    dir->cd();
    TIter nextkey(dir->GetListOfKeys());
    TKey* key;
    TKey* oldkey = 0;

    std::string sep = "__";
    std::string cafanaType = "predNE";
    std::string sampleLabel = "";
    std::string systLabel = "";
    int sigma = 0;

    while((key = (TKey*)nextkey())) {
      if(oldkey && !strcmp(oldkey->GetName(), key->GetName())) { continue; }
      if(key->ReadObj()->IsFolder()) {
        keyDir = (TDirectory*)key->ReadObj();
      }
      else {
        continue;
      }

      std::string name(key->GetName());
      std::vector<std::string> tokens;
      std::size_t pos = 0, prev = 0;
      while(pos != std::string::npos) {
        pos = name.find(sep, prev);

        if(pos != std::string::npos) {
          tokens.push_back(name.substr(prev, pos-prev));
          prev = pos+2;
        }
        else {
          tokens.push_back(name.substr(prev, name.length()-prev));
        }
      }

      if(cafanaType.compare(tokens[0]) != 0) { continue; }

      sampleLabel = tokens[1];
      if(tokens.size() > 2) {
        systLabel = tokens[2];
        sigma = std::stoi(tokens[3]);

        (*shifted)[sampleLabel][systLabel][sigma] = (PredictionNoExtrap*)PredictionExtrap::LoadFrom(keyDir).release();
      }
      else {
        (*nominal)[sampleLabel] = (PredictionNoExtrap*)PredictionExtrap::LoadFrom(keyDir).release();
      }
    }

    tmp->cd();
    return;
  }

  //----------------------------------------------------------------------------
  /// Load nominal and shifted PredictionSterile objects from a TDirectory
  /// The maps to contain these objects must already exist
  void LoadMaps(TDirectory* dir,
                std::map<std::string,                                     PredictionSterile*>    * nominal,
                std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > >* shifted)
  {
    TDirectory* tmp = gDirectory;
    TDirectory* keyDir = gDirectory;

    dir->cd();
    TIter nextkey(dir->GetListOfKeys());
    TKey* key;
    TKey* oldkey = 0;

    std::string sep = "__";
    std::string cafanaType = "predSt";
    std::string sampleLabel = "";
    std::string systLabel = "";
    int sigma = 0;

    while((key = (TKey*)nextkey())) {
      if(oldkey && !strcmp(oldkey->GetName(), key->GetName())) { continue; }
      if(key->ReadObj()->IsFolder()) {
        keyDir = (TDirectory*)key->ReadObj();
      }
      else {
        continue;
      }

      std::string name(key->GetName());
      std::vector<std::string> tokens;
      std::size_t pos = 0, prev = 0;
      while(pos != std::string::npos) {
        pos = name.find(sep, prev);

        if(pos != std::string::npos) {
          tokens.push_back(name.substr(prev, pos-prev));
          prev = pos+2;
        }
        else {
          tokens.push_back(name.substr(prev, name.length()-prev));
        }
      }

      if(cafanaType.compare(tokens[0]) != 0) { continue; }

      sampleLabel = tokens[1];
      if(tokens.size() > 2) {
        systLabel = tokens[2];
        sigma = std::stoi(tokens[3]);

        (*shifted)[sampleLabel][systLabel][sigma] = PredictionSterile::LoadFrom(keyDir).release();
      }
      else {
        (*nominal)[sampleLabel] = PredictionSterile::LoadFrom(keyDir).release();
      }
    }

    tmp->cd();
    return;
  }

  //----------------------------------------------------------------------------
  /// \param F
  void PlotMultiSyst(TH1* h, std::vector<TH1*> hp, TDirectory* out, FILE* text, strings strs)
  {
    // Construct histogram names and titles
    std::string xlabel = strs.fXLabel;
    std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
    std::string title = strs.fDet + " " + strs.fComponent +
                        " Error for " + strs.fSystL + " Systematic";
    std::string fullTitle = ";" + xlabel + ";" + ylabel;
//    std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
    std::string name = "c" + strs.fComponent + strs.fDet +
                       strs.fSystS + strs.fSample;
    std::string ylabelRat = "Shifted / Nominal";
    std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;
 /* TH1* hRat;
  std::vector<TH1*>hp1V;
  std::vector<TH1*>hp1RatV;
  hp1V.clear();
  hp1RatV.clear();
 */ // for(const auto hp1 :hp){
    TH1* hCVF = hp.at(0);
    TH1* hCV = (TH1*)h->Clone();
    map< int,float > cv_events_in_bins;
    map<int, vector<float> > univ_events_in_bins;
    unsigned int count = 0;

    for(auto& it: hp){
            if(count!= hp.size()){
                    TH1* curhist = it;
                    int nbinsx = curhist->GetNbinsX();
                    for(int i = 1; i <= nbinsx; i++){
                            univ_events_in_bins[i].push_back(curhist->GetBinContent(i));
                    }
            }
            else if(count== hp.size()){
                    int nbinsx = hCV->GetNbinsX();
                    for(int i = 1; i <= nbinsx; i++){
                            cv_events_in_bins[i] = hCV->GetBinContent(i);
                    }
            }
            count++;
    }

    TH1* hUp = (TH1*)hCVF->Clone("hUp");
    TH1* hLow = (TH1*)hCVF->Clone("hLow");
    int binIt = 1;
    for(auto& it: univ_events_in_bins){
            double cv_val = hCVF->GetBinContent( it.first);
            double err = 0.0;
            for(unsigned int iuniv=0;iuniv<hp.size();iuniv++){
                    err += pow(it.second[iuniv] - cv_val,2);
            }
            err /= double(hp.size());
            err = sqrt(err);
            hUp->SetBinContent(binIt, cv_val+err);
            hLow->SetBinContent(binIt,cv_val-err);
            binIt++;
    }

//    TH1* hp1=hp.at(0);
    TH1* hRat = (TH1*)h->Clone();
  //  TH1* hp1Rat = (TH1*)hp1->Clone();
    TH1* hUpRat=(TH1*)hUp->Clone();
     TH1* hLowRat=(TH1*)hLow->Clone();
    int nBins = h->GetNbinsX();
    double percent_diff = 100. * std::abs(hUp->Integral(1, nBins) - h->Integral(1, nBins))/
                          h->Integral(1, nBins);
     double percent_diff1 = 100. * std::abs(hLow->Integral(1, nBins) - h->Integral(1, nBins))/
                          h->Integral(1, nBins);

    fprintf(text, "Up %s %s %.2s: +/-%6.4f%%\n",
            strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
            percent_diff);
     fprintf(text, "Low %s %s %.2s: +/-%6.4f%%\n",
            strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
            percent_diff1);


    for(int i = 0; i <= nBins+1; ++i) {
      hRat->SetBinContent(i, 1.);
    }
    hUpRat->Divide(h);
    hLowRat->Divide(h);

    double maxval = h->GetBinContent(h->GetMaximumBin());
    maxval = std::max(maxval, hUp->GetBinContent(hUp->GetMaximumBin()));

    double minval = h->GetBinContent(h->GetMinimumBin());
    minval = std::min(minval, hLow->GetBinContent(hLow->GetMinimumBin()));

    if(maxval < minval) { std::swap(maxval, minval); }

    double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
    maxvalRat = std::max(maxvalRat, hUpRat->GetBinContent(hUpRat->GetMaximumBin()));

    double minvalRat = hLowRat->GetBinContent(hLowRat->GetMinimumBin());
    for(int i = 1; i <= nBins; ++i) {
      if(hUpRat->GetBinContent(i) > 0.) {
        minvalRat = std::min(minvalRat, hLowRat->GetBinContent(i));
      }
    }

    if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
      minvalRat = 1 - std::abs(maxvalRat - 1.);
      maxvalRat = 1 + std::abs(maxvalRat - 1.);
    }
    else {
      minvalRat = 1 - std::abs(1. - minvalRat);
      maxvalRat = 1 + std::abs(1. - minvalRat);
    }

   CenterTitles(h);
    h->SetLineColor(kBlack);
    h->SetMaximum(1.1*maxval);
    h->SetMinimum(0);
    h->SetTitle(fullTitle.c_str());

    CenterTitles(hLow);
    hLow->SetLineColor(kRed);
    hLow->SetLineStyle(2);
    hLow->SetMaximum(1.1*maxval);
    hLow->SetMinimum(0);
    hLow->SetTitle(fullTitle.c_str());

    CenterTitles(hUp);
    hUp->SetLineColor(kGreen+2);
    hUp->SetLineStyle(2);
    hUp->SetMaximum(1.1*maxval);
    hUp->SetMinimum(0);
    hUp->SetTitle(fullTitle.c_str());


    CenterTitles(hRat);
    hRat->SetLineColor(kBlack);
    hRat->SetMaximum(1.1*maxvalRat);
    hRat->SetMinimum(0.9*minvalRat);
    hRat->SetTitle(fullTitleRat.c_str());

    CenterTitles(hUpRat);
    hUpRat->SetLineColor(kGreen+2);
    hUpRat->SetLineStyle(2);
    hUpRat->SetMaximum(1.1*maxvalRat);
    hUpRat->SetMinimum(0.9*minvalRat);

      CenterTitles(hLowRat);
    hLowRat->SetLineColor(kRed);
    hLowRat->SetLineStyle(2);
    hLowRat->SetMaximum(1.1*maxvalRat);
    hLowRat->SetMinimum(0.9*minvalRat);
    //hp1Rat->SetTitle(fullTitleRat.c_str());
  //  hp1RatV.push_back(hp1Rat);
  //  hp1V.push_back(hp1);
//}
    double xL = 0.6, xR = 0.85;
    double yB = 0.6, yT = 0.85;
    TLegend* leg = new TLegend(xL, yB, xR, yT);
    leg->AddEntry(h,   "Nominal",      "l");
    leg->AddEntry(hUp, "Shifted Up 34%", "l");
    leg->AddEntry(hLow, "Shifted down 34%", "l");

 // for(const auto hp1 :hp){
    TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
    TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
    TPad* pRatio = new TPad("pRatio", "", 0., 0.,    1., 0.375);
    pSpecs->Draw();
    pRatio->Draw();

    pSpecs->cd();
    h->Draw("hist");
    hUp->Draw("hist same");
    hLow->Draw("hist same");
     c->Update();
    leg->Draw();
    Simulation();

    pRatio->cd();
    hRat->Draw("hist");
//   for(auto hp1Rat:hp1RatV){
    hUpRat->Draw("hist same");
    hLowRat->Draw("hist same");
//   }
    c->Update();
//}
    out->WriteTObject(c);
 // } 
   return;
  }

  //----------------------------------------------------------------------------
  /// Plot a nominal spectrum and a shifted spectra
  /// This is for a systematic with multiple sigma levels
  /// \param h The nominal spectrum
  /// \param hp1 The shifted spectrum, +1 sigma
  /// \param hm1 The shifted spectrum, -1 sigma
  /// \param hp2 The shifted spectrum, +2 sigma
  /// \param hm2 The shifted spectrum, -2 sigma
  /// \param out The directory/file for root output
  /// \param text The file for text output
  /// \param strs Contains strings for pretty plots
  void PlotSyst(TH1* h, TH1* hp1, TH1* hm1, TH1* hp2, TH1* hm2,
                TDirectory* out, FILE* text, strings strs)
  {
    std::string xlabel = strs.fXLabel;
    std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
    std::string title = strs.fDet + " " + strs.fComponent +
                        " Error for " + strs.fSystL + " Systematic";
    std::string fullTitle = ";" + xlabel + ";" + ylabel;
//    std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
    std::string name = "c" + strs.fComponent + strs.fDet +
                       strs.fSystS + strs.fSample;
    std::string ylabelRat = "Shifted / Nominal";
    std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;

    TH1* hRat = (TH1*)h->Clone();
    TH1* hp1Rat = (TH1*)hp1->Clone();
    TH1* hm1Rat = (TH1*)hm1->Clone();
    TH1* hp2Rat = (TH1*)hp2->Clone();
    TH1* hm2Rat = (TH1*)hm2->Clone();

    int nBins = h->GetNbinsX();
    double nom_evts = h->Integral(1, nBins);
    double evts_diff = 0.;
    double percent_diff   = 100.;
    double percent_diff_p = 100.;
    double percent_diff_m = 100.;
    for(int i = 0; i <= nBins+1; ++i) {
      hRat->SetBinContent(i, 1.);

      if(i != 0 && i != nBins+1) {
        evts_diff += std::max(std::abs(hp1->GetBinContent(i) - h->GetBinContent(i)),
                              std::abs(hm1->GetBinContent(i) - h->GetBinContent(i)));
      }
    }
    percent_diff   *= (evts_diff/nom_evts);
    percent_diff_p *= ((hp1->Integral(1, nBins) - nom_evts)/nom_evts);
    percent_diff_m *= ((nom_evts - hm1->Integral(1, nBins))/nom_evts);

    fprintf(text, "%s %.2s %.2s: +/-%6.4f%%, +%6.4f%%, -%6.4f%%\n",
            strs.fSystS.c_str(), strs.fComponent.c_str(), strs.fDet.c_str(),
            percent_diff, percent_diff_p, percent_diff_m);

    hp1Rat->Divide(h);
    hm1Rat->Divide(h);
    hp2Rat->Divide(h);
    hm2Rat->Divide(h);

    double maxval = h->GetBinContent(h->GetMaximumBin());
    maxval = std::max(maxval, hp1->GetBinContent(hp1->GetMaximumBin()));
    maxval = std::max(maxval, hm1->GetBinContent(hm1->GetMaximumBin()));
    maxval = std::max(maxval, hp2->GetBinContent(hp2->GetMaximumBin()));
    maxval = std::max(maxval, hm2->GetBinContent(hm2->GetMaximumBin()));

    double minval = h->GetBinContent(h->GetMinimumBin());
    minval = std::min(minval, hp1->GetBinContent(hp1->GetMinimumBin()));
    minval = std::min(minval, hm1->GetBinContent(hm1->GetMinimumBin()));
    minval = std::min(minval, hp2->GetBinContent(hp2->GetMinimumBin()));
    minval = std::min(minval, hm2->GetBinContent(hm2->GetMinimumBin()));

    if(maxval < minval) { std::swap(maxval, minval); }

    double maxvalRat = hRat->GetBinContent(hRat->GetMaximumBin());
    maxvalRat = std::max(maxvalRat, hp1Rat->GetBinContent(hp1Rat->GetMaximumBin()));
    maxvalRat = std::max(maxvalRat, hm1Rat->GetBinContent(hm1Rat->GetMaximumBin()));
    maxvalRat = std::max(maxvalRat, hp2Rat->GetBinContent(hp2Rat->GetMaximumBin()));
    maxvalRat = std::max(maxvalRat, hm2Rat->GetBinContent(hm2Rat->GetMaximumBin()));

    double minvalRat = hRat->GetBinContent(hRat->GetMinimumBin());
    for(int i = 1; i <= nBins; ++i) {
      if(hp1Rat->GetBinContent(i) > 0.) {
        minvalRat = std::min(minvalRat, hp1Rat->GetBinContent(i));
      }
      if(hm1Rat->GetBinContent(i) > 0.) {
        minvalRat = std::min(minvalRat, hm1Rat->GetBinContent(i));
      }
      if(hp2Rat->GetBinContent(i) > 0.) {
        minvalRat = std::min(minvalRat, hp2Rat->GetBinContent(i));
      }
      if(hm2Rat->GetBinContent(i) > 0.) {
        minvalRat = std::min(minvalRat, hm2Rat->GetBinContent(i));
      }
    }

    if(std::abs(maxvalRat - 1.) > std::abs(1. - minvalRat)) {
      minvalRat = 1 - std::abs(maxvalRat - 1.);
      maxvalRat = 1 + std::abs(maxvalRat - 1.);
    }
    else {
      minvalRat = 1 - std::abs(1. - minvalRat);
      maxvalRat = 1 + std::abs(1. - minvalRat);
    }

    CenterTitles(h);
    h->SetLineColor(kBlack);
    h->SetMaximum(1.1*maxval);
    h->SetMinimum(0);
    h->SetTitle(fullTitle.c_str());

    CenterTitles(hp1);
    hp1->SetLineColor(kRed);
    hp1->SetLineStyle(2);
    hp1->SetMaximum(1.1*maxval);
    hp1->SetMinimum(0);
    hp1->SetTitle(fullTitle.c_str());

    CenterTitles(hm1);
    hm1->SetLineColor(kRed);
    hm1->SetLineStyle(2);
    hm1->SetMaximum(1.1*maxval);
    hm1->SetMinimum(0);
    hm1->SetTitle(fullTitle.c_str());

    CenterTitles(hp2);
    hp2->SetLineColor(kBlue);
    hp2->SetLineStyle(2);
    hp2->SetMaximum(1.1*maxval);
    hp2->SetMinimum(0);
    hp2->SetTitle(fullTitle.c_str());

    CenterTitles(hm2);
    hm2->SetLineColor(kBlue);
    hm2->SetLineStyle(2);
    hm2->SetMaximum(1.1*maxval);
    hm2->SetMinimum(0);
    hm2->SetTitle(fullTitle.c_str());

    CenterTitles(hRat);
    hRat->SetLineColor(kBlack);
    hRat->SetMaximum(1.1*maxvalRat);
    hRat->SetMinimum(0.9*minvalRat);
    hRat->SetTitle(fullTitleRat.c_str());

    CenterTitles(hp1Rat);
    hp1Rat->SetLineColor(kRed);
    hp1Rat->SetLineStyle(2);
    hp1Rat->SetMaximum(1.1*maxvalRat);
    hp1Rat->SetMinimum(0.9*minvalRat);
    hp1Rat->SetTitle(fullTitleRat.c_str());

    CenterTitles(hm1Rat);
    hm1Rat->SetLineColor(kRed);
    hm1Rat->SetLineStyle(2);
    hm1Rat->SetMaximum(1.1*maxvalRat);
    hm1Rat->SetMinimum(0.9*minvalRat);
    hm1Rat->SetTitle(fullTitleRat.c_str());

    CenterTitles(hp2Rat);
    hp2Rat->SetLineColor(kBlue);
    hp2Rat->SetLineStyle(2);
    hp2Rat->SetMaximum(1.1*maxvalRat);
    hp2Rat->SetMinimum(0.9*minvalRat);
    hp2Rat->SetTitle(fullTitleRat.c_str());

    CenterTitles(hm2Rat);
    hm2Rat->SetLineColor(kBlue);
    hm2Rat->SetLineStyle(2);
    hm2Rat->SetMaximum(1.1*maxvalRat);
    hm2Rat->SetMinimum(0.9*minvalRat);
    hm2Rat->SetTitle(fullTitleRat.c_str());

    double xL = 0.6, xR = 0.85;
    double yB = 0.6, yT = 0.85;
    TLegend* leg = new TLegend(xL, yB, xR, yT);
    leg->AddEntry(h,   "Nominal",      "l");
    leg->AddEntry(hp1, "+/- 1 #sigma", "l");
    leg->AddEntry(hp2, "+/- 2 #sigma", "l");

    TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
    TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
    TPad* pRatio = new TPad("pRatio", "", 0., 0.,    1., 0.375);
    pSpecs->Draw();
    pRatio->Draw();

    pSpecs->cd();
    h->Draw("hist");
    hp2->Draw("hist same");
    hm2->Draw("hist same");
    hp1->Draw("hist same");
    hm1->Draw("hist same");
    leg->Draw();
    Simulation();

    pRatio->cd();
    hRat->Draw("hist");
    hp2Rat->Draw("hist same");
    hm2Rat->Draw("hist same");
    hp1Rat->Draw("hist same");
    hm1Rat->Draw("hist same");

    c->Update();
    out->WriteTObject(c);
    return;
  }

  //----------------------------------------------------------------------------
  /// Plot ND syst shifted and nominal spectra
  /// This is for plotting an on/off style shift
  void PlotMultiSyst(TH1* nom, std::vector<IDecomp*> shifts,
                TDirectory* out, FILE* text, strings strs,
                double POT, bool NC, bool use_max)
  {

    std::vector<TH1*> hp;
   for( auto shift: shifts){
    TH1* hp1 =(NC ? GetNC(shift, POT) : GetBG(shift, POT));
    hp.push_back(hp1);
}
  //  AddErrorInQuadrature(nom, hp, use_max);
    PlotMultiSyst(nom, hp, out, text, strs);
    return;
  }

  //----------------------------------------------------------------------------
  /// Plot FD syst shifted and nominal spectra
  /// Should be able to change the input to IPrediction
  /// This is for plotting an on/off style shift
  void PlotMultiSyst(TH1* nom, std::vector<PredictionNoExtrap*> shifts,
                TDirectory* out, FILE* text, strings strs,
                double POT, bool NC, bool use_max)
  {
    std::vector<TH1*> hp;
    
for(auto shift: shifts){
    TH1* hp1=(NC ? GetNC(shift, POT) : GetBG(shift, POT));
     
    hp.push_back(hp1);
}
//    AddErrorInQuadrature(nom, hp, use_max);
    PlotMultiSyst(nom, hp, out, text, strs);
    return;
  }

  //----------------------------------------------------------------------------
  /// Plot FD syst shifted and nominal spectra
  /// Should be able to change the input to IPrediction and make this redundant
  /// This is for plotting an on/off style shift
  void PlotMultiSyst(TH1* nom, std::vector<PredictionSterile*> shifts,
                TDirectory* out, FILE* text, strings strs,
                double POT, bool NC, bool use_max)
  {
   std::vector<TH1*> hp;
     for(auto shift: shifts){
   TH1* hp1 =(NC ? GetNC(shift, POT) : GetBG(shift, POT));
    hp.push_back(hp1);
}
   // AddErrorInQuadrature(nom, hp, use_max);
    PlotMultiSyst(nom, hp, out, text, strs);
    return;
  }

  //----------------------------------------------------------------------------
  /// Plot ND syst shifted and nominal spectra
  /// This is for plotting a syst with multiple sigma levels
  void PlotSyst(TH1* nom, std::map<int, IDecomp*> shifts,
                TDirectory* out, FILE* text, strings strs,
                double POT, bool NC, bool use_max)
  {
    TH1* hp1 = (NC ? GetNC(shifts[+1], POT) : GetBG(shifts[+1], POT));
    TH1* hm1 = (NC ? GetNC(shifts[-1], POT) : GetBG(shifts[-1], POT));
    TH1* hp2 = (NC ? GetNC(shifts[+2], POT) : GetBG(shifts[+2], POT));
    TH1* hm2 = (NC ? GetNC(shifts[-2], POT) : GetBG(shifts[-2], POT));

    AddErrorInQuadrature(nom, hp1, hm1, use_max);
    PlotSyst(nom, hp1, hm1, hp2, hm2, out, text, strs);
    return;
  }

  //----------------------------------------------------------------------------
  /// Plot FD syst shifted and nominal spectra
  /// Should be able to change the input to IPrediction
  /// This is for plotting a syst with multiple sigma levels
  void PlotSyst(TH1* nom, std::map<int, PredictionNoExtrap*> shifts,
                TDirectory* out, FILE* text, strings strs,
                double POT, bool NC, bool use_max)
  {
    TH1* hp1 = (NC ? GetNC(shifts[+1], POT) : GetBG(shifts[+1], POT));
    TH1* hm1 = (NC ? GetNC(shifts[-1], POT) : GetBG(shifts[-1], POT));
    TH1* hp2 = (NC ? GetNC(shifts[+2], POT) : GetBG(shifts[+2], POT));
    TH1* hm2 = (NC ? GetNC(shifts[-2], POT) : GetBG(shifts[-2], POT));

    AddErrorInQuadrature(nom, hp1, hm1, use_max);
    PlotSyst(nom, hp1, hm1, hp2, hm2, out, text, strs);
    return;
  }

  //----------------------------------------------------------------------------
  /// Plot FD syst shifted and nominal spectra
  /// Should be able to change the input to IPrediction and make this redundant
  /// This is for plotting a syst with multiple sigma levels
  void PlotSyst(TH1* nom, std::map<int, PredictionSterile*> shifts,
                TDirectory* out, FILE* text, strings strs,
                double POT, bool NC, bool use_max)
  {
    TH1* hp1 = (NC ? GetNC(shifts[+1], POT) : GetBG(shifts[+1], POT));
    TH1* hm1 = (NC ? GetNC(shifts[-1], POT) : GetBG(shifts[-1], POT));
    TH1* hp2 = (NC ? GetNC(shifts[+2], POT) : GetBG(shifts[+2], POT));
    TH1* hm2 = (NC ? GetNC(shifts[-2], POT) : GetBG(shifts[-2], POT));

    AddErrorInQuadrature(nom, hp1, hm1, use_max);
    PlotSyst(nom, hp1, hm1, hp2, hm2, out, text, strs);
    return;
  }

  //----------------------------------------------------------------------------
  /// Plot the systematic error band around a nominal spectrum
  /// The errors of the input histogram should be squared
  void PlotSystBand(TH1* h, TDirectory* out, FILE* text, strings strs)
  {
    // Fix the errors by taking the square root!
    SqrtError(h);

    std::string xlabel = strs.fXLabel;
    std::string ylabel = "Events / " + strs.fPOT + " / 0.25 GeV";
    std::string title = strs.fDet + " " + strs.fComponent +
                        " Error Band for All " + strs.fSystType + " Systematics";
    std::string fullTitle = ";" + xlabel + ";" + ylabel;
//    std::string fullTitle = title + ";" + xlabel + ";" + ylabel;
    std::string name = "c" + strs.fComponent + strs.fDet +
                       strs.fSystType + "Systs" + strs.fSample;
    std::string ylabelRat = "Shifted / Nominal";
    std::string fullTitleRat = ";" + xlabel + ";" + ylabelRat;

    int nBins = h->GetNbinsX();
    double maxval = 0.;
    double nom_evts = h->Integral(1, nBins);
    double evts_diff = 0.;
    double percent_diff = 100.;
    for(int i = 1; i <= nBins; ++i) {
      maxval = std::max(maxval, h->GetBinContent(i) + h->GetBinError(i));
      evts_diff += std::abs(h->GetBinError(i));
    }
    percent_diff *= (evts_diff/nom_evts);

    fprintf(text, "\nAll systs %.2s %.2s: +/-%6.4f%%\n",
            strs.fComponent.c_str(), strs.fDet.c_str(),
            percent_diff);

    CenterTitles(h);
    h->SetLineColor(kTotalMCColor);
    h->SetMaximum(1.1*maxval);
    h->SetMinimum(0);
    h->SetTitle(fullTitle.c_str());

    TH1* hErr = (TH1*)h->Clone();
    CenterTitles(hErr);
    hErr->SetFillColor(kTotalMCColor-9);
    hErr->SetMaximum(1.1*maxval);
    hErr->SetMinimum(0);

    double maxvalRat = 1.;
    TH1* hRat = (TH1*)h->Clone();
    TH1* hUpp = (TH1*)h->Clone();
    TH1* hDwn = (TH1*)h->Clone();
    for(int i = 1, n = hRat->GetNbinsX(); i <= n; ++i) {
      double binVal = hRat->GetBinContent(i);
      double errVal = hRat->GetBinError(i);
      if(binVal <= 0.) {
        binVal = 1.;
        errVal = 0.;
      }

      hRat->SetBinContent(i, 1);
      hUpp->SetBinContent(i, (binVal + errVal)/binVal);
      hDwn->SetBinContent(i, (binVal - errVal)/binVal);

      double diff = (binVal + errVal)/binVal;
      if(diff > maxvalRat) {
        maxvalRat = diff;
      }
    }
    double minvalRat = 1. - (maxvalRat-1.);

    CenterTitles(hRat);
    hRat->SetLineColor(kBlack);
    hRat->SetMaximum(1.05*maxvalRat);
    hRat->SetMinimum(0.95*minvalRat);
    hRat->SetTitle(fullTitleRat.c_str());

    CenterTitles(hUpp);
    hUpp->SetLineColor(kRed);
    hUpp->SetLineStyle(2);
    hUpp->SetMaximum(1.05*maxvalRat);
    hUpp->SetMinimum(0.95*minvalRat);
    hUpp->SetTitle(fullTitleRat.c_str());

    CenterTitles(hDwn);
    hDwn->SetLineColor(kRed);
    hDwn->SetLineStyle(2);
    hDwn->SetMaximum(1.05*maxvalRat);
    hDwn->SetMinimum(0.95*minvalRat);
    hDwn->SetTitle(fullTitleRat.c_str());

    TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 800);
    TPad* pSpecs = new TPad("pSpecs", "", 0., 0.375, 1., 1.);
    TPad* pRatio = new TPad("pRatio", "", 0., 0.,    1., 0.375);
    pSpecs->Draw();
    pRatio->Draw();

    pSpecs->cd();
    hErr->Draw("e2");
    h->Draw("hist same");
    Simulation();

    pRatio->cd();
    hRat->Draw("hist");
    hUpp->Draw("hist same");
    hDwn->Draw("hist same");

    c->Update();
    out->WriteTObject(c);
    return;
  }

  //----------------------------------------------------------------------------
  /// Save all of the objects in the input maps to the out directory/file
  void SaveMaps(
    TDirectory* out,
    std::map<std::string, IDecomp*> decomps_nominal,
    std::map<std::string, std::map<std::string, std::map<int, IDecomp*> > > decomps_shifted,
    std::map<std::string, PredictionNoExtrap*> predsNE_nominal,
    std::map<std::string, std::map<std::string, std::map<int, PredictionNoExtrap*> > > predsNE_shifted,
    std::map<std::string, PredictionSterile*> predsSt_nominal,
    std::map<std::string, std::map<std::string, std::map<int, PredictionSterile*> > > predsSt_shifted
  )
  {
    std::map<std::string, TDirectory*> savedirs;
    std::string dir;
    std::string sep = "__";

    for(const auto& decomp : decomps_nominal) {
      dir = "decomp" + sep + decomp.first; // sample name
      savedirs[dir] = out->mkdir(dir.c_str());
      decomp.second->SaveTo(savedirs[dir]);
    }

    for(const auto& sample : decomps_shifted) {
      for(const auto& shiftlabel : sample.second) {
        for(const auto& decomp : shiftlabel.second) {
          dir = "decomp" + sep + sample.first +
                sep + shiftlabel.first + sep +
                StringFromInt(decomp.first, false);
          savedirs[dir] = out->mkdir(dir.c_str());
          decomp.second->SaveTo(savedirs[dir]);
        }
      }
    }

    for(const auto& pred : predsNE_nominal) {
      dir = "predNE" + sep + pred.first;
      savedirs[dir] = out->mkdir(dir.c_str());
      pred.second->SaveTo(savedirs[dir]);
    }

    for(const auto& sample : predsNE_shifted) {
      for(const auto& shiftlabel : sample.second) {
        for(const auto& pred : shiftlabel.second) {
          dir = "predNE" + sep + sample.first +
                sep + shiftlabel.first + sep +
                StringFromInt(pred.first, false);
          savedirs[dir] = out->mkdir(dir.c_str());
          pred.second->SaveTo(savedirs[dir]);
        }
      }
    }

    for(const auto& pred : predsSt_nominal) {
      dir = "predSt" + sep + pred.first;
      savedirs[dir] = out->mkdir(dir.c_str());
      pred.second->SaveTo(savedirs[dir]);
    }

    for(const auto& sample : predsSt_shifted) {
      for(const auto& shiftlabel : sample.second) {
        for(const auto& pred : shiftlabel.second) {
          dir = "predSt" + sep + sample.first +
                sep + shiftlabel.first + sep +
                StringFromInt(pred.first, false);
          savedirs[dir] = out->mkdir(dir.c_str());
          pred.second->SaveTo(savedirs[dir]);
        }
      }
    }

    return;
  }
}