resolution2018.C

Go to the documentation of this file.

#pragma once
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Analysis/Prod4Loaders.h"
#include "CAFAna/Prediction/PredictionGenerator.h"
#include "CAFAna/Decomp/CheatDecomp.h"
#include "CAFAna/Core/Binning.h"
#include "CAFAna/Cuts/Cuts.h"

#include "3FlavorAna/Cuts/QuantileCuts.h"

#include "CAFAna/Analysis/Exposures.h"
#include "3FlavorAna/Systs/EnergySysts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "3FlavorAna/Cuts/BeamNueCuts.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/Binnings.h"
#include "3FlavorAna/Vars/HistAxes.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Systs/Systs.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Prediction/PredictionScaleComp.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "CAFAna/Extrap/ModularExtrap.h"
#include "3FlavorAna/Decomp/BENDecomp.h"
#include "3FlavorAna/Decomp/MichelDecomp.h"
#include "3FlavorAna/Decomp/NumuDecomp.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TGaxis.h"
#include "TFile.h"
#include "TLegend.h"
#include "TLine.h"
#include "TStyle.h"
#include "TLatex.h"
#include "TROOT.h"
#include "TColor.h"
#include "THStack.h"
#include "TPaveText.h"
#include "TGraphErrors.h"
#include "TMultiGraph.h"

#include "Utilities/rootlogon.C"

using namespace ana;
#include "OscLib/OscCalcDumb.h"
#include "OscLib/IOscCalc.h"
#include <vector>

//void Simulation(bool normal = false)
//{
//  if(normal){
//    TLatex* prelim = new TLatex(.9, .95, "NOvA Simulation");
//    prelim->SetTextColor(kGray+1);
//    prelim->SetNDC();
//    prelim->SetTextSize(2/30.);
//    prelim->SetTextAlign(32);
//    prelim->Draw();
//  }
//  else{
//   TLatex* prelim = new TLatex(.95, .92, "NOvA Simulation");
//   prelim->SetTextColor(kGray+1);
//   prelim->SetNDC();
//   prelim->SetTextSize(2/20.);
//   prelim->SetTextAlign(32);
//   prelim->Draw();
//  }
//}

void Beam(bool isRHC){
  if(isRHC){
    TLatex* prelim = new TLatex(.276, .95, "Antineutrino Beam");
    prelim->SetTextColor(kBlack);
    prelim->SetNDC();
    prelim->SetTextSize(1/30.);
    prelim->SetTextAlign(32);
    prelim->Draw();
  }
  else{
    TLatex* prelim = new TLatex(.241, .95, "Neutrino Beam");
    prelim->SetTextColor(kBlack);
    prelim->SetNDC();
    prelim->SetTextSize(1/30.);
    prelim->SetTextAlign(32);
    prelim->Draw();
  }
}

TCanvas * ResolutionPanels(TString cname, const std::vector<TString> labels, std::string epoch, bool, double, bool, bool);
TGraph * GetBiasErrorLine (TMultiGraph* multig, TString opt, int color);
void DrawResBin(PredictionExtrap * pred, 
                osc::IOscCalc* calc, double fdpot,
                double center,
                TString sigBkg,
                int color,
                bool normalized , bool , TMultiGraph* = NULL, bool draw = true, bool tags = true);

    void ResolutionBars(std::vector<TH1D*> allResol,const std::vector<int> idc,
                        const std::vector<TString> idl, bool resid=false, bool isRHC=false);
void DifferenceOverlay(std::vector<TH1D*> allRes, const std::vector<int> idc,
                       const std::vector<TString> idl, double maxy, bool isResolution=false, bool isRHC=false);

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

void resolution2018(bool createFile = false, bool isRHC = false, std::string epoch="period3",
                   TString suffix = "test",
                   bool unoscillated = true, bool useTrue = true)
{

  TString filename = "pred_resol_2018_"+suffix+"_"+epoch+".root";

  const Var resVar = (kCCE - kTrueE)/kTrueE;
  const HistAxis resAxis ("Energy Resolution (Reco - True) / True", 
                          Binning::Simple(200,-2.5,2.5),resVar);

  const Var residVar = (kCCE - kTrueE);
  const HistAxis residAxis ("Energy Difference (Reco - True)", 
                            Binning::Simple(200,-2.5,2.5),residVar);

  const HistAxis trueAxis("True Neutrino Energy (GeV)",kNumuEnergyBinning ,kTrueE);

//  std::vector <Cut> binCuts;
//  std::vector <Cut> binCuts_truth;
//
//  for(int i=0;i<20;i++){
//    double low = i*0.25;
//    double hig = (i+1.)*0.25;
//    Cut tempCut = ((kCCE >= low) && (kCCE <hig));
//    binCuts.emplace_back(tempCut);
//    Cut tempCut_truth = ((kTrueE >= low) && (kTrueE <hig));
//    binCuts_truth.emplace_back(tempCut_truth);
//  }

  std::vector <Cut> intModes = {kIsQE, kIsRes, kIsDIS, kIsCoh, kIsDytmanMEC,kNoCut};
  std::vector <TString> intNames = {"QE", "Res", "DIS", "Coh", "DytmanMEC","allModes"};

  // --- Quantiles
  std::string fdspecfile = "/pnfs/nova/persistent/analysis/numu/Ana2018/provisional/quantiles/quantiles__rhc_full__numu2018.root";
  TFile* inFile = TFile::Open( pnfs2xrootd(fdspecfile.c_str()).c_str() );
  TH2 *FDSpec2D = (TH2*)inFile->FindObjectAny( "FDSpec2D" ); //TString("FD_full_") + TString(sFHC_RHC) );
  const int NHadEFracQuantiles = 4;
  std::vector<Cut> HadEFracQuantCuts = QuantileCutsFromTH2(FDSpec2D, kNumuCCOptimisedAxis, kHadEFracAxis, NHadEFracQuantiles);
  HadEFracQuantCuts.emplace_back( kNoCut );
  std::vector<TString> QuantNames = { "Quant1", "Quant2", "Quant3", "Quant4", "AllQuant" };


  struct ShiDef {
    TString name;
    SystShifts shift;
  };
  std::vector <ShiDef> gshifts;

  for(auto idx:{+1,-1}){
    TString sigs = (idx>0? "plus":"minus");
    // for now, just adding in the neutron systematic

    gshifts.push_back({"NeutronVisEScalePrimariesSyst2018_"+sigs, 
            SystShifts(&kNeutronVisEScalePrimariesSyst2018, idx)});


//    gshifts.emplace_back({"MaCCQEshape_"+sigs,
//        SystShifts( new GenieRwgtTableSyst(rwgt::fReweightMaCCQEshape), idx)});
//    gshifts.emplace_back({"NormCCQE_"+sigs,
//        SystShifts((ISyst*) new GenieRwgtTableSyst(rwgt::fReweightNormCCQE), idx)});
//    gshifts.emplace_back({"MvNCRES_"+sigs,
//        SystShifts((ISyst*) new GenieRwgtTableSyst(rwgt::fReweightMvNCRES), idx)});
//    gshifts.emplace_back({"MaCCRES_"+sigs,
//        SystShifts((ISyst*)new GenieRwgtTableSyst(rwgt::fReweightMaCCRES), idx)});
//    gshifts.emplace_back({"MaNCRES_"+sigs,
//        SystShifts(new GenieRwgtTableSyst(rwgt::fReweightMaNCRES), idx)});
//    gshifts.emplace_back({"MvCCRES_"+sigs,
//        SystShifts(new GenieRwgtTableSyst(rwgt::fReweightMvCCRES), idx)});  
//    gshifts.emplace_back({"CCQEPauliSupViaKF_"+sigs,
//        SystShifts(new GenieRwgtTableSyst(rwgt::fReweightCCQEPauliSupViaKF), idx)});  
//    
//    gshifts.emplace_back({"RPACCQEEnh_"+sigs,
//          SystShifts(kRPACCQEEnhSyst2018, idx)});
//    gshifts.emplace_back({"kRPACCQESupp_"+sigs,
//            SystShifts(kRPACCQESuppSyst2018, idx)});
//
//    gshifts.emplace_back({"MECq0Shape_"+sigs,
//            SystShifts(kMECq0ShapeSyst2018, idx)});
//    gshifts.emplace_back({"MECInitStateNPFrac_"+sigs,
//            SystShifts(kMECInitStateNPFracSyst2018, idx)});
//    gshifts.emplace_back({"MECEnuShape_"+sigs,
//            SystShifts(kMECEnuShapeSyst2018, idx)});
//
  }
  

  
  if(createFile){  
    ECAFType ltype = kNumuConcat; 
    Loaders::FluxType ftype = Loaders::kFHC;
    if(isRHC) ftype = Loaders::kRHC;
    Prod4NomLoaders loaders(ltype, ftype, epoch, epoch);

    std::vector <PredictionNoExtrap *> predsNXP;
    std::vector <CheatDecomp *> decomps;
    
    std::vector <TString> names;

    if(!suffix.Contains("ndonly")){    
      for(auto axis:{kNumuCCOptimisedAxis,resAxis,residAxis,trueAxis}){
        for(auto cut:intModes){
          for(auto quant:HadEFracQuantCuts){
            predsNXP.emplace_back(new PredictionNoExtrap (loaders, axis, kNumuCutFD2018 && cut && quant, 
                                                     kNoShift, 
                                                     kPPFXFluxCVWgt*kXSecCVWgt2018));
          }
        }//end modes
        for(auto shift:gshifts){
          for(auto quant:HadEFracQuantCuts){
            predsNXP.emplace_back(new PredictionNoExtrap (loaders, axis, kNumuCutFD2018 && quant, 
                                                     shift.shift, 
                                                     kPPFXFluxCVWgt*kXSecCVWgt2018));
          }
        }//end shifts
      }
    }
    /*
    if(suffix.Contains("ndonly")){
      for(auto axis:{kNumuCCOptimisedAxis,resAxis,residAxis,trueAxis}){
        decomps.emplace_back(new CheatDecomp(loaders.GetLoader(caf::kNEARDET,Loaders::kMC),
                                          axis,kNumuCutND2018,kNoShift,kPPFXFluxCVWgt*kXSecCVWgt2018
                            ));
      }
      }*/
    for(TString name:{"numu_recoE_all","numu_resol_all",
          "numu_resid_all","numu_trueE_all"}){
      for(auto intN:intNames){
        for(auto quantN:QuantNames){
          names.emplace_back(name+"_" + intN + "_" + quantN);
        }
      }
      for(auto shift:gshifts){ 
        for(auto quantN:QuantNames){
          names.emplace_back(name+"_" + shift.name + "_" + quantN);
        }
      }
    }
    
//    for (int i=0;i<20;i++){ 
//      if(!suffix.Contains("ndonly")){
//      for(auto axis:{kNumuCCOptimisedAxis,resAxis,residAxis}){
//        for(auto cut:intModes){
//            for(auto quant:HadEFracQuantCuts){
//            predsNXP.emplace_back(new PredictionNoExtrap(loaders, axis, 
//                                                    kNumuCutFD2018 && binCuts[i] && cut && quant,
//                                                    kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018));
//          }
//          }
//        for(auto shift:gshifts){
//            for(auto quant:HadEFracQuantCuts){
//            predsNXP.emplace_back(new PredictionNoExtrap(loaders, axis, 
//                                                    kNumuCutFD2018 && binCuts[i] && quant,
//                                                    shift.shift, kPPFXFluxCVWgt*kXSecCVWgt2018));
//            }
//        }
//      }
//      for(auto axis:{trueAxis,resAxis,residAxis}){
//        for(auto cut:intModes){
//          for(auto quant:HadEFracQuantCuts){
//              predsNXP.emplace_back(new PredictionNoExtrap(loaders, axis, 
//                                                    kNumuCutFD2018 && binCuts_truth[i] && cut,
//                                                    kNoShift, kPPFXFluxCVWgt*kXSecCVWgt2018));
//            }
//        }
//        for(auto shift:gshifts){
//            for(auto quant:HadEFracQuantCuts){
//               predsNXP.emplace_back(new PredictionNoExtrap(loaders, axis, 
//                                                    kNumuCutFD2018 && binCuts_truth[i],
//                                                    shift.shift, kPPFXFluxCVWgt*kXSecCVWgt2018));
//            }
//        }
//      }
//      }
//      
//      /* if(suffix.Contains("ndonly")){
//       for(auto axis:{kNumuCCOptimisedAxis,resAxis,residAxis}){
//        decomps.emplace_back(new CheatDecomp(loaders.GetLoader(caf::kNEARDET,Loaders::kMC),
//        axis,kNumuCutND2018 && binCuts[i],kNoShift,kPPFXFluxCVWgt*kXSecCVWgt2018
//        ));
//      }
//      for(auto axis:{trueAxis,resAxis,residAxis}){
//      decomps.emplace_back(new CheatDecomp(loaders.GetLoader(caf::kNEARDET,Loaders::kMC),
//      axis,kNumuCutND2018 && binCuts_truth[i],kNoShift,kPPFXFluxCVWgt*kXSecCVWgt2018
//      ));
//      }      
//      }*/
//
//      for(TString name:{"numu_recoE_bin_","numu_resol_bin_","numu_resid_bin_"}){
//      for(auto intN:intNames){
//          for(auto quantN:QuantNames){
//          names.emplace_back(name+std::to_string(i)+"_" + intN + "_" + quantN);
//          }
//      }
//      for(auto shift:gshifts){
//          for(auto quantN:QuantNames){
//          names.emplace_back(name+std::to_string(i)+"_" + shift.name + "_" + quantN);
//          }
//      }
//      }
//      
//      for(TString name:{"numu_trueE_tbin_","numu_resol_tbin_","numu_resid_tbin_"}){
//      for(auto intN:intNames){
//          for(auto quantN:QuantNames){
//          names.emplace_back(name+std::to_string(i)+"_" + intN + "_" + quantN);
//          }
//      }
//      for(auto shift:gshifts){
//          for(auto quantN:QuantNames){
//          names.emplace_back(name+std::to_string(i)+"_" + shift.name + "_" + quantN);
//          }
//      }
//      }
//    }//end bins
    
    loaders.Go();
    
    TFile * file = new TFile(filename,"recreate");
    
    
    for(int i=0;i<(int)names.size();++i){
      if(!suffix.Contains("ndonly")){
        predsNXP[i]->SaveTo(file, "pred_nxp_"+names[i]);
      }
//      if(suffix.Contains("extrap") || suffix.Contains("ndonly")){
//      decomps[i]->SaveTo(file, "decomp_"+names[i]);
      //    }
    }
    return;
  }
  
  
  else{
    //make some plots
    TFile * file = new TFile(filename,"read");
    double fdpot = kAna2018FHCPOT;
    if(isRHC) double fdpot = kAna2018RHCPOT;

    osc::NoOscillations nos;
    auto   calc = nos.Copy();
    if(!unoscillated) calc = (osc::IOscCalc*) DefaultOscCalc();
    TString oscStr ="_";
    if(unoscillated) oscStr = "_noosc_";
    if(useTrue) oscStr="_true" + oscStr;
    
    TString pstr = "pred_nxp_numu_";

    for(int qIdx = 0; qIdx < QuantNames.size(); ++qIdx){
      auto pred = PredictionExtrap::LoadFrom(file,
                                             pstr+(useTrue?"trueE":"recoE") + "_all_allModes_"+QuantNames[qIdx]);
      auto hp = pred->Predict(calc).ToTH1(fdpot);

      bool normalized = true;
      std::vector <TString> ids = {"sig","bkg","tot"};
      std::vector <int> idc = {kRed,kBlue,kBlack};
      std::vector <TString> idl =  {"#nu_{#mu} CC","Bkgd.","Total"};

      //all bins
      auto c1_h = new TCanvas();
      TString ps_all = pstr + "resid_all_allModes_" + QuantNames[qIdx];
      auto resid =  PredictionExtrap::LoadFrom(file, ps_all);
      std::vector <TH1D*> allRes;
      allRes.emplace_back((resid->ComponentCC(14,14).Oscillated(calc,14,14) +
                          resid->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));     
      allRes.emplace_back((resid->Predict(calc) -  
                          resid->ComponentCC(14,14).Oscillated(calc,14,14) - resid->ComponentCC(-14,-14).Oscillated(calc,-14,-14))
                         .ToTH1(fdpot)); 
      allRes.emplace_back(resid->Predict(calc).ToTH1(fdpot));     
      
      DifferenceOverlay(allRes,idc,idl, 0.1, false, isRHC);
      
      auto c1_hu = new TCanvas();
      std::vector <TH1D*> allResu = {(resid->ComponentCC(14,14).Oscillated(calc,14,14)+resid->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot)};
      DifferenceOverlay(allResu,{idc[0]},{idl[0]},0.15, false, isRHC);
      
      
      ps_all = pstr + "resol_all_allModes_"+QuantNames[qIdx];
      auto resol =  PredictionExtrap::LoadFrom(file, ps_all);
      std::vector <TH1D*> allResol;
      allResol.emplace_back((resol->ComponentCC(14,14).Oscillated(calc,14,14) +
                          resol->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));     
      allResol.emplace_back((resol->Predict(calc) -  
                          resol->ComponentCC(14,14).Oscillated(calc,14,14) - resol->ComponentCC(-14,-14).Oscillated(calc,-14,-14))
                         .ToTH1(fdpot));     
      allResol.emplace_back(resol->Predict(calc).ToTH1(fdpot));     

      auto c1_h2 = new TCanvas;
      DifferenceOverlay(allResol,idc,idl, 0.2, true, isRHC);

      auto c1_h2u = new TCanvas();
      std::vector <TH1D*> allResolu = {(resol->ComponentCC(14,14).Oscillated(calc,14,14) +
                          resol->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot)};
      DifferenceOverlay(allResolu,{idc[0]},{idl[0]},0.1);

      auto c1_mg = new TCanvas();
      ResolutionBars(allResol, idc,idl, false, isRHC);

      auto c1_mg2 = new TCanvas();
      ResolutionBars(allRes, idc,idl, true, isRHC);


      //all bins
      //repeat for numucc mec
      //WTF Erika
  //    std::vector <TString> oldN = {"DIS","Res","QE","DytmanMEC"};
      std::vector <TString> oldN = {"All","DytmanMEC","QE","Res","DIS"};
      std::vector <TString> interN = {"All", "MEC","QE","RES","DIS"};
      std::vector <int> interC = {kBlack, kOrange-1,kAzure+7, kGreen+2,kGray+2};
      
      std::vector <TString> cuteLabel;
      for (int i =0; i<5; ++i){
        TString tstr; tstr.Form("#color[%d]{%s}",interC[i],interN[i].Data());
        cuteLabel.emplace_back(tstr);
      }

      std::vector <TH1D*> allResM;
      std::vector <TH1D*> allResolM;

      allResM.emplace_back((resid->ComponentCC(14,14).Oscillated(calc,14,14)+resid->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
      allResolM.emplace_back((resol->ComponentCC(14,14).Oscillated(calc,14,14)+resol->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
      for(int i=1;i<5;++i){
        ps_all = pstr + "resid_all_"+oldN[i]+"_"+QuantNames[qIdx];
        auto residM =  PredictionExtrap::LoadFrom(file, ps_all);
        allResM.emplace_back((residM->ComponentCC(14,14).Oscillated(calc,14,14)+residM->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
        
        ps_all = pstr + "resol_all_"+oldN[i]+"_"+QuantNames[qIdx];
        auto resolM =  PredictionExtrap::LoadFrom(file, ps_all);
        allResolM.emplace_back((resolM->ComponentCC(14,14).Oscillated(calc,14,14)+resolM->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
      }

      auto c2_h = new TCanvas();
      DifferenceOverlay(allResM,interC,{"All","MEC","QE","RES","DIS"},0.15, false, isRHC);

      auto c2_h2 = new TCanvas;
      DifferenceOverlay(allResolM,interC,interN, 0.25, true, isRHC);

      auto c2_mg = new TCanvas();
      ResolutionBars(allResolM, interC,{"All","MEC","QE","RES","DIS"}, false, isRHC);
      auto c2_mg2 = new TCanvas();
      ResolutionBars(allResM, interC,{"All", "MEC","QE","RES","DIS"}, true, isRHC);




  //////////////////////////////////////////////////////////////////////////////////
  //repeat for genie
  //all bins
      //repeat for numucc mec
      std::vector<TString> gOld = {"NeutronVisEScalePrimariesSyst2018_"};

      std::vector<TString> gShort = {"NeutronSyst"};

  //    std::vector <TString> gOld = {"NormCCQE_",
  //                              "MaCCQEshape_",
  //                              "CCQEPauliSupViaKF_",
  //                              "RPACCQE_",
  //                              "MECScale_",
  //                              "MaCCRES_",
  //                              "MaNCRES_",
  //                              "MvCCRES_",
  //                              "MvNCRES_"  };
  //
  //    std::vector <TString> gShort= {"Norm^{CCQE}",  
  //                               "M_{A}^{CCQE,sc}",
  //                               "Fermi^{CCQE}",
  //                               "RPA^{CCQE}",
  //                               "MEC^{sc}",
  //                               "M_{A}^{CCRES}",
  //                               "M_{A}^{NCRES}",
  //                               "M_{v}^{CCRES}",
  //                               "M_{v}^{NCRES}"
  //    };
  //
  //      
  //      
        for (auto part:{1,1}){
          int ming = 0; int maxg = 1;
  //    if(part==2) {ming=5; maxg=9;}

          std::vector <TString> gOldSig;
          std::vector <TString> gShortSig;
          std::vector <int> gColors;
        
          for(int i=ming;i<maxg;i++){
            gOldSig.emplace_back(gOld[i]+"plus");
            gShortSig.emplace_back(gShort[i]+" (#plus1#sigma)");
            gColors.emplace_back(kPink -i);

            gOldSig.emplace_back("Nominal");
            gShortSig.emplace_back("Nominal");
            gColors.emplace_back(kBlack);
            
            gOldSig.emplace_back(gOld[i]+"minus");
            gShortSig.emplace_back(gShort[i]+" (#minus1#sigma)");
            gColors.emplace_back(kAzure -i);
          }
          
          std::vector <TH1D*> allResG;
          std::vector <TH1D*> allResolG;

          // this is a hack for the neutron systematic only
//          for(int i=0;i<(int)gOldSig.size();i++){ 
          for(int i=0;i<1;i++){ 

            ps_all = pstr + "resid_all_"+gOldSig[i]+"_"+QuantNames[qIdx];
            std::cout <<ps_all << std::endl;
            auto residM =  PredictionExtrap::LoadFrom(file, ps_all);
            allResG.emplace_back((residM->ComponentCC(14,14).Oscillated(calc,14,14)+residM->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
            ps_all = pstr + "resol_all_"+gOldSig[i]+"_"+QuantNames[qIdx];
            auto resolM =  PredictionExtrap::LoadFrom(file, ps_all);
            allResolG.emplace_back((resolM->ComponentCC(14,14).Oscillated(calc,14,14)+resolM->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
          }

            // Nominal
            allResG.emplace_back((resid->ComponentCC(14,14).Oscillated(calc,14,14) +
                          resid->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));     
            allResolG.emplace_back((resol->ComponentCC(14,14).Oscillated(calc,14,14) +
                          resol->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));     

          for(int i=2;i<3;i++){ 

            ps_all = pstr + "resid_all_"+gOldSig[i]+"_"+QuantNames[qIdx];
            std::cout <<ps_all << std::endl;
            auto residM =  PredictionExtrap::LoadFrom(file, ps_all);
            allResG.emplace_back((residM->ComponentCC(14,14).Oscillated(calc,14,14)+residM->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
            ps_all = pstr + "resol_all_"+gOldSig[i]+"_"+QuantNames[qIdx];
            auto resolM =  PredictionExtrap::LoadFrom(file, ps_all);
            allResolG.emplace_back((resolM->ComponentCC(14,14).Oscillated(calc,14,14)+resolM->ComponentCC(-14,-14).Oscillated(calc,-14,-14)).ToTH1(fdpot));
          }


//          }
          
          auto c3_h = new TCanvas();
          DifferenceOverlay(allResG,gColors, gShortSig,0.1, false, isRHC);

          auto c3_h2 = new TCanvas();
          DifferenceOverlay(allResolG,gColors, gShortSig,0.2,true, isRHC);

          auto c3_mg = new TCanvas();
          ResolutionBars(allResolG, gColors, gShortSig, false, isRHC);
          auto c3_mg2 = new TCanvas();
          ResolutionBars(allResG, gColors, gShortSig, true, isRHC);
          
          c3_mg->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bars_xsec"+std::to_string(part)+"_"+QuantNames[qIdx]+".pdf");
          c3_mg2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bars2_xsec"+std::to_string(part)+"_"+QuantNames[qIdx]+".pdf");
          c3_h->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlay_xsec"+std::to_string(part)+"_"+QuantNames[qIdx]+".pdf");
          c3_h2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlayf_xsec"+std::to_string(part)+"_"+QuantNames[qIdx]+".pdf");
        }
  //      //      return;
  //      
  //      ////Need to do those before gstyle change im an idiot
  //
  //
  //      for (auto part:{1,1}){
  //    int ming = 0; int maxg = 1;
  //    if(part==2) {ming=5; maxg=9;}
  //
  //    std::vector <TString> gOldSig;
  //    std::vector <TString> gShortSig;
  //    std::vector <int> gColors;
  //      
  //    for(int i=ming;i<maxg;i++){
  //      gOldSig.emplace_back(gOld[i]+"plus");
  //      gShortSig.emplace_back(gShort[i]+" (#plus1#sigma)");
  //      gColors.emplace_back(kPink -i);
  //      
  //      gOldSig.emplace_back(gOld[i]+"minus");
  //      gShortSig.emplace_back(gShort[i]+" (#minus1#sigma)");
  //      gColors.emplace_back(kAzure -i);
  //    }
  //    
  //    auto c3_f =  ResolutionPanels("genie", gShortSig, epoch, useTrue, 
  //                                  1.2*hp->GetMaximum(),true);
  //    
  //
  //    TLegend * legG = new TLegend(0.08,0.58,0.95,0.84);
  //    legG->SetFillStyle(0);
  //    legG->SetNColumns(4);
  //    legG->SetTextSize(0.06);
  //    
  //    for(int idx = 0; idx <gShortSig.size(); ++idx)
  //    {
  //      TMultiGraph* gr_resol = new TMultiGraph();
  //      
  //      for(int i =0; i<20;i++){
  //        TString ps_bin = (TString)pstr+(useTrue?+"resol_tbin_":"resol_bin_")
  //          + std::to_string(i)+"_";
  //        auto mbin = PredictionExtrap::LoadFrom(file, ps_bin+gOldSig[idx]);
  //        
  //        DrawResBin(mbin.release(),calc,fdpot,(double)i*0.25+.125,"sig",kBlack,
  //                   normalized, unoscillated,gr_resol,false, false);
  //      }
  //      
  //      //Transfotm pretty multigraph into boring lines
  //      c3_f->cd(2);
  //      auto bias = GetBiasErrorLine(gr_resol,"bias",gColors[idx]);
  //      bias->Draw("Lp same");
  //      legG->AddEntry(bias,gShortSig[idx],"lp");
  //      gPad->RedrawAxis();
  //      c3_f->cd(1);
  //      GetBiasErrorLine(gr_resol,"error",gColors[idx])->Draw("lp same");
  //      gPad->RedrawAxis();
  //    }
  //
  //    c3_f->cd(1);
  //    legG->Draw();
  //    
  //    c3_f->cd(2);
  //    legG->Draw();
  //    
  //
  //    c3_f->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"fracs_xsec"+
  //                std::to_string(part)+"_"+QuantNames[qIdx]+".pdf");
  //    
  //      }//end parts genie


        
  //      auto c1 = ResolutionPanels("total", {"#color[2]{#nu_{#mu}CC}","#color[4]{Bkgd.}","#color[1]{Total}"}, epoch,useTrue, 1.2*hp->GetMaximum(),false);
  //      auto c1_f =  ResolutionPanels("total", {"True #nu_{#mu} CC signal","Background","Total (Sig. + Bkg.)"}, epoch, useTrue, 1.2*hp->GetMaximum(),true);
  //      auto c1_fu =  ResolutionPanels("numu", {"True #nu_{#mu} CC signal"},epoch, useTrue, 1.2*hp->GetMaximum(),true);
  //
  //
  //    TLegend * leg3 = new TLegend(0.08,0.74,0.68,0.84);
  //    leg3->SetFillStyle(0);
  //    leg3->SetNColumns(3);
  //
  //    int mycolors[]={kPink-5,kSpring-7,kOrange+7,kAzure+2,kPink+10,/*kTeal-5,*/kViolet-6,kCyan+3};
  //    auto bkgcolor = kGray+3;
  //
  //    //Draw the spectrum
  //    //Pretty bin colors
  //    c1->cd(1);
  //    for(int i =0; i<20;i++){ 
  //      auto pbin =  PredictionExtrap::LoadFrom(
  //    file->GetDirectory((TString)pstr+(useTrue?+"trueE_tbin_":"recoE_bin_")+ std::to_string(i)+"_allModes"));  
  //      int color = mycolors[i%7]; 
  //      auto htot_bin = pbin->Predict(calc).ToTH1(fdpot,color);
  //      auto hbkg_bin = (pbin->Predict(calc) - pbin->ComponentCC(14,14).Oscillated(calc,14,14)).ToTH1(fdpot,bkgcolor);
  //      htot_bin->SetFillColorAlpha(color,.1);
  //      hbkg_bin->SetFillColor(bkgcolor);
  //      htot_bin->Draw("hist same");
  //      hbkg_bin->Draw("hist f same");
  //      double offset = 0.1 * hp->GetMaximum();
  //      TLatex * ltx = new TLatex();
  //      ltx->SetTextAlign(22);
  //      ltx->SetTextSize(1.2*ltx->GetTextSize());
  //      TString binStr; binStr.Form("#color[%d]{#bf{B%d}}",color,i+1);
  //      ltx->DrawLatex((double)i*.25+0.125, htot_bin->Integral()+offset,binStr);
  //    }
  //    
  //    for(int idx = 0; idx <3; ++idx)
  //    {
  //      TMultiGraph* gr_resid = new TMultiGraph();
  //      TMultiGraph* gr_resol = new TMultiGraph();
  //
  //      for(int i =0; i<20;i++){
  //    int color = mycolors[i%7]; 
  //    //True - reco graphs
  //    auto pbin =  PredictionExtrap::LoadFrom(
  //      file->GetDirectory((TString)pstr+(useTrue?+"resid_tbin_":"resid_bin_")+ std::to_string(i)+"_allModes")).release();  
  //    c1->cd(idx+2);
  //    DrawResBin(pbin,calc,fdpot,(double)i*0.25+.125,ids[idx],color,normalized,unoscillated,gr_resid, true, true);
  //        
  //        //Get the resol info; false since we are not drawing these anymore   
  //    pbin =  PredictionExtrap::LoadFrom(
  //      file->GetDirectory((TString)pstr+(useTrue?+"resol_tbin_":"resol_bin_")+ std::to_string(i)+"_allModes")).release();  
  //    DrawResBin(pbin,calc,fdpot,(double)i*0.25+.125,ids[idx],color,normalized,unoscillated,gr_resol,false, false);
  //      gPad->RedrawAxis();
  //      }
  //      
  //      //Transfotm pretty multigraph into boring lines
  //      c1_f->cd(2);
  //      auto bias = GetBiasErrorLine(gr_resol,"bias",idc[idx]);
  //      bias->Draw("Lp same");
  //      leg3->AddEntry(bias,idl[idx],"lp");
  //      gPad->RedrawAxis();
  //      c1_f->cd(1);
  //      GetBiasErrorLine(gr_resol,"error",idc[idx])->Draw("lp same");
  //      gPad->RedrawAxis();
  //
  //      if(idx==0){
  //    c1_fu->cd(1);
  //    GetBiasErrorLine(gr_resol,"error",idc[idx])->Draw("lp same");
  //    leg3->Draw();
  //    gPad->RedrawAxis();
  //    c1_fu->cd(2);
  //    GetBiasErrorLine(gr_resol,"bias",idc[idx])->Draw("lp same");
  //    leg3->Draw();
  //    gPad->RedrawAxis();
  //    c1_fu->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"fracs_numu"+"_"+QuantNames[qIdx]+.pdf");
  //      }
  //    }
  //
  //    c1_f->cd(1);
  //    leg3->Draw();
  //
  //    c1_f->cd(2);
  //    leg3->Draw();
  //
  //    c1->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bins"+"_"+QuantNames[qIdx]+".pdf");
  //    c1_f->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"fracs"+"_"+QuantNames[qIdx]+".pdf");
      c1_mg->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bars"+"_"+QuantNames[qIdx]+".pdf");
      c1_mg2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bars2"+"_"+QuantNames[qIdx]+".pdf");
      c1_h->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlay"+"_"+QuantNames[qIdx]+".pdf");
      c1_h2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlayf"+"_"+QuantNames[qIdx]+".pdf");
      c1_hu->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlay_numu"+"_"+QuantNames[qIdx]+".pdf");
      c1_h2u->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlayf_numu"+"_"+QuantNames[qIdx]+".pdf");

    
    ////////////////////////////////////////////////////////////////////////////////
      
  //    auto c2 = ResolutionPanels("mode", cuteLabel,epoch,useTrue, 1.2*hp->GetMaximum(),false);
  //    auto c2_f = ResolutionPanels("mode", interN,epoch,useTrue, 1.2*hp->GetMaximum(),true);

      
  //    c2->cd(1);
  //    for(int i =0; i<20;i++){ 
  //      TString ps_bin = (TString)pstr+(useTrue?+"trueE_tbin_":"recoE_bin_")
  //    + std::to_string(i)+"_";
  //      auto pbin =  PredictionExtrap::LoadFrom(file, ps_bin+"allModes");  
  //      auto cbin =  PredictionExtrap::LoadFrom(file, ps_bin+"Coh");  
  //
  //      THStack *hs = new THStack("hs","");      
  //      //Other is bkg + coh
  //      int color = mycolors[i%7]; 
  //      auto htot_bin = pbin->Predict(calc).ToTH1(fdpot,/*color*/kBlack);
  //      auto hoth_bin = (pbin->Predict(calc) - 
  //                   pbin->ComponentCC(14,14).Oscillated(calc,14,14) + 
  //                   cbin->ComponentCC(14,14).Oscillated(calc,14,14))
  //                                          .ToTH1(fdpot,bkgcolor);
  //      hoth_bin->SetFillColorAlpha(bkgcolor, 1);
  //      hoth_bin->SetLineColorAlpha(bkgcolor, 0);
  //      hs->Add(hoth_bin,"hist f");
  //
  //      for (int idx=3;idx>=0;--idx){
  //    auto mbin =  PredictionExtrap::LoadFrom(file, ps_bin+oldN[idx]);
  //    auto hmod_bin = mbin->ComponentCC(14,14).Oscillated(calc,14,14).ToTH1(fdpot,interC[idx]);
  //    hmod_bin->SetFillColorAlpha(interC[idx], 0.3);
  //    hmod_bin->SetLineColorAlpha(interC[idx], 0);
  //    hs->Add(hmod_bin,"hist f");
  //      }
  //      hs->Draw("f same");
  //      htot_bin->Draw("hist same");
  //  
  //      double offset = 0.1 * hp->GetMaximum();
  //      TLatex * ltx = new TLatex();
  //      ltx->SetTextAlign(22);
  //      ltx->SetTextSize(1.2*ltx->GetTextSize());
  //      TString binStr; binStr.Form("#color[%d]{#bf{B%d}}",color,i+1);
  //      ltx->DrawLatex((double)i*.25+0.125, htot_bin->Integral()+offset,binStr);
  //    }//end bins
  //
  //    TLegend * legM = new TLegend(0.08,0.74,0.68,0.84);
  //    legM->SetFillStyle(0);
  //    legM->SetNColumns(4);
  //
  //    for(int idx = 0; idx <4; ++idx)
  //    {
  //      TMultiGraph* gr_resid = new TMultiGraph();
  //      TMultiGraph* gr_resol = new TMultiGraph();
  //      
  //      for(int i =0; i<20;i++){
  //    //True - reco graphs
  //    int color = mycolors[i%7];
  //    TString ps_bin = (TString)pstr+(useTrue?+"resid_tbin_":"resid_bin_")
  //      + std::to_string(i)+"_";
  //    auto mbin = PredictionExtrap::LoadFrom(file, ps_bin+oldN[idx]);
  //    c2->cd(idx+2);
  //    DrawResBin(mbin.release(),calc,fdpot,(double)i*0.25+.125,"sig",color,
  //               normalized,unoscillated,gr_resid, true, true);
  //    gPad->RedrawAxis();
  //        
  //    ps_bin = (TString)pstr+(useTrue?+"resol_tbin_":"resol_bin_")
  //      + std::to_string(i)+"_";
  //    mbin = PredictionExtrap::LoadFrom(file, ps_bin+oldN[idx]);
  //    
  //    DrawResBin(mbin.release(),calc,fdpot,(double)i*0.25+.125,"sig",color,
  //               normalized, unoscillated,gr_resol,false, false);
  //      }
  //      
  //      //Transfotm pretty multigraph into boring lines
  //      c2_f->cd(2);
  //      auto bias = GetBiasErrorLine(gr_resol,"bias",interC[idx]);
  //      bias->Draw("Lp same");
  //      legM->AddEntry(bias,interN[idx],"lp");
  //      gPad->RedrawAxis();
  //      c2_f->cd(1);
  //      GetBiasErrorLine(gr_resol,"error",interC[idx])->Draw("lp same");
  //      gPad->RedrawAxis();
  //    }
  //
  //    c2_f->cd(1);
  //    legM->Draw();
  //
  //    c2_f->cd(2);
  //    legM->Draw();
  //
  //
  //
  //    c2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bins_mode"+"_"+QuantNames[qIdx]+".pdf");
  //    c2_f->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"fracs_mode"+"_"+QuantNames[qIdx]+".pdf");
      c2_mg->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bars_mode"+"_"+QuantNames[qIdx]+".pdf");
      c2_mg2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"bars2_mode"+"_"+QuantNames[qIdx]+".pdf");
      c2_h->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlay_mode"+"_"+QuantNames[qIdx]+".pdf");
      c2_h2->Print("plot4_resol_2018_"+suffix+"_"+epoch + oscStr +"overlayf_mode"+"_"+QuantNames[qIdx]+".pdf");
    } // end quant
  }//end plots  
}

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

TGraph * GetBiasErrorLine (TMultiGraph* multig, TString opt, int color){
  TGraph * ret = new TGraph();
  auto ngraphs = multig->GetListOfGraphs()->GetSize();


  for (int idx=0; idx<ngraphs; ++idx){
    auto gr = (TGraphErrors*) multig->GetListOfGraphs()->At(idx); 
    auto x = gr->GetX()[0];
    auto y = gr->GetY()[0];
    auto erry = gr->GetEY()[0];
    if(opt=="bias") ret->SetPoint(ret->GetN(),x,y);
    if(opt=="error") ret->SetPoint(ret->GetN(),x,erry);
  }
  ret->SetLineWidth(3);
  ret->SetLineColor(color);
  ret->SetMarkerStyle(kFullCircle);
  ret->SetMarkerColor(color);
  
  return ret;

}


void FormatResHist(TH1D* hist, TString xaxis, TString yaxis,
                   bool choplabel,
                   double xmin, double xmax){
  hist->SetNdivisions(124);
  hist->GetXaxis()->SetTitle(xaxis);
  hist->GetYaxis()->SetTitle(yaxis);
  hist->GetYaxis()->CenterTitle();
  hist->GetYaxis()->SetDecimals();  

  if(choplabel) {
    hist->GetXaxis()->SetLabelSize(0);
    hist->GetYaxis()->SetLabelSize(0);
    hist->GetXaxis()->SetTitleOffset(1.9);
  }
  else {
    hist->GetXaxis()->SetLabelSize(15);
    hist->GetYaxis()->SetLabelSize(15);
    hist->GetXaxis()->SetTitleOffset(3);
  }
  hist->GetYaxis()->SetRangeUser(xmin,xmax);
    
}



TCanvas * ResolutionPanels(TString cname, const std::vector<TString> labels, std::string epoch, bool useTrue, double max1, bool fracOnly, bool isRHC){
  
  gStyle->SetTitleSize(17, "xyz");
  gStyle->SetTitleFont(43, "xyz");
  gStyle->SetLabelSize(15, "xyz");
  gStyle->SetLabelFont(43, "xyz");
  gStyle->SetTitleOffset(0.8, "x");
  gStyle->SetTitleOffset(1.4, "yz");

  TCanvas * c ;
  TH1D* hist = new TH1D (UniqueName().c_str(),"",230,-0.25,5.5);
  TString xaxis,yaxis;
  xaxis = "Reconstructed neutrino energy (GeV)";
  if(useTrue) xaxis = "True neutrino energy (GeV)";
  int npads= labels.size();
    
  //Divide the canvas and construct the panels
  if(fracOnly){
    c = new TCanvas (cname+"f",cname,700,550);
    c->Divide(1,2);
    for(int i=1;i<=2;++i){
      c->cd(i);
      if(i==1) FormatResHist(hist, xaxis,  "Frac. E resolution", false,   0, 0.18);
      if(i==2) FormatResHist(hist, xaxis,  "Frac. shift in mean",  false, -0.18, 0.18);

      gPad->SetTopMargin(0.15);    
      gPad->SetBottomMargin(0.18);
      gPad->SetRightMargin(0.02*5./7.);
      gPad->SetLeftMargin(0.1 *5./7.);

      hist->GetXaxis()->SetTitleOffset(2);
      hist->GetYaxis()->SetTitleOffset(1.2);

      hist->DrawCopy();
      gPad->SetGridx();
      gPad->SetFillStyle(0);
      gPad->SetGridy();
    }
  }
  else{
    c = new TCanvas(cname,cname,700,700);
    double del = 0.66/npads;
    c->Divide(1,1+npads);
    c->GetPad(1)->SetPad(0,0.66,1,1);
    for(int i=0;i<npads; ++i){
      c->GetPad(2+i)->SetPad(0,0.66-(i+1)*del,1,0.66-i*del);
    }
    for(int i=1; i<= 1+npads ; ++i){
      c->cd(i);
      gPad->SetBottomMargin((i>1?0.18:0.16));
      gPad->SetRightMargin(0.02);
      gPad->SetLeftMargin(0.08);
      if(i>1) gPad->SetTopMargin(0.05);
      else gPad->SetTopMargin(0.15);
      
      if(i==1) FormatResHist(hist, xaxis,  "Events/bin", false,    0, max1);
      else FormatResHist(hist, "Energy  difference (True-Reco) (GeV)", 
                         "A.U.", true, 0, 0.77);
      hist->DrawCopy();
      gPad->SetGridx();
      gPad->SetFillStyle(0);
    }
  }
  
  
  TLatex * ltx = new TLatex();
  ltx->SetTextSize(19);
  ltx->SetTextFont(43);
  ltx->SetTextAlign(11);
  
  c->cd();
  gPad->SetFillStyle(0);

  if(!fracOnly) {
    ltx->DrawLatexNDC(0.13,0.9,"#splitline{FD MC}{6 #times 10^{20} POT}");      
    for(int i = 0; i < npads; ++i){
      c->cd(i+2);
      ltx->DrawLatexNDC(0.13,0.6, "#bf{"+labels[i]+"}");
    }
  }
  c->cd(1); Simulation(); Beam(isRHC);
  if(fracOnly)  c->cd(2); Simulation(); Beam(isRHC);
  
  return c;
  
}



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



void DrawResBin(PredictionExtrap * pred, 
                osc::IOscCalc* calc, double fdpot,
                double center,
                TString sigBkg,
                int color,
                bool normalized,
                bool unoscillated  ,
                TMultiGraph * graph, 
                bool draw, 
                bool tags
  ){
  gPad->SetFillStyle(0);

  TH1D* hist = new TH1D (UniqueName().c_str(),"temp",280,-1,6);
  TH1D* hpred = new TH1D (UniqueName().c_str(),"p",10,0,1);

  if(sigBkg=="tot") hpred = (*pred).Predict(calc).ToTH1(fdpot,color);
  else if(sigBkg=="sig") hpred= ((*pred).ComponentCC(14,14).Oscillated(calc,14,14)+(*pred).ComponentCC(-14,-14).Oscillated(calc,-14,-14)).
                           ToTH1(fdpot,color);
  else if(sigBkg=="bkg") hpred = ((*pred).Predict(calc)-(*pred).ComponentCC(14,14).Oscillated(calc,14,14)
                           -(*pred).ComponentCC(-14,-14).Oscillated(calc,-14,-14)).
                           ToTH1(fdpot,color);
  else {std::cerr <<"imconfused\n";}

  hpred->ResetStats();
  
  if(hpred->Integral()>0.01)
  {  
    if(tags && draw){
      
      TLatex * ltx = new TLatex(); 
      ltx->SetTextSize(9);
      ltx->SetTextFont(43);
      ltx->SetTextSize(1.2*ltx->GetTextSize());
      double y, lower1,lower2;
      if((int)((center-0.125)/0.25)%2==0){  
        y=0.61;
        lower1=0;
        lower2=-0.3;
      }
      else{
        y = 0.47;
        lower1=0.5;
        lower2=-0.3;
      }
      TString label; 
      label.Form
        ("#color[%d]{#bf{#splitline{#lower[%.2f]{#mu %.2f}}{#lower[%.2f]{#sigma %.2f}}}}",
         color,lower1,hpred->GetMean(),lower2,hpred->GetRMS());
      ltx->DrawLatex(center-0.125, y,label);
    }
    if(graph){
      TGraphErrors * mygraph = new TGraphErrors();
      mygraph->SetPoint(mygraph->GetN(),center, hpred->GetMean()) ;
      mygraph->SetPointError(mygraph->GetN()-1, 0, hpred->GetRMS());   
      mygraph->SetMarkerStyle(kFullCircle);
      mygraph->SetLineColor(color);
      mygraph->SetFillColor(color);
      mygraph->SetMarkerColor(color);
      mygraph->SetLineWidth(3);
      graph->Add(mygraph, "lp same");
    }
  }
  
  for(int i=1;i<=hpred->GetNbinsX();++i){
    int bin = hist->FindBin(hpred->GetBinCenter(i)+center);
    hist->SetBinContent(bin,hpred->GetBinContent(i));
  }

  hist->SetLineColor(color);
//  hist->SetFillColorAlpha(color,.10);
  hist->Rebin();
  if(draw && hpred->Integral()>0.001){
    if(normalized) hist->DrawNormalized("hist same");
    else hist->DrawCopy("hist same");
  }

  delete hist, hpred;
}



   void DifferenceOverlay(std::vector<TH1D*> allRes, const std::vector<int> idc,
                          const std::vector<TString> idl, double maxy, bool isResolution, bool isRHC){
     gPad->SetFillStyle(0);
     TH1D * haux = new TH1D(UniqueName().c_str(),
                            ";Energy difference (Reco - True) (GeV); A.U. (Area normalized)",100,
                             -1,0.5);
     if(isResolution) haux->GetXaxis()->SetTitle("(Reco - True)/True");
     haux->GetYaxis()->SetRangeUser(0,maxy);
      haux->GetYaxis()->SetLabelSize(0);
      haux->GetXaxis()->CenterTitle();
      haux->GetYaxis()->CenterTitle();
      haux->DrawClone();
      TLegend* leg = new TLegend(0.12,0.2,0.4,0.89);
      leg->SetFillStyle(0); 
      leg->SetTextSize(haux->GetXaxis()->GetTitleSize());
      //hack
//      if(allRes.size()>7) {
//      leg->SetTextAlign(31);
//      leg->SetX2(0.4);
//      leg->SetTextSize(0.9*haux->GetXaxis()->GetTitleSize());
//      }
      for (int i=0;i< (int)allRes.size();++i){
//        if(allRes[i]->GetMaximum()>maxy) maxy = allRes[i]->GetMaximum();
        allRes[i]->SetLineColor(idc[i]);
        allRes[i]->SetLineWidth(3);
//      allRes[i]->Rebin();
        allRes[i]->DrawNormalized("hist same");
        leg->AddEntry(allRes[i],idl[i],"l");
      }
       
//     haux->GetYaxis()->SetRangeUser(0,maxy);
      leg->SetHeader("FD MC");
      leg->Draw();
//      gPad->SetGridx();
      TLine *line = new TLine(0, 0, 0, maxy);
      line->SetVertical();
      line->SetLineColor(13);
      line->SetLineStyle(7);
      line->Draw("same");
      gPad->RedrawAxis();
      Simulation(); Beam(isRHC);
      delete haux;
    }


    void ResolutionBars(std::vector<TH1D*> allResol,const std::vector<int> idc,
                        const std::vector<TString> idl, bool resid, bool isRHC){
      gPad->SetFillStyle(0);
      int nbins = allResol.size();
      TH2D* haux2D = new TH2D(UniqueName().c_str(),"; (True-Reco)/True;",1000,-0.5,0.5,
                              nbins,0,nbins);
      if(resid) haux2D ->GetXaxis()->SetTitle("(True-Reco) (GeV)");
      haux2D->GetXaxis()->CenterTitle();
      auto yaxis =  haux2D->GetYaxis();
      gPad->SetLeftMargin(.25);
      for(int i=0;i<nbins;++i){
        int bin = nbins - i;
        std::cout << "Test labels" << idl[i] << std::endl;
        yaxis->SetBinLabel(bin,idl[i]);
      }
      yaxis->SetLabelSize(0.07);
      haux2D->DrawClone();
  
      for(int i=0;i<nbins;++i){ 
        int bin = nbins - i;
        auto  mygraph = new TGraphErrors();
        int color = idc[i];
        allResol[i]->ResetStats();
        mygraph->SetPoint(mygraph->GetN(),allResol[i]->GetMean(), yaxis->GetBinCenter(bin)) ;
        mygraph->SetPointError(mygraph->GetN()-1, allResol[i]->GetRMS(),0);   
        mygraph->SetMarkerStyle(kFullCircle);
        mygraph->SetLineColor(color);
        mygraph->SetFillColor(color);
        mygraph->SetMarkerColor(color);
        mygraph->SetLineWidth(4);
        mygraph->DrawClone("pl");
        
        TLatex * ltx = new TLatex();
        ltx->SetTextAlign(22);
        ltx->SetTextSize(0.04);
        if(resid){
          ltx->DrawLatex(0.3,yaxis->GetBinCenter(bin),
                         TString::Format("#color[%d]{#splitline{Mean %.1f MeV}{RMS  %.1f MeV}}",
                                         color,
                                       allResol[i]->GetMean()*1000, allResol[i]->GetRMS()*1000));

        }else{
        ltx->DrawLatex(0.3,yaxis->GetBinCenter(bin),
                       TString::Format("#color[%d]{#splitline{Mean %.1f \%}{RMS  %.1f \%}}",
                                       color,
                                       allResol[i]->GetMean()*100, allResol[i]->GetRMS()*100));

        }
      }
      gPad->SetGridx();
      Simulation(); Beam(isRHC);
      delete haux2D;
    }