systematics_summary_from_pred_interp.C

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////////////
//
// Macro for making the nue systematics bar plots
//
////////////////////////////////////////////////////////////////////////////////

//#include "CAFAna/Analysis/Calcs.h"
//#include "CAFAna/Analysis/Exposures.h"
//#include "CAFAna/Analysis/Plots.h"

//#include "OscLib/IOscCalc.h"
//#include "OscLib/OscCalcDumb.h"

//#include "3FlavorAna/Ana2020/joint_fit_style_tools.h"
//#include "3FlavorAna/Ana2020/joint_fit_2020_loader_tools.h"
//#include "CAFAna/Analysis/3FlavorAna2020Systs.h"

//#include "3FlavorAna/Systs/DummySystStorage.h"
#include "systematics_summary_utils.h"

#include "TCanvas.h"
#include "TGraphAsymmErrors.h"

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

using namespace ana;

void systematics_summary_from_pred_interp(std::string beam = "rhc", bool IsNuMu = false, bool sum = true, bool extrapComp = false, bool quantsComp = false) 
{
  std::vector<ColumnDef> cols;
  double pot = 0;
  double livetime = 0;
  const bool beamLab = beam == "fhc" ? true : false;
  std::string sampLab = "";
  std::vector<std::vector<const ISyst *>> systs;
  std::vector<std::string> extrapStr;
  std::vector<std::string> extrapLab;
  
  // Do you want to make extrap/noextrap comparisons?
  if(extrapComp) 
  {
    extrapStr.push_back("noextrap_noPt");
    extrapLab.push_back("Not Extrapolated");
  }
  std::string sNuMu = IsNuMu? "NuMu":"NuE";
  
  // Definitions of table columns follow, use the convention or you'll get wrong numbers:
  // 0 = Total Pred
  // 1 = Signal
  // 2 = Total Bkg
  // 3 = Beam Bkg
  // last = Cosmics
  if( beam == "fhc" ) {
    pot   = kAna2020FHCPOT;
    livetime = kAna2020FHCLivetime;
    if (IsNuMu) {
      cols = {
          {"Total Pred",         Flavors::kAll,      Current::kBoth, Sign::kBoth},
          {"$\\nu_\\mu$ Signal",  Flavors::kNuMuToNuMu,  Current::kCC,   Sign::kBoth},
          {"Total Bkg",         Flavors::kAll,      Current::kBoth, Sign::kBoth},
              {"Beam Bkg",          Flavors::kAll,      Current::kBoth, Sign::kBoth},
          {"$\\nu_\\mu$ App.",  Flavors::kNuEToNuMu,Current::kCC,   Sign::kBoth},
              {"NC",                Flavors::kAll,      Current::kNC,   Sign::kBoth},
              {"$\\nu_e$ CC",       Flavors::kAllNuE,   Current::kCC,   Sign::kBoth},
              {"$\\nu_\\tau$ CC",   Flavors::kAllNuTau, Current::kCC,   Sign::kBoth},
          {"Cosmics",           Flavors::kAll,      Current::kBoth, Sign::kBoth}
           
      };
      if(extrapComp) systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNumuAna, true, BeamType2020::kFHC, false, false ) );
      systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNumuAna, true, BeamType2020::kFHC, false, true  ) );
      extrapStr.push_back("extrap");
      extrapLab.push_back("Extrapolated");
    } else {
      cols = {
          {"Total Pred",         Flavors::kAll,      Current::kBoth, Sign::kBoth},
          {"$\\nu_e$ Signal",     Flavors::kNuMuToNuE, Current::kCC,   Sign::kNu},
          {"Total Bkg",         Flavors::kAll,      Current::kBoth, Sign::kBoth},
              {"Beam Bkg",            Flavors::kAll,       Current::kBoth, Sign::kBoth},
              {"$\\bar\\nu_e$ WS",    Flavors::kNuMuToNuE, Current::kCC,   Sign::kAntiNu},
              {"beam $\\nu_e$ CC",    Flavors::kNuEToNuE,  Current::kCC,   Sign::kBoth},
              {"NC",                  Flavors::kAll,       Current::kNC,   Sign::kBoth},
              {"$\\nu_\\mu$ CC",      Flavors::kAllNuMu,   Current::kCC,   Sign::kBoth},
              {"$\\nu_\\tau$ CC",     Flavors::kAllNuTau,  Current::kCC,   Sign::kBoth},
          {"Cosmics",             Flavors::kAll,       Current::kBoth, Sign::kBoth}
      };
      if(extrapComp)
      {
        systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNueAna, true, BeamType2020::kFHC, false, false ) );
        systs.back().push_back(&kRockScaleSyst);
      }
      systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNueAna, true, BeamType2020::kFHC, false, true ) );
      systs.back().push_back(&kRockScaleSyst);
      extrapStr.push_back("combo");
      extrapLab.push_back("Extrapolated");
    }
  } else { // if RHC
    pot   = kAna2020RHCPOT;
    livetime = kAna2020RHCLivetime;
    if (IsNuMu) {
      cols = {
          {"Total Pred",         Flavors::kAll,      Current::kBoth, Sign::kBoth},
          {"$\\bar\\nu_\\mu$ Signal",  Flavors::kNuMuToNuMu,  Current::kCC,   Sign::kBoth},
          {"Total Bkg",         Flavors::kAll,      Current::kBoth, Sign::kBoth},
              {"Beam Bkg",          Flavors::kAll,      Current::kBoth, Sign::kBoth},
          {"$\\nu_\\mu$ App.",  Flavors::kNuEToNuMu,Current::kCC,   Sign::kBoth},
              {"NC",                Flavors::kAll,      Current::kNC,   Sign::kBoth},
              {"$\\nu_e$ CC",       Flavors::kAllNuE,   Current::kCC,   Sign::kBoth},
              {"$\\nu_\\tau$ CC",   Flavors::kAllNuTau, Current::kCC,   Sign::kBoth},
          {"Cosmics",           Flavors::kAll,      Current::kBoth, Sign::kBoth}
      };
      if(extrapComp) systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNumuAna, true, BeamType2020::kRHC, false, false ) );
      systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNumuAna, true, BeamType2020::kRHC, false, true ) );
      extrapStr.push_back("extrap");
      extrapLab.push_back("Extrapolated");
    } else {
      cols = {
          {"Total Pred",         Flavors::kAll,      Current::kBoth, Sign::kBoth},
          {"$\\bar\\nu_e$ Signal",Flavors::kNuMuToNuE, Current::kCC,   Sign::kAntiNu},
          {"Total Bkg",           Flavors::kAll,      Current::kBoth, Sign::kBoth},
              {"Beam Bkg",            Flavors::kAll,       Current::kBoth, Sign::kBoth},
              {"$\\nu_e$ WS",         Flavors::kNuMuToNuE, Current::kCC,   Sign::kNu},
              {"beam $\\nu_e$ CC",    Flavors::kNuEToNuE,  Current::kCC,   Sign::kBoth},
              {"NC",                  Flavors::kAll,       Current::kNC,   Sign::kBoth},
              {"$\\nu_\\mu$ CC",      Flavors::kAllNuMu,   Current::kCC,   Sign::kBoth},
              {"$\\nu_\\tau$ CC",     Flavors::kAllNuTau,  Current::kCC,   Sign::kBoth},
          {"Cosmics",             Flavors::kAll,       Current::kBoth, Sign::kBoth}
      };
      if(extrapComp)
      {
        systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNueAna, true, BeamType2020::kRHC, false, false ) );
        systs.back().push_back(&kRockScaleSyst);
      }
      systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNueAna, true, BeamType2020::kRHC, false, true ) );
      systs.back().push_back(&kRockScaleSyst);
      extrapStr.push_back("prop");
      extrapLab.push_back("Extrapolated");
    }
  } 
  gSystem->MakeDirectory( "plots" );
  
  const unsigned int cosmicsIdx = cols.size()-1;
  
  std::cout << std::fixed;
  std::cout << std::setprecision(2);

  // 2020 Best Fit
  const double dCP    = .82;
  const double ssth23 = .568;
  const double ss2th13 = .085;
  const double dmsq32 = 2.406e-3;

  osc::IOscCalcAdjustable * calcosc= DefaultOscCalc();
  calcosc->SetdCP   ( dCP*(TMath::Pi()) );
  calcosc->SetTh23  ( asin(sqrt(ssth23)) );
  calcosc->SetDmsq32( dmsq32 );
  calcosc->SetTh13  ( .5*asin(sqrt(ss2th13)) );
  
  std::vector<std::map<std::string,std::pair<double,double>>> tots, sigs, bkgs;
  std::vector<std::pair<double,double>> totQuad, sigQuad, bkgQuad;

  for(unsigned int ext = 0; ext < extrapStr.size(); ++ext){
    std::vector< std::vector<Spectrum> > noms;
    std::vector<const IPrediction*> pred;
    std::vector< std::pair <Spectrum*, double> > cosmics;
    if (IsNuMu) {
      pred = GetNumuPredictions2020( 4, extrapStr[ext], calcosc, true, beam );
      cosmics = GetNumuCosmics2020( 4, beam );
    } else {
      pred.push_back( GetNuePrediction2020( extrapStr[ext], calcosc, true, beam, false ) );
      cosmics.push_back( GetNueCosmics2020( beam ) );
    }
    systs[ext].push_back(&kCosmicBkgScaleSyst);
    
    for (size_t pr=0; pr<pred.size(); ++pr) {
      std::vector<Spectrum> Tempnom;
      for(ColumnDef col: cols) {
        if(col.name == "Cosmics") 
        {
          Tempnom.push_back(*cosmics[pr].first);
          std::cerr << "Cosmics integral is " << Tempnom.back().Integral(livetime, 0, kLivetime) << std::endl;
        }
        else 
        {
          Tempnom.push_back(pred[pr]->PredictComponent(calcosc, col.flav, col.curr, col.sign));
          std::cerr << "Tempnom integral is " << Tempnom.back().Integral(pot) << std::endl;
        }
      }
      noms.push_back(Tempnom);
    }
    
    // Sum of numu quartiles
    if (IsNuMu)
    {
      std::cerr << "Summing numu nominal quantiles!" << std::endl;
      std::vector<Spectrum> Tempnom;
      for(ColumnDef col: cols) {
        if(col.name == "Cosmics") 
        {
          Spectrum sumCosmics = *cosmics[0].first + *cosmics[1].first + *cosmics[2].first + *cosmics[3].first;
          Tempnom.push_back(sumCosmics);
          std::cerr << "Cosmics integral is " << Tempnom.back().Integral(livetime, 0, kLivetime) << std::endl;
        }
        else 
        {
          Spectrum sumBeam = pred[0]->PredictComponent(calcosc, col.flav, col.curr, col.sign) +
                             pred[1]->PredictComponent(calcosc, col.flav, col.curr, col.sign) +
                             pred[2]->PredictComponent(calcosc, col.flav, col.curr, col.sign) +
                             pred[3]->PredictComponent(calcosc, col.flav, col.curr, col.sign);
          Tempnom.push_back(sumBeam);
          std::cerr << "Tempnom integral is " << Tempnom.back().Integral(pot) << std::endl;
        }
      }
      noms.push_back(Tempnom);
    }
    
    documentHeader(IsNuMu);
    std::cout<<"\\section{Best Fit Point}"<<std::endl;
    std::cout<<"$\\delta / \\pi = " << dCP << "$\\\\" << std::endl;
    std::cout<<"$\\sin^{2} \\theta_{23} = " << ssth23 << "$\\\\" << std::endl;
    std::cout<<"$\\Delta m_{32}^{2} = " << dmsq32 << "$\\\\" << std::endl;
    std::cout<<"$\\sin^{2} 2\\theta_{13} = " << ss2th13 << "$\\\\"<<std::endl;
    std::cout<<"\\pagebreak"<<std::endl;

    unsigned int predSize = pred.size();
    std::cerr << "predSize = " << predSize << std::endl;
    if(IsNuMu) predSize += 1; // for summed quartiles prediction, add one more summing step
    if(!IsNuMu && quantsComp) predSize += 3; // for low/high PID and peripheral bin, mimic the three additional bins
    // --- And now I want to loop through my vectors...
    for (unsigned int prIdx=0; prIdx<predSize; ++prIdx) {
      unsigned int pr = prIdx;
      if(!IsNuMu && quantsComp) pr = 0;
      
      TH1* hNomTot = noms[pr][0].ToTH1(pot);
      // Add cosmics
      hNomTot->Add(noms[pr][cosmicsIdx].ToTH1(livetime, kLivetime));
      
      TH1* hNomSig = noms[pr][1].ToTH1(pot);
      
      if(IsNuMu || prIdx==0)
      {
        // Total bkg - signal
        noms[pr][2] -= noms[pr][1];
        // Beam bkg - signal
        noms[pr][3] -= noms[pr][1];
      }
      
      TH1* hNomBkg = noms[pr][2].ToTH1(pot);
      hNomBkg->Add(noms[pr][cosmicsIdx].ToTH1(livetime, kLivetime));
      
      std::map<std::string, TH1*> hUpsTot, hUpsSig, hUpsBkg;
      std::map<std::string, TH1*> hDnsTot, hDnsSig, hDnsBkg;
      
      // need to add cosmic to total prediction and bkg prediction
      double totErr = 100/sqrt(Integral(noms[pr][0], pot)+Integral(noms[pr][cosmicsIdx], livetime, 0, kLivetime));
      double sigErr = 100/sqrt(Integral(noms[pr][1], pot));
      double bkgErr = 100/sqrt(Integral(noms[pr][2], pot)+Integral(noms[pr][cosmicsIdx], livetime, 0, kLivetime));
    
      std::string sSamp = beam + " " + sNuMu;

      std::cout<<"\n\n"<<std::endl;
      // -- Table Header.
      if (IsNuMu) {
        if(pr==4) {
          sSamp += " All Quartiles";
          std::cout<< "\\section{All Quartiles}"<<std::endl;
        }
        else {
          sSamp += " Quartile " + std::to_string(pr+1);
          std::cout<< "\\section{Ehad Quartile "<<std::to_string(pr+1)<<"}"<<std::endl;
        }
        std::cout<<"\\hspace*{-1.5cm}\\resizebox{1.25\\textwidth}{!}{"<<std::endl;
        std::cout<<"\\small"<<std::endl;
        std::cout<<"\\begin{tabular}{l c c c c c c c c c c c c c c c c c c}"<<std::endl;
        std::cout<<"\\multicolumn{19}{c}{\\bf{Systematic Summary for " << sSamp << " (\\%)}} \\\\"<<std::endl;
      } else {
        std::string sBin = "All Prediction";
        if(prIdx==1) sBin = "Low PID";
        else if(prIdx==2) sBin = "High PID";
        else if(prIdx==3) sBin = "Peripheral";
        sSamp += " "+sBin;
        std::cout<<"\\section{"<< sBin <<"}"<<std::endl;
        std::cout<<"\\hspace*{-1.5cm}\\resizebox{1.25\\textwidth}{!}{"<<std::endl;
        std::cout<<"\\small"<<std::endl;
        std::cout<<"\\begin{tabular}{l c c c c c c c c c c c c c c c c c c c c}"<<std::endl;
        std::cout<<"\\multicolumn{21}{c}{\\bf{Systematic Summary for " << sSamp << " (\\%)}} \\\\"<<std::endl;
      }
      // -- Top line of table
      std::cout << "{Systematic}";
      for(ColumnDef col: cols) std::cout << " & \\multicolumn{2}{c}{" << col.name << "}";
      std::cout << " \\\\" << std::endl;
      std::cout << "\\hline" << std::endl;

      std::pair<double,double> quads[cols.size()];
      for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx) quads[colIdx] = {0,0};

      std::map<std::string, std::pair<double,double>> barsTot, barsSig, barsBkg;
      
      for(const ISyst* syst: systs[ext]){

        std::cout << FixSystName(syst->ShortName());
        
        std::vector<Spectrum> shiftsUp, shiftsDn;
        for(ColumnDef col: cols){
          // not summing the quantiles (nue or not the last pr)
          if(!IsNuMu || pr != pred.size())
          {
            if(col.name != "Cosmics")
            {
              if(syst == &kCosmicBkgScaleSyst)
              {
                shiftsUp.push_back(pred[pr]->PredictComponent(calcosc, col.flav, col.curr, col.sign));
                shiftsDn.push_back(pred[pr]->PredictComponent(calcosc, col.flav, col.curr, col.sign));
              }
              else
              {
                shiftsUp.push_back(pred[pr]->PredictComponentSyst(calcosc, SystShifts(syst, +1), col.flav, col.curr, col.sign));
                shiftsDn.push_back(pred[pr]->PredictComponentSyst(calcosc, SystShifts(syst, -1), col.flav, col.curr, col.sign));
              }
            }
            else
            {
              shiftsUp.push_back(*cosmics[pr].first);
              shiftsDn.push_back(*cosmics[pr].first);
          
              if(syst == &kCosmicBkgScaleSyst)
              {
                shiftsUp.back() = *ShiftedCosmics(cosmics[pr].first, +1);
                shiftsDn.back() = *ShiftedCosmics(cosmics[pr].first, -1);
              }
            }   
          }
          // summing the numu quantiles
          else
          {
            if(col.name != "Cosmics")
            {
              if(syst == &kCosmicBkgScaleSyst)
              {
                shiftsUp.push_back(pred[0]->PredictComponent(calcosc, col.flav, col.curr, col.sign)+
                                   pred[1]->PredictComponent(calcosc, col.flav, col.curr, col.sign)+
                                   pred[2]->PredictComponent(calcosc, col.flav, col.curr, col.sign)+
                                   pred[3]->PredictComponent(calcosc, col.flav, col.curr, col.sign));
                shiftsDn.push_back(pred[0]->PredictComponent(calcosc, col.flav, col.curr, col.sign)+
                                   pred[1]->PredictComponent(calcosc, col.flav, col.curr, col.sign)+
                                   pred[2]->PredictComponent(calcosc, col.flav, col.curr, col.sign)+
                                   pred[3]->PredictComponent(calcosc, col.flav, col.curr, col.sign));
              }
              else
              {
                shiftsUp.push_back(pred[0]->PredictComponentSyst(calcosc, SystShifts(syst, +1), col.flav, col.curr, col.sign)+
                                   pred[1]->PredictComponentSyst(calcosc, SystShifts(syst, +1), col.flav, col.curr, col.sign)+
                                   pred[2]->PredictComponentSyst(calcosc, SystShifts(syst, +1), col.flav, col.curr, col.sign)+
                                   pred[3]->PredictComponentSyst(calcosc, SystShifts(syst, +1), col.flav, col.curr, col.sign));
                shiftsDn.push_back(pred[0]->PredictComponentSyst(calcosc, SystShifts(syst, -1), col.flav, col.curr, col.sign)+
                                   pred[1]->PredictComponentSyst(calcosc, SystShifts(syst, -1), col.flav, col.curr, col.sign)+
                                   pred[2]->PredictComponentSyst(calcosc, SystShifts(syst, -1), col.flav, col.curr, col.sign)+
                                   pred[3]->PredictComponentSyst(calcosc, SystShifts(syst, -1), col.flav, col.curr, col.sign));
              }
            }
            else
            {          
              if(syst == &kCosmicBkgScaleSyst)
              {
                Spectrum cosmicsUp = *ShiftedCosmics(cosmics[0].first, +1) + 
                                     *ShiftedCosmics(cosmics[1].first, +1) + 
                                     *ShiftedCosmics(cosmics[2].first, +1) + 
                                     *ShiftedCosmics(cosmics[3].first, +1);
                
                Spectrum cosmicsDn = *ShiftedCosmics(cosmics[0].first, -1) + 
                                     *ShiftedCosmics(cosmics[1].first, -1) + 
                                     *ShiftedCosmics(cosmics[2].first, -1) + 
                                     *ShiftedCosmics(cosmics[3].first, -1);
                                    
                shiftsUp.push_back(cosmicsUp);
                shiftsDn.push_back(cosmicsDn);
              }
              else
              {
                shiftsUp.push_back(*cosmics[0].first+*cosmics[1].first+*cosmics[2].first+*cosmics[3].first);
                shiftsDn.push_back(*cosmics[0].first+*cosmics[1].first+*cosmics[2].first+*cosmics[3].first);
              }
            }
          }
        }
        
        hUpsTot[ (std::string)FixSystName(syst->ShortName()) ] = shiftsUp[0].ToTH1(pot);
        hDnsTot[ (std::string)FixSystName(syst->ShortName()) ] = shiftsDn[0].ToTH1(pot);
        hUpsSig[ (std::string)FixSystName(syst->ShortName()) ] = shiftsUp[1].ToTH1(pot);
        hDnsSig[ (std::string)FixSystName(syst->ShortName()) ] = shiftsDn[1].ToTH1(pot);
        hUpsBkg[ (std::string)FixSystName(syst->ShortName()) ] = (shiftsUp[0] - shiftsUp[1]).ToTH1(pot);
        hDnsBkg[ (std::string)FixSystName(syst->ShortName()) ] = (shiftsDn[0] - shiftsDn[1]).ToTH1(pot);
        if(syst == &kCosmicBkgScaleSyst)
        {
          TH1D* cosUp;
          TH1D* cosDn;
          if(IsNuMu && pr==4) 
          {
            cosUp = ShiftedCosmics(cosmics[0].first,+1)->ToTH1(livetime,kLivetime);
            cosUp -> Add(ShiftedCosmics(cosmics[1].first,+1)->ToTH1(livetime,kLivetime));
            cosUp -> Add(ShiftedCosmics(cosmics[2].first,+1)->ToTH1(livetime,kLivetime));
            cosUp -> Add(ShiftedCosmics(cosmics[3].first,+1)->ToTH1(livetime,kLivetime));
            
            cosDn = ShiftedCosmics(cosmics[0].first,-1)->ToTH1(livetime,kLivetime);
            cosDn -> Add(ShiftedCosmics(cosmics[1].first,-1)->ToTH1(livetime,kLivetime));
            cosDn -> Add(ShiftedCosmics(cosmics[2].first,-1)->ToTH1(livetime,kLivetime));
            cosDn -> Add(ShiftedCosmics(cosmics[3].first,-1)->ToTH1(livetime,kLivetime));
          } else {
            cosUp = ShiftedCosmics(cosmics[pr].first,+1)->ToTH1(livetime,kLivetime);
            cosDn = ShiftedCosmics(cosmics[pr].first,-1)->ToTH1(livetime,kLivetime);
          }
          hUpsTot[ (std::string)FixSystName(syst->ShortName()) ]->Add(cosUp);
          hDnsTot[ (std::string)FixSystName(syst->ShortName()) ]->Add(cosDn);
          hUpsBkg[ (std::string)FixSystName(syst->ShortName()) ]->Add(cosUp);
          hDnsBkg[ (std::string)FixSystName(syst->ShortName()) ]->Add(cosDn);
        }
        else
        {
          hUpsTot[ (std::string)FixSystName(syst->ShortName()) ]->Add(noms[pr][cosmicsIdx].ToTH1(livetime, kLivetime));
          hDnsTot[ (std::string)FixSystName(syst->ShortName()) ]->Add(noms[pr][cosmicsIdx].ToTH1(livetime, kLivetime));
          hUpsBkg[ (std::string)FixSystName(syst->ShortName()) ]->Add(noms[pr][cosmicsIdx].ToTH1(livetime, kLivetime));
          hDnsBkg[ (std::string)FixSystName(syst->ShortName()) ]->Add(noms[pr][cosmicsIdx].ToTH1(livetime, kLivetime));
        }
        
        // Shifted Total bkg = Total pred with cosmics - signal
        shiftsUp[2] -= shiftsUp[1];
        shiftsDn[2] -= shiftsDn[1];
      
        // Shifted Beam bkg = Total pred - signal
        shiftsUp[3] -= shiftsUp[1];
        shiftsDn[3] -= shiftsDn[1];
      
        for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx){
              const Spectrum& nom = noms[pr][colIdx];
          const Spectrum& cosNom = noms[pr][cosmicsIdx];
        
          double ratioUp, errUp, ratioDn, errDn;
          if(colIdx == 0 || colIdx == 2)
          {
            ratioUp = 1.0*( Integral(shiftsUp[colIdx], pot) + Integral(shiftsUp[cosmicsIdx], livetime, 0, kLivetime) ) / (1.0*(Integral(nom, pot)+Integral(cosNom,livetime,0,kLivetime)));
            errUp = 100.0*(ratioUp-1.0);
          
            ratioDn = 1.0*( Integral(shiftsDn[colIdx], pot) + Integral(shiftsDn[cosmicsIdx], livetime, 0, kLivetime) ) / (1.0*(Integral(nom, pot)+Integral(cosNom,livetime,0,kLivetime)));
            errDn = 100.0*(ratioDn-1.0);
          }
          else if(colIdx == cosmicsIdx)
          {
            ratioUp = 1.0*( Integral(shiftsUp[colIdx],livetime,0,kLivetime)) / (1.0*(Integral(nom,livetime,0,kLivetime)));
            errUp = 100.0*(ratioUp-1.0);
          
            ratioDn = 1.0*( Integral(shiftsDn[colIdx],livetime,0,kLivetime)) / (1.0*(Integral(nom,livetime,0,kLivetime)));
            errDn = 100.0*(ratioDn-1.0);
          }
          else
          {
                ratioUp = 1.0*Integral(shiftsUp[colIdx], pot) / (1.0*Integral(nom, pot));
                errUp = 100.0*(ratioUp-1.0);
      
                ratioDn = 1.0*Integral(shiftsDn[colIdx], pot) / (1.0*Integral(nom, pot));
                errDn = 100.0*(ratioDn-1.0);
          }
          
          // if investigating the individual nue bins, we need to integrate over the particular bins of the predicted spectrum
          if(!IsNuMu && quantsComp && prIdx != 0)
          {
            std::pair<int,int> Int = {1,9}; // low PID
            if(prIdx==2) Int = {10,18}; // high PID
            if(prIdx==3) Int = {19,hNomTot->GetNbinsX()}; // peripheral
             
            if(colIdx == 0 || colIdx == 2)
            {
              ratioUp = 1.0*( shiftsUp[colIdx].ToTH1(pot)->Integral(Int.first,Int.second) + shiftsUp[cosmicsIdx].ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.ToTH1(pot)->Integral(Int.first,Int.second) + cosNom.ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) );
              errUp = 100.0*(ratioUp-1.0);
              
              ratioDn = 1.0*( shiftsDn[colIdx].ToTH1(pot)->Integral(Int.first,Int.second) + shiftsDn[cosmicsIdx].ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.ToTH1(pot)->Integral(Int.first,Int.second) + cosNom.ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) );
              errDn = 100.0*(ratioDn-1.0);
            }
            else if(colIdx == cosmicsIdx)
            {
              ratioUp = 1.0*( shiftsUp[colIdx].ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) );
              errUp = 100.0*(ratioUp-1.0);
          
              ratioDn = 1.0*( shiftsDn[colIdx].ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) );
              errDn = 100.0*(ratioDn-1.0);
            }
            else
            {
                  ratioUp = 1.0*( shiftsUp[colIdx].ToTH1(pot)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.ToTH1(pot)->Integral(Int.first,Int.second) ) );
                  errUp = 100.0*(ratioUp-1.0);
      
                  ratioDn = 1.0*( shiftsDn[colIdx].ToTH1(pot)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.ToTH1(pot)->Integral(Int.first,Int.second) ) );
                  errDn = 100.0*(ratioDn-1.0);
            }
          }
        
              double min = std::min(errUp,errDn);
              double max = std::max(errUp,errDn);
              if(max <  0.01) max = 0;
              if(min > -0.01) min = 0;
      
              quads[colIdx].first += min*min;
              quads[colIdx].second += max*max;
      
        
              std::string minstr = "white";
              std::string maxstr = "white";
          /*
              if( fabs(min) >= 1 && fabs(min) < 2) minstr="Apricot";
              if( fabs(min) >= 2 && fabs(min) < 3) minstr="YellowOrange";
              if( fabs(min) >= 3                 ) minstr="RedOrange";

              if( fabs(max) >= 1 && fabs(max) < 2) maxstr="Apricot";
              if( fabs(max) >= 2 && fabs(max) < 3) maxstr="YellowOrange";
              if( fabs(max) >= 3                 ) maxstr="RedOrange";
          */

          std::cout << " & \\cellcolor{" << minstr << "}" << (min == 0 ? "-" : "") << min << " & \\cellcolor{"<< maxstr << "}+" << max;

              min = fabs(min);
              max = fabs(max);
              if(cols[colIdx].name == "$\\nu_e$ Signal" || 
                 cols[colIdx].name == "$\\bar\\nu_e$ Signal" || 
                 cols[colIdx].name == "$\\nu_\\mu$ Signal" ||
                 cols[colIdx].name == "$\\bar\\nu_\\mu$ Signal")
                barsSig[ (std::string)FixSystName(syst->ShortName()) ] = {min,max};
            
          if(cols[colIdx].name == "Total Pred")
            barsTot[ (std::string)FixSystName(syst->ShortName()) ] = {min,max};
            
              if(cols[colIdx].name == "Total Bkg")
                barsBkg[ (std::string)FixSystName(syst->ShortName()) ] = {min,max};
        } // end for colIdx

        std::cout << " \\\\" << std::endl;
      } // end for syst

      std::cout << "\\hline" << std::endl;

      std::cout << "Quadrature";
      for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx) std::cout << " & -" << sqrt(quads[colIdx].first) << " & +" << sqrt(quads[colIdx].second);
      std::cout << std::endl;
      std::cout<<"\\end{tabular}"<<std::endl;
      std::cout<<"}"<<std::endl;

      if(sum){
        barsTot = SumSysts(barsTot);
        barsSig = SumSysts(barsSig);
        barsBkg = SumSysts(barsBkg);

        std::cout<<std::endl;
        std::cout<<"\\vfill" << std::endl;
        std::cout<<"Summing..."<<std::endl;
        std::cout<<std::endl;

        std::cout<<"\\begin{tabular}{l c c c c c c}"<<std::endl;
        std::cout<<"\\multicolumn{7}{c}{\\bf{Summed Systematic for " << sSamp << " (\\%)}} \\\\"<<std::endl;
        std::cout<<"{Systematic}";
        for(int colIdx = 0; colIdx < 3; ++colIdx) std::cout << " & \\multicolumn{2}{c}{" << cols[colIdx].name << "}";
        std::cout<<" \\\\"<<std::endl;
        std::cout << "\\hline" << std::endl;

        for(int i = 0;i < kNumSystTypes; ++i)
        {
          std::cout << systTypes[i] << " & "
            << "-" << barsTot[systTypes[i]].first << " & +" << barsTot[systTypes[i]].second << " & "
                    << "-" << barsSig[systTypes[i]].first << " & +" << barsSig[systTypes[i]].second << " & "
                    << "-" << barsBkg[systTypes[i]].first << " & +" << barsBkg[systTypes[i]].second <<" \\\\"
                    << std::endl;
        }
    
        std::cout << "\\hline" << std::endl;
        std::cout << "Quadrature";
        for(int colIdx = 0; colIdx < 3; ++colIdx) std::cout << " & -" << sqrt(quads[colIdx].first) << " & +" << sqrt(quads[colIdx].second);
        std::cout<<std::endl;
        std::cout<<"\\end{tabular}"<<std::endl;
        std::cout<<std::endl;
      }
      
      // tots, sigs and bkg have following struct:
      // nue: extrap1_all, extrap1_lowPID, extrap1_highPID, extrap1_peripheral, extrap2_all...
      // numu: extrap1_Q1, extrap1_Q2, extrap1_Q3, extrap1_Q4, extrap1_allQuants, extrap2_Q1...
      tots.push_back(barsTot);
      sigs.push_back(barsSig);
      bkgs.push_back(barsBkg);
      
      for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx){
        quads[colIdx].first = sqrt(quads[colIdx].first);
        quads[colIdx].second = sqrt(quads[colIdx].second);
      }

      totQuad.push_back(quads[0]);
      sigQuad.push_back(quads[1]);
      bkgQuad.push_back(quads[2]);

      std::string  sumStr  = sum  == true  ? "short" : "full";
      std::string  sQuant  = "";
      if (IsNuMu) { 
        sQuant  = "_Q"+std::to_string(pr+1); 
        if(pr==pred.size()) sQuant = "_AllQuartiles";
      }
      else
      {
        sQuant = "_AllPrediction";
        if(prIdx==1) sQuant = "_lowPID";
        else if(prIdx==2) sQuant = "_highPID";
        else if(prIdx==3) sQuant = "_peripheral";
      }
      
      std::cout << "\n\n\n And now to save things..." << std::endl;
      
      sampLab = "#nu_{e} Selection";
      if(beam == "rhc") sampLab = "#bar{#nu}_{e} Selection";
      if(IsNuMu)
      {
        std::string quart = "All Quartiles";
        if(pr<4) quart ="Quartile "+std::to_string(pr+1);
        sampLab = "#nu_{#mu} Selection "+quart;
        if(beam == "rhc") sampLab = "#bar{#nu}_{#mu} Selection "+quart;
      }
      else
      {
         if(prIdx==1) sampLab += " Low PID";
         if(prIdx==2) sampLab += " High PID";
         if(prIdx==3) sampLab += " Peripheral";
      }
      
      // Energy prediction ratios
      TH1* hD = (TH1*)hNomTot->Clone();
      TH1* hNomTotRatio = (TH1*)hNomTot->Clone();
      hNomTotRatio->Divide(hNomTot);
      std::vector<TH1*> tempUp, tempDn;
      for(auto it: hUpsTot) tempUp.push_back(it.second);
      for(auto it: hDnsTot) tempDn.push_back(it.second);
      
      std::vector<TH1*> hUpsTotRatios = MakeRatios(hNomTot, tempUp);
      std::vector<TH1*> hDnsTotRatios = MakeRatios(hNomTot, tempDn);
      
      std::map<std::string, std::vector<TH1*>> hUpsTotType = SortSystHists(hUpsTot);
      std::map<std::string, std::vector<TH1*>> hDnsTotType = SortSystHists(hDnsTot);
      
      std::vector<int> types = {0,1,2,4,5};
      if(IsNuMu) types = {1,2,4,5,6};
      TLegend *lTypes = new TLegend(.1,.1,.9,.9);
      lTypes->SetFillStyle(0);
      
      TCanvas *c = new TCanvas(TString("c")+TString(sQuant),"c");
      MyPlotWithSystErrorBand(hNomTotRatio, hUpsTotRatios, hDnsTotRatios, kBlack, kTotalMCErrorBandColor, 2, true);
      hD->SetLineColor(kTotalMCErrorBandColor);
      lTypes->AddEntry(hD,"Total syst. uncertainty","l");
      
      for(auto j: types)
      {
        std::vector<TH1*> hUpsRatios = MakeRatios(hNomTot, hUpsTotType[systTypes[j]]);
        std::vector<TH1*> hDnsRatios = MakeRatios(hNomTot, hDnsTotType[systTypes[j]]);
        
        MyPlotWithSystErrorBand(hNomTotRatio, hUpsRatios, hDnsRatios, kBlack, systCols[j], 1, false);
        TH1* hDummy = (TH1*)hNomTot->Clone();
        hDummy->SetLineColor(systCols[j]);
        lTypes->AddEntry(hDummy,systTypes[j].c_str(),"l");
      }
      MyPlotWithSystErrorBand(hNomTotRatio, hUpsTotRatios, hDnsTotRatios, kBlack, kTotalMCErrorBandColor, 2, false);
      drawBeamLabel(beamLab,0.1);
      drawSampleLabel(sampLab);
      preliminary();
      gPad->Print(("plots/ratioband_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_tot.pdf").c_str());
      
      TCanvas *cLeg = new TCanvas(TString("cLeg")+TString(sQuant),"cLeg");
      
      lTypes->Draw();
      gPad->Print(("plots/ratioband_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_legend.pdf").c_str());
      
      // Classic bar-plots
      TCanvas *cTot = new TCanvas(TString("cTot")+TString(sQuant),"cTot");
      ErrorBarChart(barsTot,{totErr,totErr},"Total Prediction Uncertainty (%)");
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();
      gPad->Print(("plots/syst_summary_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_tot.pdf").c_str());

      TCanvas *cTot2 = new TCanvas(TString("cTot2")+TString(sQuant),"cTot2");
      ErrorBarChart(barsTot,{0,0},"Total Prediction Uncertainty (%)");
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();
      gPad->Print(("plots/syst_summary_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_tot_nostat.pdf").c_str());
      
      TCanvas *cSig = new TCanvas(TString("cSig")+TString(sQuant),"cSig");
      ErrorBarChart(barsSig,{sigErr,sigErr},"Signal Uncertainty (%)");
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();
      gPad->Print(("plots/syst_summary_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_sig.pdf").c_str());

      TCanvas *cSig2 = new TCanvas(TString("cSig2")+TString(sQuant),"cSig2");
      ErrorBarChart(barsSig,{0,0},"Signal Uncertainty (%)");
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();
      gPad->Print(("plots/syst_summary_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_sig_nostat.pdf").c_str());

      TCanvas *cBkg = new TCanvas(TString("cBkg")+TString(sQuant),"cBkg");
      ErrorBarChart(barsBkg,{bkgErr,bkgErr},"Background Uncertainty (%)");
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();
      gPad->Print(("plots/syst_summary_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_bkg.pdf").c_str());

      TCanvas *cBkg2 = new TCanvas(TString("cBkg2")+TString(sQuant),"cBkg2");
      ErrorBarChart(barsBkg,{0,0},"Background Uncertainty (%)");
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();
      gPad->Print(("plots/syst_summary_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_bkg_nostat.pdf").c_str());
    } // Pred loop
    std::cout<<"\n"<<std::endl;

    if (!IsNuMu) std::cout<<"\\end{landscape}"<<std::endl;
    std::cout<<"\\end{document}"<<std::endl;
  } // extrapStr loop
  
  // numu quartiles comparison
  if(IsNuMu && quantsComp)
  {
    unsigned int nQuants = 4; // 4
    
    for(unsigned int i = 0; i < extrapStr.size(); ++i)
    {
      std::vector<std::string> quantLab = {"Quartile 1", "Quartile 2", "Quartile 3", "Quartile 4"};
      std::string sumStr = "full";
      std::vector<std::string> labels;
      std::vector<std::vector<std::pair<double,double>>> fTot;
      std::vector<std::vector<std::pair<double,double>>> fSig;
      std::vector<std::vector<std::pair<double,double>>> fBkg;
      std::vector<std::pair<double,double>> fTotQuad;
      std::vector<std::pair<double,double>> fSigQuad;
      std::vector<std::pair<double,double>> fBkgQuad;
      
      for(unsigned int q = 0; q < nQuants; ++q)
      {
        std::vector<std::pair<double,double>> tempTot;
        std::vector<std::pair<double,double>> tempSig;
        std::vector<std::pair<double,double>> tempBkg;
        int idx = i*(nQuants+1) + q;
        std::cerr << "Quantile idx: " << idx << " " << extrapStr[i] << std::endl;
        
        if(sum)
        {
          sumStr = "short";
          for(int j = 0; j < kNumSystTypes; ++j){
                if(q == 0) labels.push_back(systTypes[j]);
            tempTot.push_back(tots[idx][systTypes[j]]);
                tempSig.push_back(sigs[idx][systTypes[j]]);
                tempBkg.push_back(bkgs[idx][systTypes[j]]);
          }
        }
        else
        {
          for(const ISyst* syst: systs[0]){
            if(q == 0) labels.push_back( (std::string)FixSystName(syst->ShortName()) );
            tempTot.push_back(tots[idx][ (std::string)FixSystName(syst->ShortName()) ]);
                tempSig.push_back(sigs[idx][ (std::string)FixSystName(syst->ShortName()) ]);
                tempBkg.push_back(bkgs[idx][ (std::string)FixSystName(syst->ShortName()) ]);
          }
        }

        
        fTot.push_back(tempTot);
        fSig.push_back(tempSig);
        fBkg.push_back(tempBkg);
        fTotQuad.push_back(totQuad[idx]);
        fSigQuad.push_back(sigQuad[idx]);
        fBkgQuad.push_back(bkgQuad[idx]);
      }
      
      sampLab = "#nu_{#mu} Selection";
      if(beam == "rhc") sampLab = "#bar{#nu}_{#mu} Selection";
      
      std::vector<int> colors = {kRed-7,kPink+6,kMagenta-6,kViolet-3};
      new TCanvas;
      myBarChart("Total Prediction Uncertainty (%)", labels,
               fTot, fTotQuad, quantLab, colors);
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();
  
      gPad->Print(("plots/syst_quant_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_tot_"+extrapStr[i]+".pdf").c_str());
      
      new TCanvas;
      myBarChart("Signal Uncertainty (%)", labels,
               fSig, fSigQuad, quantLab, colors);
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();
  
      gPad->Print(("plots/syst_quant_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_sig_"+extrapStr[i]+".pdf").c_str());
  
      new TCanvas;
      myBarChart("Background Uncertainty (%)", labels,
               fBkg, fBkgQuad, quantLab, colors);
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();

      gPad->Print(("plots/syst_quant_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_bkg_"+extrapStr[i]+".pdf").c_str());
    }
  }
  
  // nue bins comparison
  if(!IsNuMu && quantsComp)
  {
    unsigned int nBins = 3; // 4
    
    for(unsigned int i = 0; i < extrapStr.size(); ++i)
    {
      std::vector<std::string> quantLab = {"Low PID", "High PID", "Peripheral"};
      std::string sumStr = "full";
      std::vector<std::string> labels;
      std::vector<std::vector<std::pair<double,double>>> fTot;
      std::vector<std::vector<std::pair<double,double>>> fSig;
      std::vector<std::vector<std::pair<double,double>>> fBkg;
      std::vector<std::pair<double,double>> fTotQuad;
      std::vector<std::pair<double,double>> fSigQuad;
      std::vector<std::pair<double,double>> fBkgQuad;
      
      for(unsigned int b = 1; b <= nBins; ++b)
      {
        std::vector<std::pair<double,double>> tempTot;
        std::vector<std::pair<double,double>> tempSig;
        std::vector<std::pair<double,double>> tempBkg;
        int idx = i*(nBins+1) + b;
        std::cerr << "Bin idx: " << idx << " " << extrapStr[i] << std::endl;
        
        if(sum)
        {
          sumStr = "short";
          for(int j = 0; j < kNumSystTypes; ++j){
                if(b == 1) labels.push_back(systTypes[j]);
            tempTot.push_back(tots[idx][systTypes[j]]);
                tempSig.push_back(sigs[idx][systTypes[j]]);
                tempBkg.push_back(bkgs[idx][systTypes[j]]);
          }
        }
        else
        {
          for(const ISyst* syst: systs[0]){
            if(b == 1) labels.push_back( (std::string)FixSystName(syst->ShortName()) );
            tempTot.push_back(tots[idx][ (std::string)FixSystName(syst->ShortName()) ]);
                tempSig.push_back(sigs[idx][ (std::string)FixSystName(syst->ShortName()) ]);
                tempBkg.push_back(bkgs[idx][ (std::string)FixSystName(syst->ShortName()) ]);
          }
        }

        
        fTot.push_back(tempTot);
        fSig.push_back(tempSig);
        fBkg.push_back(tempBkg);
        fTotQuad.push_back(totQuad[idx]);
        fSigQuad.push_back(sigQuad[idx]);
        fBkgQuad.push_back(bkgQuad[idx]);
      }
      
      sampLab = "#nu_{e} Selection";
      if(beam == "rhc") sampLab = "#bar{#nu}_{e} Selection";
      
      std::vector<int> colors = {kViolet-4,kRed-7,kViolet+8};
      new TCanvas;
      myBarChart("Total Prediction Uncertainty (%)", labels,
               fTot, fTotQuad, quantLab, colors);
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();
  
      gPad->Print(("plots/syst_bins_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_tot_"+extrapStr[i]+".pdf").c_str());
      
      new TCanvas;
      myBarChart("Signal Uncertainty (%)", labels,
               fSig, fSigQuad, quantLab, colors);
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();
  
      gPad->Print(("plots/syst_bins_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_sig_"+extrapStr[i]+".pdf").c_str());
  
      new TCanvas;
      myBarChart("Background Uncertainty (%)", labels,
               fBkg, fBkgQuad, quantLab, colors);
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();

      gPad->Print(("plots/syst_bins_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_bkg_"+extrapStr[i]+".pdf").c_str());
    }
  }
  
  if(extrapComp)
  {
    unsigned int nQuants = 1; //nue
    if(!IsNuMu && quantsComp) nQuants = 4; // low, high PID, peripheral
    if(IsNuMu) nQuants = 5; // 4 + summed

    for(unsigned int q = 0; q < nQuants; ++q)
    {
      std::string sumStr = "full";
      std::vector<std::string> labels;
      std::vector<std::vector<std::pair<double,double>>> fTot;
      std::vector<std::vector<std::pair<double,double>>> fSig;
      std::vector<std::vector<std::pair<double,double>>> fBkg;
      std::vector<std::pair<double,double>> fTotQuad;
      std::vector<std::pair<double,double>> fSigQuad;
      std::vector<std::pair<double,double>> fBkgQuad;
    
      for(unsigned int i = 0; i < extrapStr.size(); ++i)
      {
        std::vector<std::pair<double,double>> tempTot;
        std::vector<std::pair<double,double>> tempSig;
        std::vector<std::pair<double,double>> tempBkg;
        int idx = i*nQuants + q;
        std::cerr << "Quartile idx: " << idx << " " << extrapStr[i] << std::endl;

        if(sum)
        {
          sumStr = "short";
          for(int j = 0; j < kNumSystTypes; ++j){
                if(i == 0) labels.push_back(systTypes[j]);
            tempTot.push_back(tots[idx][systTypes[j]]);
                tempSig.push_back(sigs[idx][systTypes[j]]);
                tempBkg.push_back(bkgs[idx][systTypes[j]]);
          }
        }
        else
        {
          // just for now hardcoded because the noextrap and normal extrap does not have Lepton Angle Syst and different accept syst
          for(const ISyst* syst: systs[1]){
                if(i == 0) labels.push_back( (std::string)FixSystName(syst->ShortName()) );
            if ( i==0 && FixSystName(syst->ShortName()).Contains("LeptonAngle") )
            {
              tempTot.push_back({.0,.0});
              tempSig.push_back({.0,.0});
                  tempBkg.push_back({.0,.0});
            }
            else if ( i==0 && FixSystName(syst->ShortName()).Contains("accept") )
            {
              if (beam == "fhc") 
              {
                tempTot.push_back(tots[idx][ (std::string)"accept signalkin FHC 2020" ]);
                tempSig.push_back(sigs[idx][ (std::string)"accept signalkin FHC 2020" ]);
                tempBkg.push_back(bkgs[idx][ (std::string)"accept signalkin FHC 2020" ]);
              }
              else
              {
                tempTot.push_back(tots[idx][ (std::string)"accept signalkin RHC 2020" ]);
                tempSig.push_back(sigs[idx][ (std::string)"accept signalkin RHC 2020" ]);
                tempBkg.push_back(bkgs[idx][ (std::string)"accept signalkin RHC 2020" ]);
              }
            }
            else
            {
              tempTot.push_back(tots[idx][ (std::string)FixSystName(syst->ShortName()) ]);
                  tempSig.push_back(sigs[idx][ (std::string)FixSystName(syst->ShortName()) ]);
                  tempBkg.push_back(bkgs[idx][ (std::string)FixSystName(syst->ShortName()) ]);
            }
          }
        }

        fTot.push_back(tempTot);
        fSig.push_back(tempSig);
        fBkg.push_back(tempBkg);
        fTotQuad.push_back(totQuad[idx]);
        fSigQuad.push_back(sigQuad[idx]);
        fBkgQuad.push_back(bkgQuad[idx]);
      }
      
      std::string sQuant = "";
      if (IsNuMu) { 
        sQuant  = "_Q"+std::to_string(q+1); 
        if(q==4) sQuant = "_AllQuartiles";
      }
      else {
        sQuant = "_AllPrediction";
        if(q==1) sQuant = "_lowPID";
        if(q==2) sQuant = "_highPID";
        if(q==3) sQuant = "_peripheral";
      }
      
      sampLab = "#nu_{e} Selection";
      if(beam == "rhc") sampLab = "#bar{#nu}_{e} Selection";
      if(IsNuMu)
      {
        std::string quart = "All Quartiles";
        if(q<4) quart ="Quartile "+std::to_string(q+1);
        sampLab = "#nu_{#mu} Selection "+quart;
        if(beam == "rhc") sampLab = "#bar{#nu}_{#mu} Selection "+quart;
      }
      else
      {
        if(q==1) sampLab += " Low PID";
        if(q==2) sampLab += " High PID";
        if(q==3) sampLab += " Peripheral";
      }
      
      new TCanvas;
      myBarChart("Total Prediction Uncertainty (%)", labels,
               fTot, fTotQuad, extrapLab, {-1});
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();
  
      gPad->Print(("plots/syst_extrap_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_tot"+sQuant+".pdf").c_str());
      
      new TCanvas;
      myBarChart("Signal Uncertainty (%)", labels,
               fSig, fSigQuad, extrapLab, {-1});
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();
  
      gPad->Print(("plots/syst_extrap_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_sig"+sQuant+".pdf").c_str());
  
      new TCanvas;
      myBarChart("Background Uncertainty (%)", labels,
               fBkg, fBkgQuad, extrapLab, {-1});
      gPad->SetLeftMargin(.35);
      drawBeamLabel(beamLab);
      drawSampleLabel(sampLab);
      preliminary();

      gPad->Print(("plots/syst_extrap_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_bkg"+sQuant+".pdf").c_str());
    }
  }
  
}