doTemplateFit.C

Go to the documentation of this file.

#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/ISyst.h"
#include "CAFAna/Core/SystShifts.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Systs/XSecSystLists.h"
#include "3FlavorAna/Cuts/NueCutsSecondAna.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/XSec/GenieMultiverseSyst.h"
#include "CAFAna/XSec/FluxMultiverseSyst.h"
#include "CAFAna/Vars/XsecTunes.h"

#include "TMath.h"
#include "Math/ProbFunc.h"

#include "NDAna/nuecc_inc/NueCCIncTemplateFit.h"
#include "NDAna/nuecc_inc/NueCCIncVars.h"
#include "NDAna/nuecc_inc/NueCCIncCuts.h"
#include "NDAna/nuecc_inc/NueCCIncExtra.h"
#include "NDAna/nuecc_inc/NueCCIncCrossSectionAnalysis.h"
#include "NDAna/nuecc_inc/NueCCIncTemplateFitter.h"

using namespace ana;

const double pot = 8.09e20;

const std::string sDir       = "/nova/ana/users/mjudah/CrossSection_Feb2019/";
const std::string sTemplates = "fAnalysisSpectra.root";

const bool plotFitTemplateResults = false;
const bool plotFitAnalysisResults = false;
const bool plotSignalPrediction   = true;


struct axis_info
{
  int color;
  int marker_style;
  std::string title;
  std::string cutName;
};

const axis_info kChnsInfo[kcNumChns+1] = 
  {
    {1,     20,"Data",               "data"},
    {632+2, 1,"Total",              "mc"},
    {600+2, 1,"#nu_{e} CC",         "nuecc"},
    {616,   1,"#bar{#nu}_{e} CC",   "nuebarcc"},
    {416-3, 1,"#nu_{#mu} CC",       "numucc"},
    {416+3, 1,"#bar{#nu}_{#mu} CC", "numubarcc"},
    {800,   1,"NC",                 "nc"},
    {400-3, 1,"Other",              "other"},
    {416-3, 1,"Non-Fiducial Signal","nonfidnue"},
  };

//sort Vector
bool sort_in_ascending_order (int i,int j) { return (i<j); }


float BinSigma(std::vector<double> events, float nsigmas, float pivot)
  {
    int    pivotbin = 0;
    double pivotbincenter = 0;

    std::sort(events.begin(), events.end());

    for(unsigned int i = 0; i < events.size()-1; i++)
      if(pivot >= events[i] && pivot < events[i+1]){
        pivotbin = i;
        break;
      }
    pivotbincenter = pivotbin+0.5;

    double count_fraction = 
      2.0 * (ROOT::Math::normal_cdf(nsigmas)-ROOT::Math::normal_cdf(0));

    int nsideevents = 0;
    int lastbinindex = (int)events.size() - 1;

    if (nsigmas >= 0) nsideevents = lastbinindex - pivotbin;
    else nsideevents = pivotbin;

    int boundIdx = pivotbincenter + count_fraction*(double)nsideevents;
    int index = 0;

    if(nsigmas >= 0) index = min(boundIdx, (int)events.size()-1);
    else index = max(boundIdx, 0);

    return events.at(index);
  }



//Return 1 sigma bands from Nominal From Multiverse
//1 sigma corresponds to 68% of all universes 
std::vector<TH3F*> CalculateOneSigmaBandFromMultiverse(TH3F* nominal,
                                                       std::vector<TH3F*> vHist)
{
  TH3F* hResultUpSigma   = (TH3F*)nominal->Clone(ana::UniqueName().c_str());
  TH3F* hResultDownSigma = (TH3F*)nominal->Clone(ana::UniqueName().c_str());  
  
  for(int  xbin= 1; xbin <= hResultUpSigma->GetXaxis()->GetNbins(); xbin++){
    for(int  ybin= 1; ybin <= hResultUpSigma->GetYaxis()->GetNbins(); ybin++){
      for(int zbin = 1; zbin <= hResultUpSigma->GetZaxis()->GetNbins(); zbin++){
        
        double cv_value = nominal->GetBinContent(xbin,ybin,zbin);
        std::vector<double> multiverse_vals;

        for(uint iUniv = 0; iUniv < vHist.size(); iUniv++){
          multiverse_vals.push_back(vHist[iUniv]->GetBinContent
                                    (xbin,ybin,zbin));
        }
        
        std::sort(multiverse_vals.begin(),multiverse_vals.end(),
                  sort_in_ascending_order);

        
        float one_sigma_up   =  BinSigma(multiverse_vals,1,cv_value);
        float one_sigma_down =  BinSigma(multiverse_vals,-1,cv_value);

        hResultUpSigma->SetBinContent(xbin,ybin,zbin,one_sigma_up);
        hResultDownSigma->SetBinContent(xbin,ybin,zbin,one_sigma_down);
      }
    }
  } 

  return {(TH3F*)hResultUpSigma->Clone(ana::UniqueName().c_str()),
      (TH3F*)hResultUpSigma->Clone(ana::UniqueName().c_str())};

}

std::vector<TH2F*> CalculateOneSigmaBandFromMultiverse(TH2F* nominal,
                                                       std::vector<TH2F*> vHist)
{
  TH2F* hResultUpSigma   = (TH2F*)nominal->Clone(ana::UniqueName().c_str());
  TH2F* hResultDownSigma = (TH2F*)nominal->Clone(ana::UniqueName().c_str());  
  
  for(int  xbin= 1; xbin <= hResultUpSigma->GetXaxis()->GetNbins(); xbin++){
    for(int  ybin= 1; ybin <= hResultUpSigma->GetYaxis()->GetNbins(); ybin++){
        double cv_value = nominal->GetBinContent(xbin,ybin);
        std::vector<double> multiverse_vals;

        for(uint iUniv = 0; iUniv < vHist.size(); iUniv++){
          multiverse_vals.push_back(vHist[iUniv]->GetBinContent
                                    (xbin,ybin));
        }
        
        std::sort(multiverse_vals.begin(),multiverse_vals.end(),
                  sort_in_ascending_order);

        
        float one_sigma_up   =  BinSigma(multiverse_vals,1,cv_value);
        float one_sigma_down =  BinSigma(multiverse_vals,-1,cv_value);

        hResultUpSigma->SetBinContent(xbin,ybin,one_sigma_up);
        hResultDownSigma->SetBinContent(xbin,ybin,one_sigma_down);
    }
  } 

  return {(TH2F*)hResultUpSigma->Clone(ana::UniqueName().c_str()),
      (TH2F*)hResultUpSigma->Clone(ana::UniqueName().c_str())};
}



void MakeTemplatePlotsSyst(TH1F* hData,
                           std::vector<TH1F*> hTemplates,
                           std::vector<TH1F*> vUp,std::vector<TH1F*> vDown,
                           std::string plotType, std::string xTitle,
                           std::string dataName, bool isFitted,
                           int xbin, int ybin)
{
  hData->SetMarkerStyle(kChnsInfo[0].marker_style);
  hData->SetLineColor(kChnsInfo[0].color);
  hData->SetMarkerColor(kChnsInfo[0].color);

  for(int i = 0; i < (int)hTemplates.size(); i++){
    hTemplates[i]->SetLineColor(kChnsInfo[i+1].color);
  }


  hTemplates.front()->GetXaxis()->SetTitle(xTitle.c_str());
  hTemplates.front()->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
  hTemplates.front()->GetXaxis()->CenterTitle();
  hTemplates.front()->GetYaxis()->CenterTitle();

  hData->SetTitle("");
  hData->GetXaxis()->SetTitle(xTitle.c_str());
  hData->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
  hData->GetXaxis()->CenterTitle();
  hData->GetYaxis()->CenterTitle();
  
  TH1F* nominal_clone = (TH1F*)hTemplates[0]->Clone(ana::UniqueName().c_str());
  
  TGraphAsymmErrors* g = PlotSystErrorBand(nominal_clone,
                                           vUp, vDown,-1,-1,1.3,false);
  
  // std::cout << g->GetErrorYhigh(2) << "\t" << g->GetErrorYlow(2) << std::endl;
  //std::cout << nominal_clone->GetBinContent(2) << std::endl;



  TH1F* hRatio = (TH1F*)hData->Clone(ana::UniqueName().c_str());
  hRatio->Divide(hTemplates[0]);
  hRatio->SetMarkerStyle(20);
  hRatio->SetLineColor(kBlack);
  hRatio->GetYaxis()->SetTitle("Data/MC");

  //Systematic Ratio
  std::vector<TH1F*> vUpRatio;
  std::vector<TH1F*> vDownRatio;  
  TH1F* nominal_clone_2 = (TH1F*)hTemplates[0]->Clone(ana::UniqueName().c_str());
  nominal_clone_2->Divide(nominal_clone);
  for(uint i = 0; i < vUp.size(); i++){
    vUpRatio.push_back( (TH1F*)vUp[i]->Clone(ana::UniqueName().c_str()));
    vDownRatio.push_back( (TH1F*)vDown[i]->Clone(ana::UniqueName().c_str()));
    
    vUpRatio[i]->Divide(nominal_clone);
    vDownRatio[i]->Divide(nominal_clone);
  }

  TGraphAsymmErrors* g_ratio = PlotSystErrorBand(nominal_clone_2,
                                                 vUpRatio, vDownRatio,
                                                 -1,-1,1.3,false);
  //Signal Like Templates
  hTemplates[1]->Add(hTemplates[2]);
  hTemplates[1]->Add(hTemplates[7]);

  //NumuCC Like Templates
  hTemplates[3]->Add(hTemplates[4]);
  

  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
  c1->SetLeftMargin(0.10);
  c1->SetRightMargin(0.10);
  TPad *pad1 = new TPad(ana::UniqueName().c_str(), "pad1",0,0,1,1);
  TPad *pad2 = new TPad(ana::UniqueName().c_str(), "pad2",0,0,1,1);
  pad1->SetFillStyle(0);
  pad2->SetFillStyle(0);
  pad1->SetBottomMargin(0.30);
  pad2->SetTopMargin(1-0.30);
  pad2->SetGridx();
  pad2->SetGridy();
  c1->cd();
  pad1->Draw();
  pad1->cd();
  hTemplates[0]->GetXaxis()->SetLabelColor(kWhite);
  hTemplates[0]->GetYaxis()->SetLabelSize(0.03);
  hTemplates[0]->GetYaxis()->SetTitleSize(0.035);
  hTemplates[0]->Draw("hist");
  g->Draw("e2 same");
  hTemplates[0]->Draw("hist same");
  hTemplates[1]->Draw("hist same");
  hTemplates[3]->Draw("hist same");
  hTemplates[5]->Draw("hist same");
  hTemplates[6]->Draw("hist same");
  hData->Draw("same");
  TLegend* leg = ana::AutoPlaceLegend(0.2,0.2);
  if(isFitted)leg->SetHeader("Template Weights Applied");
  leg->AddEntry(hData, "Fake Data", "p");
  leg->AddEntry(hTemplates[0], "Nominal MC", "l");
  leg->AddEntry(hTemplates[1], "Signal-Like Template", "l");
  leg->AddEntry(hTemplates[3], "#nu_{#mu} CC-Like Template", "l");
  leg->AddEntry(hTemplates[5], "NC Template", "l");
  leg->AddEntry(hTemplates[6], "Other Template", "l");
  leg->Draw();
  SimulationSide();
  gPad->RedrawAxis();
  c1->cd();
  pad2->Draw();
  pad2->cd();
  hRatio->GetYaxis()->SetLabelSize(0.03);
  hRatio->GetYaxis()->SetTitleSize(0.035); 
  hRatio->GetYaxis()->SetRangeUser(0.4,1.5);
  hRatio->Draw("same");
  g_ratio->Draw("e2 same");
  hRatio->Draw("same");
  pad2->Draw("same");
  gPad->RedrawAxis();
  c1->cd();
  c1->Update();
  std::string fitted_string = "nominal_prediction";
  if(isFitted)fitted_string = "template_fit_prediction";
  char outname[50];
  sprintf(outname, "%s%s_%s_%s_%i_%i.png","Plots/", dataName.c_str(),
          plotType.c_str(), fitted_string.c_str(), xbin,ybin);
  c1->SaveAs(outname);
  c1->Close();

}

void PlotFrom2D(TH2F* hData2D,TH2F* hNominal2D,
                std::vector<TH2F*> hFitResults2D,
                std::vector<TH2F*> hSystsUp2D,
                std::vector<TH2F*> hSystsDown2D,                   
                std::string plotType, std::string xTitle,
                std::string dataName)
{
  for(int xbin = 1; xbin <= hData2D->GetXaxis()->GetNbins(); xbin++){
    TH1F* hData = 
      (TH1F*)hData2D->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);

    if(hData->Integral() == 0) continue;

    TH1F* hNominal = 
      (TH1F*)hNominal2D->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
    std::vector<TH1F*> hFitResults1D;
    for(uint i = 0; i < hFitResults2D.size(); i++){
      hFitResults1D.push_back((TH1F*)hFitResults2D[i]->ProjectionY(ana::UniqueName().c_str(),xbin,xbin));
    }
    std::vector<TH1F*> hSystsUp1D;
    std::vector<TH1F*> hSystsDown1D;
    for(uint i = 0; i < hSystsUp2D.size(); i++){
      hSystsUp1D.push_back((TH1F*)hSystsUp2D[i]->ProjectionY(ana::UniqueName().c_str(),xbin,xbin));
      hSystsDown1D.push_back((TH1F*)hSystsDown2D[i]->ProjectionY(ana::UniqueName().c_str(),xbin,xbin));
    }

  hData->SetMarkerStyle(kChnsInfo[0].marker_style);
  hData->SetLineColor(kChnsInfo[0].color);
  hData->SetMarkerColor(kChnsInfo[0].color);
  
  hFitResults1D.front()->SetLineColor(kBlue);
  hFitResults1D.front()->SetMarkerColor(kBlue);

  hNominal->SetLineColor(kRed);
  hNominal->SetMarkerColor(kRed);

  hFitResults1D.front()->GetXaxis()->SetTitle(xTitle.c_str());
  hFitResults1D.front()->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
  hFitResults1D.front()->GetXaxis()->CenterTitle();
  hFitResults1D.front()->GetYaxis()->CenterTitle();

  hData->SetTitle("");
  hData->GetXaxis()->SetTitle(xTitle.c_str());
  hData->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
  hData->GetXaxis()->CenterTitle();
  hData->GetYaxis()->CenterTitle();

  TH1F* nominal_clone = (TH1F*)hNominal->Clone(ana::UniqueName().c_str());
  
  TGraphAsymmErrors* g = PlotSystErrorBand(nominal_clone,
                                           hSystsUp1D, hSystsDown1D,
                                           -1,-1,1.3,false);

  TH1F* hRatio = (TH1F*)hData->Clone(ana::UniqueName().c_str());
  hRatio->Divide(hNominal);
  hRatio->SetMarkerStyle(20);
  hRatio->SetLineColor(kRed);
  hRatio->GetYaxis()->SetTitle("Data/MC");

  TH1F* hRatio2 = (TH1F*)hData->Clone(ana::UniqueName().c_str());
  hRatio->Divide(hFitResults1D[0]);
  hRatio->SetMarkerStyle(20);
  hRatio->SetLineColor(kBlue);
  hRatio->GetYaxis()->SetTitle("Data/MC");

  
  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
  c1->SetLeftMargin(0.10);
  c1->SetRightMargin(0.10);
  TPad *pad1 = new TPad(ana::UniqueName().c_str(), "pad1",0,0,1,1);
  TPad *pad2 = new TPad(ana::UniqueName().c_str(), "pad2",0,0,1,1);
  pad1->SetFillStyle(0);
  pad2->SetFillStyle(0);
  pad1->SetBottomMargin(0.30);
  pad2->SetTopMargin(1-0.30);
  pad2->SetGridx();
  pad2->SetGridy();
  c1->cd();
  pad1->Draw();
  pad1->cd();
  hData->GetXaxis()->SetLabelColor(kWhite);
  hData->GetYaxis()->SetLabelSize(0.03);
  hData->GetYaxis()->SetTitleSize(0.035);
  hData->GetXaxis()->SetRangeUser(1.0,10.0);
  hData->Draw();
  hNominal->Draw("hist same");
  g->Draw("e2 same");
  hNominal->Draw("hist same");
  hFitResults1D.front()->Draw("same");
  hData->Draw("same");
  hFitResults1D.front()->Draw("same");
  TLegend* leg = ana::AutoPlaceLegend(0.2,0.2);
  leg->AddEntry(hData, "Fake Data", "p");
  leg->AddEntry(hNominal, "Nominal MC", "l");
  leg->AddEntry(hFitResults1D[0], "Weighted MC", "pl");
  leg->Draw();
  SimulationSide();
  gPad->RedrawAxis();
  c1->cd();
  pad2->Draw();
  pad2->cd();
  hRatio->GetXaxis()->SetRangeUser(1.0,10.0);
  hRatio->GetYaxis()->SetLabelSize(0.03);
  hRatio->GetYaxis()->SetTitleSize(0.035); 
  hRatio->GetYaxis()->SetRangeUser(0.4,1.5);
  hRatio->Draw("same");
  hRatio2->Draw("same");
  pad2->Draw("same");
  gPad->RedrawAxis();
  c1->cd();
  c1->Update();
  char outname[50];
  sprintf(outname, "%s%s_%s_%i.png","Plots/", dataName.c_str(),
          plotType.c_str(), xbin);
  c1->SaveAs(outname);
  c1->Close();
  }//loop over xbins 
}

void PlotWithRatio2D(TH2F* hResultA, TH2F* hTruthA, TH2F* hNominalA,
                     std::string ytitle, std::string xtitle,
                     float xmin, float xmax,
                     std::string dataName, 
                     std::string plotType, std::string plotName, 
                     std::string varName, 
                     std::string legend_results = "Fake Data",
                     std::string legend_truth   = "Fake Data Truth",
                     std::string legend_nominal = "Nominal MC")
{


  for(int xbin = 1; xbin <= hResultA->GetXaxis()->GetNbins(); xbin++){

    TH1F* hResult = 
      (TH1F*)hResultA->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
    TH1F* hTruth = 
      (TH1F*)hTruthA->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);
    TH1F* hNominal =
      (TH1F*)hNominalA->ProjectionY(ana::UniqueName().c_str(),xbin,xbin);      

    if(hTruth->Integral() == 0) continue;

    std::stringstream low_X;
    low_X << std::fixed << std::setprecision(2) << 
      hResultA->GetXaxis()->GetBinLowEdge(xbin);
    std::stringstream hi_X;
    hi_X << std::fixed << std::setprecision(2) << 
      hResultA->GetXaxis()->GetBinUpEdge(xbin);
    
    std::string caption = low_X.str() + " #leq cos #theta_{e} < "+
      hi_X.str();

    
    hResult->GetYaxis()->SetTitle(ytitle.c_str());
    hResult->GetXaxis()->SetTitle(xtitle.c_str());
    hResult->SetTitle(caption.c_str());
    hResult->GetYaxis()->CenterTitle();
    hResult->GetXaxis()->CenterTitle();
    hResult->GetYaxis()->SetTitleOffset(1.15);
    hResult->SetMarkerStyle(20);
    hResult->SetMarkerColor(kBlack);
    hTruth->SetLineColor(kRed);
    hTruth->SetMarkerColor(kRed);

    hNominal->SetLineColor(kBlue);
    hNominal->SetMarkerColor(kBlue);

    TH1F* hResultClone = (TH1F*)hResult->Clone(ana::UniqueName().c_str());
    TH1F* hNominalClone = (TH1F*)hNominal->Clone(ana::UniqueName().c_str());
    hResultClone->SetTitle("");
    hResultClone->Divide(hTruth);
    hNominalClone->Divide(hTruth);

    hResult->GetXaxis()->SetRangeUser(xmin,xmax);
    hResultClone->GetXaxis()->SetRangeUser(xmin,xmax);
    

    //Find Y limits
    double ymax = hResult->GetMaximum();

    if(ymax < hTruth->GetMaximum())
      hResult->GetYaxis()->SetRangeUser(0,hTruth->GetMaximum()*1.3);
    else
      hResult->GetYaxis()->SetRangeUser(0,ymax*1.3);


    TCanvas *c1 = new TCanvas(ana::UniqueName().c_str(), 
                              ana::UniqueName().c_str());
    c1->SetLeftMargin(0.10);
    c1->SetRightMargin(0.10);
    c1->SetBottomMargin(0.10);
    TPad *pad1 = new TPad(ana::UniqueName().c_str(), "pad1",0,0,1,1);
    TPad *pad2 = new TPad(ana::UniqueName().c_str(), "pad2",0,0,1,1);
    pad1->SetFillStyle(0);
    pad2->SetFillStyle(0);
    pad1->SetBottomMargin(0.30);
    pad2->SetTopMargin(1-0.30);
    pad2->SetGridx();
    pad2->SetGridy();
    c1->cd();
    pad1->Draw();
    pad1->cd();
    hResult->GetXaxis()->SetLabelColor(kWhite);
    hResult->GetYaxis()->SetLabelSize(0.03);
    hResult->GetYaxis()->SetTitleSize(0.035);
    hResult->Draw();
    //    PrintChiSq(caption_chi.c_str());
    hResult->Draw("same");
    hTruth->Draw("hist same");
    hNominal->Draw("same");
    TLegend* leg = ana::AutoPlaceLegend(0.23,0.23);
    leg->AddEntry(hResult, legend_results.c_str(), "p");
    leg->AddEntry(hTruth,  legend_truth.c_str(), "l");
    leg->AddEntry(hNominal, legend_nominal.c_str(), "p");
    leg->Draw();
    SimulationSide();
    gPad->RedrawAxis();
    c1->cd();
    pad2->Draw();
    pad2->cd();
    hResultClone->GetYaxis()->SetTitle("Ratio to Truth");
    hResultClone->GetYaxis()->SetLabelSize(0.03);
    hResultClone->GetYaxis()->SetTitleSize(0.035); 
    hResultClone->GetYaxis()->SetRangeUser(0.4,1.5);
    hResultClone->Draw("same");
    hNominalClone->Draw("same");
    pad2->Draw("same");
    gPad->RedrawAxis();
    c1->cd();
    c1->Update();
    char outname[50];
    sprintf(outname, "%s%s_%s_%s_%s_%i.png","Plots/", dataName.c_str(),
            plotType.c_str(), plotName.c_str(), varName.c_str(),xbin);
    c1->SaveAs(outname);
    c1->Close();

    delete c1;
    delete leg;
    delete hResult;
    delete hTruth;
    delete hResultClone;

  }//loop over xbins

}


void MakeCovarianceCorrelationPlots(TH2F* hCov, TH2F* hCorr, 
                                    int xbin, int ybin)
{
  hCov->GetXaxis()->CenterTitle();
  hCov->GetYaxis()->CenterTitle();
  hCov->GetZaxis()->CenterTitle();

  hCorr->GetXaxis()->CenterTitle();
  hCorr->GetYaxis()->CenterTitle();
  hCorr->GetZaxis()->CenterTitle();

  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
  c1->SetLeftMargin(0.15);
  c1->SetRightMargin(0.15);
  c1->SetBottomMargin(0.15);
  gStyle->SetTitleAlign(33);
  hCov->Draw("COLZ");
  Simulation();
  char outname[50];
  sprintf(outname, "%s%s_%s_%i_%i.png","Plots/",
          "covariance", "template_bins", xbin,ybin);
  c1->SaveAs(outname);
  c1->Close();
  delete c1;

  TCanvas *c2 = new TCanvas(ana::UniqueName().c_str());
  c2->SetLeftMargin(0.15);
  c2->SetRightMargin(0.15);
  c2->SetBottomMargin(0.15);
  gStyle->SetTitleAlign(33);
  hCorr->Draw("COLZ");
  Simulation();
  sprintf(outname, "%s%s_%s_%i_%i.png","Plots/",
          "correlation", "template_bins", xbin,ybin);
  c2->SaveAs(outname);
  c2->Close();
  delete c2;

}



void doTemplateFit(std::string dataName = "nominal",
                   bool shouldSave = true)
{
  TFile* fTemplates = new TFile((sDir+sTemplates).c_str(),"read");

  //////////////////////////////////////////////////////////////////////////
  //////// Grab Systematic Templates For Covariance Matrix /////////////////
  //////////////////////////////////////////////////////////////////////////

  NueCCIncCrossSectionAnalysis nominal_ana =
    *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("Nominal/Analysis"));
  NueCCIncCrossSectionAnalysis light_dwn_ana =
    *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("LightDwn/Analysis"));
  NueCCIncCrossSectionAnalysis light_up_ana =
    *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("LightUp/Analysis"));
  NueCCIncCrossSectionAnalysis calib_dwn_ana =
    *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("CalibDwn/Analysis"));
  NueCCIncCrossSectionAnalysis calib_up_ana =
    *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("CalibUp/Analysis"));
  NueCCIncCrossSectionAnalysis calib_shape_ana =
    *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("CalibShape/Analysis"));
  NueCCIncCrossSectionAnalysis ckv_ana =
    *NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory("Ckv/Analysis"));

  std::vector<std::unique_ptr<NueCCIncCrossSectionAnalysis>> genie_ana;
  std::vector<std::unique_ptr<NueCCIncCrossSectionAnalysis>> ppfx_ana;

  for(int i = 0; i < 100; i++){
    genie_ana.push_back(NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory(Form("GENIE/%s_%i","Universe",i))));
    ppfx_ana.push_back(NueCCIncCrossSectionAnalysis::LoadFrom(fTemplates->GetDirectory(Form("PPFX/%s_%i","Universe",i))));
  }

  TH3F* nominal_template   = nominal_ana.getTemplate3D(costhetabins,
                                                       eelecbins,
                                                       pidbins);

  TH3F* light_dwn_template = light_dwn_ana.getTemplate3D(costhetabins,
                                                         eelecbins,
                                                         pidbins);

  TH3F* light_up_template  = light_up_ana.getTemplate3D(costhetabins,
                                                        eelecbins,
                                                        pidbins);

  TH3F* calib_dwn_template = calib_dwn_ana.getTemplate3D(costhetabins,
                                                         eelecbins,
                                                         pidbins);

  TH3F* calib_up_template  = calib_up_ana.getTemplate3D(costhetabins,
                                                        eelecbins,
                                                        pidbins);

  TH3F* calib_shape_template = calib_shape_ana.getTemplate3D(costhetabins,
                                                             eelecbins,
                                                             pidbins);

  TH3F* ckv_template  = ckv_ana.getTemplate3D(costhetabins,
                                              eelecbins,
                                              pidbins);

  std::vector<TH3F*> genie_templates;
  std::vector<TH3F*> ppfx_templates;

  for(int i = 0; i < 100; i++){
    genie_templates.push_back(genie_ana[i]->getTemplate3D(costhetabins,
                                                          eelecbins,
                                                          pidbins));
    ppfx_templates.push_back(ppfx_ana[i]->getTemplate3D(costhetabins,
                                                        eelecbins,
                                                        pidbins));
  }
  
  std::vector<TH3F*> vMultiverse;
  for(int i = 0; i < 100; i++){
    vMultiverse.push_back((TH3F*)genie_templates[i]->Clone
                          (ana::UniqueName().c_str()));
  }
  for(int i = 0; i < 100; i++){
    vMultiverse.push_back((TH3F*)ppfx_templates[i]->Clone
                          (ana::UniqueName().c_str()));
  }


  //////////////////////////////////////////////////////////////////////////
  /////////////////// Grab Templates From Nominal///////////////////////////
  //////////////////////////////////////////////////////////////////////////
  NueCCIncCrossSectionTemplates fit_templates = 
    *NueCCIncCrossSectionTemplates::LoadFrom(fTemplates->GetDirectory("Nominal/Templates"));
   
  std::vector<TH3F*>nominal_fit_templates = 
    fit_templates.getAllTemplates3D(costhetabins,
                                    eelecbins,
                                    pidbins);

  std::vector<TH3F*>nominal_fit_analysis = 
    fit_templates.getAllAnalysis3D(costhetabins,
                                   eelecbins,
                                   enubins);

  std::vector<TH3F*>nominal_analysis_comparison;
  for(uint i = 0; i < nominal_fit_analysis.size(); i++){
    nominal_analysis_comparison.push_back((TH3F*)
                                          nominal_fit_analysis[i]->Clone
                                          (ana::UniqueName().c_str()));
  }
  
  
  NueCCIncTemplateFitter
    fitter((TH3F*)nominal_template->Clone(ana::UniqueName().c_str()),
           nominal_fit_templates,
           nominal_fit_analysis,
           (TH3F*)nominal_template->Clone(ana::UniqueName().c_str()),
           {(TH3F*)light_up_template->Clone(ana::UniqueName().c_str()),
               (TH3F*)calib_up_template->Clone(ana::UniqueName().c_str())},
           {(TH3F*)light_dwn_template->Clone(ana::UniqueName().c_str()),
               (TH3F*)calib_dwn_template->Clone(ana::UniqueName().c_str())},
           {(TH3F*)ckv_template->Clone(ana::UniqueName().c_str()),
               (TH3F*)calib_shape_template->Clone(ana::UniqueName().c_str())},
           vMultiverse,
           21);

  std::vector<TH3F*> hFitResults = fitter.doFullFit();


  //////////////////////////////////////////////////////////////////////////
  /////////////////// Plot Template Results ////////////////////////////////
  //////////////////////////////////////////////////////////////////////////

  if(plotFitTemplateResults){
    std::vector<TH3F*> hGenieSysts3D = 
      CalculateOneSigmaBandFromMultiverse(nominal_template,
                                          genie_templates);
    
    std::vector<TH3F*> hPPFXSysts3D = 
      CalculateOneSigmaBandFromMultiverse(nominal_template,
                                          ppfx_templates);
    
    //Plot Template Fit Results
    for(int xbin = 1; xbin <= hFitResults[0]->GetXaxis()->GetNbins(); xbin++){
      for(int ybin = 1; ybin <= hFitResults[0]->GetYaxis()->GetNbins(); ybin++){

        std::vector<TH1F*> hPID_Nominal = fitter.getPIDTemplates(xbin,ybin);
        std::vector<TH1F*> hPID_Fitted = fitter.getPIDFitTemplates(xbin,ybin);
        TH2F* hCovariance = fitter.getCovarianceMatrixTemplateBins(xbin,ybin);
        TH2F* hCorrelation = fitter.getCorrelationMatrixTemplateBins(xbin,ybin);
        
        if(hPID_Nominal[0]->Integral() == 0) continue;
        
        std::vector<TH1F*> hSystUpHolder1D = {
          (TH1F*)light_up_template->ProjectionZ(ana::UniqueName().c_str(),
                                                xbin,xbin,ybin,ybin),
          (TH1F*)calib_up_template->ProjectionZ(ana::UniqueName().c_str(),
                                                xbin,xbin,ybin,ybin),
          (TH1F*)ckv_template->ProjectionZ(ana::UniqueName().c_str(),
                                           xbin,xbin,ybin,ybin),
          (TH1F*)calib_shape_template->ProjectionZ(ana::UniqueName().c_str(),
                                                   xbin,xbin,ybin,ybin),
          (TH1F*)hGenieSysts3D[0]->ProjectionZ(ana::UniqueName().c_str(),
                                               xbin,xbin,ybin,ybin),
          (TH1F*)hPPFXSysts3D[0]->ProjectionZ(ana::UniqueName().c_str(),
                                              xbin,xbin,ybin,ybin)
        };
        
        std::vector<TH1F*> hSystDownHolder1D = {
          (TH1F*)light_dwn_template->ProjectionZ(ana::UniqueName().c_str(),
                                                 xbin,xbin,ybin,ybin),
          (TH1F*)calib_dwn_template->ProjectionZ(ana::UniqueName().c_str(),
                                                 xbin,xbin,ybin,ybin),
          (TH1F*)nominal_template->ProjectionZ(ana::UniqueName().c_str(),
                                               xbin,xbin,ybin,ybin),
          (TH1F*)nominal_template->ProjectionZ(ana::UniqueName().c_str(),
                                               xbin,xbin,ybin,ybin),
          (TH1F*)hGenieSysts3D[1]->ProjectionZ(ana::UniqueName().c_str(),
                                               xbin,xbin,ybin,ybin),
          (TH1F*)hPPFXSysts3D[1]->ProjectionZ(ana::UniqueName().c_str(),
                                              xbin,xbin,ybin,ybin)
        };
        
        
        MakeTemplatePlotsSyst((TH1F*)nominal_template->ProjectionZ
                              (ana::UniqueName().c_str(),xbin,xbin,ybin,ybin),
                              hPID_Nominal,
                              hSystUpHolder1D,
                              hSystDownHolder1D,       
                              "templates", "CVNe 2017", "Nominal",
                              false,xbin,ybin);
        
        MakeTemplatePlotsSyst((TH1F*)nominal_template->ProjectionZ
                              (ana::UniqueName().c_str(),xbin,xbin,ybin,ybin),
                              hPID_Fitted,
                              hSystUpHolder1D,
                              hSystDownHolder1D,       
                              "templates", "CVNe 2017", "Nominal",
                              true,xbin,ybin);
        
        MakeCovarianceCorrelationPlots(hCovariance, hCorrelation, xbin,ybin);
        
        /*
        delete hCovariance;
        delete hCorrelation;

        for (int i =0; i< hPID_Nominal.size();i++)
          {
            delete (hPID_Nominal[i]);
            delete (hPID_Fitted[i]);
          } 
        hPID_Nominal.clear();
        hPID_Fitted.clear();
        */
      }
    }
  } //if plot fit results

  //////////////////////////////////////////////////////////////////////////
  /////////////////// Plot Fit Results in Electron Kinematic Space//////////
  //////////////////////////////////////////////////////////////////////////
  //Project all analysis plots down to "yx" axis
  std::vector<TH2F*> hFitResults2D;
  for(uint i = 0; i < hFitResults.size();i++){
    TH2F* hHolder = (TH2F*)hFitResults[i]->Project3D("yx");
    hHolder->SetName(ana::UniqueName().c_str());
    hFitResults2D.push_back((TH2F*)hHolder->Clone(ana::UniqueName().c_str()));
  }

  TH2F* data_analysis_2D = nominal_ana.getAnalysis2D(costhetabins,
                                                     eelecbins,
                                                     enubins);


  TH2F* nominal_analysis_2D = nominal_ana.getAnalysis2D(costhetabins,
                                                        eelecbins,
                                                        enubins);
  if(plotFitAnalysisResults){
  TH2F* light_dwn_analysis_2D = light_dwn_ana.getAnalysis2D(costhetabins,
                                                            eelecbins,
                                                            enubins);

  TH2F* light_up_analysis_2D = light_up_ana.getAnalysis2D(costhetabins,
                                                            eelecbins,
                                                            enubins);

  TH2F* calib_dwn_analysis_2D = calib_dwn_ana.getAnalysis2D(costhetabins,
                                                            eelecbins,
                                                            enubins);

  TH2F* calib_up_analysis_2D = calib_up_ana.getAnalysis2D(costhetabins,
                                                          eelecbins,
                                                          enubins);

  TH2F* calib_shape_analysis_2D = calib_shape_ana.getAnalysis2D(costhetabins,
                                                                eelecbins,
                                                                enubins);

  TH2F* ckv_analysis_2D = ckv_ana.getAnalysis2D(costhetabins,
                                                eelecbins,
                                                enubins);

  std::vector<TH2F*> genie_analysis_2D;
  std::vector<TH2F*> ppfx_analysis_2D;

  for(int i = 0; i < 100; i++){
    genie_analysis_2D.push_back(genie_ana[i]->getAnalysis2D(costhetabins,
                                                           eelecbins,
                                                           enubins));
    ppfx_analysis_2D.push_back(ppfx_ana[i]->getAnalysis2D(costhetabins,
                                                         eelecbins,
                                                         enubins));
  }
  //if(plotFitAnalysisResults){

    std::vector<TH2F*> hGenieSysts2D = 
      CalculateOneSigmaBandFromMultiverse(nominal_analysis_2D,
                                          genie_analysis_2D);
    
    std::vector<TH2F*> hPPFXSysts2D = 
      CalculateOneSigmaBandFromMultiverse(nominal_analysis_2D,
                                          ppfx_analysis_2D);

        
        std::vector<TH2F*> hSystUpHolder2D = {
          (TH2F*)light_up_analysis_2D->Clone(ana::UniqueName().c_str()),
          (TH2F*)calib_up_analysis_2D->Clone(ana::UniqueName().c_str()),
          (TH2F*)ckv_analysis_2D->Clone(ana::UniqueName().c_str()),
          (TH2F*)calib_shape_analysis_2D->Clone(ana::UniqueName().c_str()),
          (TH2F*)hGenieSysts2D[0]->Clone(ana::UniqueName().c_str()),      
          (TH2F*)hPPFXSysts2D[0]->Clone(ana::UniqueName().c_str()),       
        };

        std::vector<TH2F*> hSystDownHolder2D = {
          (TH2F*)light_dwn_analysis_2D->Clone(ana::UniqueName().c_str()),
          (TH2F*)calib_dwn_analysis_2D->Clone(ana::UniqueName().c_str()),
          (TH2F*)nominal_analysis_2D->Clone(ana::UniqueName().c_str()),
          (TH2F*)nominal_analysis_2D->Clone(ana::UniqueName().c_str()),
          (TH2F*)hGenieSysts2D[1]->Clone(ana::UniqueName().c_str()),      
          (TH2F*)hPPFXSysts2D[1]->Clone(ana::UniqueName().c_str()),       
        };


    PlotFrom2D(data_analysis_2D,
               nominal_analysis_2D,
               hFitResults2D,
               hSystUpHolder2D,
               hSystDownHolder2D,                       
               "analysis", 
               "Reconstructed Electron Energy, E_{e} (GeV)", 
               "Nominal");
    
  }//if plot analysis fit results

  //////////////////////////////////////////////////////////////////////////
  /////////////////// Unfold Signal Estimate to True Space /////////////////
  //////////////////////////////////////////////////////////////////////////
  TH2F* hSignalEst = (TH2F*)hFitResults2D[1]->Clone("hSignalEst");

  TH2F* hUnfoldedSignalEst = nominal_ana.doUnfolding2D(hSignalEst, 6);


  TH2F* hTrueRecoSignal = 
    (TH2F*)nominal_ana.getSelectedSignalPrediction2D(costhetabins,
                                                     eelecbins);
  
  TH2F* hTrueRecoSignal_nominal = 
    (TH2F*)nominal_ana.getSelectedSignalPrediction2D(costhetabins,
                                                     eelecbins);



  if(plotSignalPrediction){
    PlotWithRatio2D(hUnfoldedSignalEst, hTrueRecoSignal,hTrueRecoSignal_nominal,
                    "Signal Events/8.09 #times 10^{20} POT", 
                    "True Electron Energy, E_{e} (GeV)",
                    1.0, 10.0,
                    dataName, 
                    "reco_signal", "comparison", 
                    "double", "Weighted MC - Signal Prediction",
                    "True Signal from Fake Data",
                    "Nominal MC - Signal Prediction");
  }

  //////////////////////////////////////////////////////////////////////////
  /////////////////// Efficiency  //////////////////////////////////////////
  //////////////////////////////////////////////////////////////////////////
  TH2F* nominal_efficiency_2D = 
    nominal_ana.getEfficiency2D(costhetabins,eelecbins);
  TH2F* light_dwn_efficiency_2D = 
    light_dwn_ana.getEfficiency2D(costhetabins,eelecbins);

  TH2F* light_up_efficiency_2D = 
    light_up_ana.getEfficiency2D(costhetabins,eelecbins);

  TH2F* calib_dwn_efficiency_2D = 
    calib_dwn_ana.getEfficiency2D(costhetabins,eelecbins);

  TH2F* calib_up_efficiency_2D = 
    calib_up_ana.getEfficiency2D(costhetabins,eelecbins);

  TH2F* calib_shape_efficiency_2D =
    calib_shape_ana.getEfficiency2D(costhetabins,eelecbins);

  TH2F* ckv_efficiency_2D = ckv_ana.getEfficiency2D(costhetabins,eelecbins);

  std::vector<TH2F*> genie_efficiency_2D;
  std::vector<TH2F*> ppfx_efficiency_2D;

  for(int i = 0; i < 100; i++){
    genie_efficiency_2D.push_back(genie_ana[i]->getEfficiency2D(costhetabins,
                                                           eelecbins));
    ppfx_efficiency_2D.push_back(ppfx_ana[i]->getEfficiency2D(costhetabins,
                                                              eelecbins));
  }

  //////////////////////////////////////////////////////////////////////////
  /////////////////// Efficiency - MRE Correction  /////////////////////////
  //////////////////////////////////////////////////////////////////////////




  std::cout << "DONE" << std::endl;
  fTemplates->Close();
}