MakePlots.C

Go to the documentation of this file.

//Script that does background estimation in pi^{0} kinetic energy  bins
//It reads the spectra generated from the getSpectra_ForFitting.C macro

#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 "3FlavorAna/Cuts/NueCutsSecondAna.h"
#include "3FlavorAna/Vars/NueVars.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"

#include "CAFAna/Systs/XSecSystLists.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "CAFAna/XSec/GenieMultiverseSyst.h"
#include "CAFAna/XSec/FluxMultiverseSyst.h"
#include "CAFAna/XSec/MultiverseCorrelation.h"


#include "NDAna/ncpi0_SemiInclusive/NDNCPi0Xsec.h"
#include "NDAna/ncpi0_SemiInclusive/NCPi0SemiIncCuts.h"
#include "NDAna/ncpi0_SemiInclusive/NCPi0SemiIncVars.h"
#include "NDAna/ncpi0_SemiInclusive/NCPi0SemiIncBins.h"

#include "Utilities/rootlogon.C"

#include "TH1.h"
#include "TH2.h"
#include "THStack.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TStyle.h"
#include "TFile.h"
#include <vector> 
#include <math.h>
#include "TLine.h"

#include <iostream>
#include <fstream>
#include <iomanip> 
#include <vector>
#include <cmath>
#include <math.h> 
#include <string>
#include <stdint.h>
#include "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#include <sstream>
#include <math.h>
#include <string>
#include <vector>
#include "TF2.h"
#include "TH2.h"
#include "TCutG.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TRandom.h"
#include "TGraph.h"
#include "TColor.h"
#include "TF3.h"
#include "TH3.h"
#include "TH3F.h"
#include "TMinuit.h"
#include "TLegend.h"
#include "THStack.h"
#include "TMatrix.h"
#include "TGraphErrors.h"
#include "TGraph.h"
//#include "TDecompLU.h"
#include "TDecompSVD.h"

using namespace ana;

const int xbins = 20;
const double pot = 8.09e20; //data pot

float fDA[xbins];
float fSig[xbins];
float fCCBkg[xbins];
float fNCBkg[xbins];

float eDA[xbins];
float eSig[xbins];
float eCCBkg[xbins];
float eNCBkg[xbins];

float covMatrix[xbins][xbins];
float stat_covMatrix[xbins][xbins];

const std::vector<int> kNueCCColorDef = {
  632, //kRed, Total MC
  600,  //kBlue Signal
  616,  //kMagenta nuebar CC
  416-3, //kGrenn - 3 NumuCC
};

void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){

  double chisquare  = 0;
  float data    = 0;
  float signal  = 0;
  float ccbkg  = 0;
  float ncbkg    = 0;

 
  float data_j    = 0;
  float signal_j  = 0;
  float ccbkg_j  = 0;
  float ncbkg_j    = 0;
 
  float error   = 0;

  float delta     = 0;
  float firstterm = 0;
  int dof_my = 10;
 
  //  Covariance Matrix Inversion First
  TMatrixF Vij_stat = TMatrix(xbins,xbins);
  TMatrixF Vij = TMatrix(xbins,xbins);
  for(int i =0; i < xbins; i++){
    for(int j = 0; j < xbins; j++){
      Vij_stat[i][j] = stat_covMatrix[i][j];
      Vij[i][j] = covMatrix[i][j];
    }
  } 

  //inversion
  TMatrixD H3 =  Vij_stat + Vij;
  TDecompSVD svd(H3);
  TMatrixD VijInv = svd.Invert();

  float delta_i = 0;
  float delta_j = 0;

  TMatrixD FitVector = TMatrixD(1,xbins);
  TMatrixD FitVector2 = TMatrixD(xbins,1);

  //Calculate Chi Sq
  for(int i =0; i < xbins; i++){
    data    = fDA[i];
    signal  = fSig[i];
    ccbkg    = fCCBkg[i];
    ncbkg    = fNCBkg[i];
    
    if(data < 20){
      FitVector[0][i] = 0;
      FitVector2[i][0] = 0;
      continue;

    }
    
    FitVector[0][i] = ( data - (par[0]*signal + par[1]*ccbkg + par[2]*ncbkg)); 
    FitVector2[i][0] = ( data -(par[0]*signal+ par[1]*ccbkg + par[2]*ncbkg));     
  }
  
  TMatrixD Chi = FitVector * VijInv;
  TMatrixD Chi2 = Chi * FitVector2;
  float secondterm = Chi2(0,0);
  
  chisquare = secondterm; 
  f = chisquare;
}

void fcn2Var(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){

  double chisquare  = 0;
  float data    = 0;
  float signal  = 0;
  float ccbkg  = 0;
  float ncbkg    = 0;
 

  float data_j    = 0;
  float signal_j  = 0;
  float ccbkg_j  = 0;
  float ncbkg_j    = 0;
 

  float error   = 0;

  float delta     = 0;
  float firstterm = 0;

  //  Covariance Matrix Inversion First
  TMatrixF Vij_stat = TMatrix(xbins,xbins);
  TMatrixF Vij = TMatrix(xbins,xbins);
  for(int i =0; i < xbins; i++){
    for(int j = 0; j < xbins; j++){
      Vij_stat[i][j] = stat_covMatrix[i][j];
      Vij[i][j] = covMatrix[i][j];
    }
  } 

 
  //Inversion
  TMatrixD H3 =  Vij_stat + Vij;
  TDecompSVD svd(H3);
  TMatrixD VijInv = svd.Invert();


  float delta_i = 0;
  float delta_j = 0;
  
  TMatrixD FitVector = TMatrixD(1,xbins);
  TMatrixD FitVector2 = TMatrixD(xbins,1);

  //Calculate Chi Sq
  for(int i =0; i < xbins; i++){
    data    = fDA[i];
    signal  = fSig[i];
    ccbkg    = fCCBkg[i];
    ncbkg    = fNCBkg[i]; 
    
    if(data < 20){
      FitVector[0][i] = 0;
      FitVector2[i][0] = 0;
      continue;

    }

    FitVector[0][i] = ( data - (par[0]*signal + par[1]*(ccbkg+ncbkg))); 
    FitVector2[i][0] = ( data -(par[0]*signal+ par[1]*(ccbkg+ncbkg)));  
    //dof_my++;   
  }
  
  TMatrixD Chi = FitVector * VijInv;
  TMatrixD Chi2 = Chi * FitVector2;
  float secondterm = Chi2(0,0);
  
  
  chisquare = secondterm; 

  f = chisquare;
 
}


void PrintChiSq(TString str)
{

  TLatex* prelim = new TLatex(.55, .80, str);
  prelim->SetTextColor(kBlack);
  prelim->SetNDC();
  prelim->SetTextSize(0.05);
  //prelim->SetTextAlign(32);
  prelim->Draw();
}
void PrintChiSq2(TString str)
{

  TLatex* prelim = new TLatex(.20, .75, str);
  prelim->SetTextColor(kBlack);
  prelim->SetNDC();
  prelim->SetTextSize(0.05);
  //prelim->SetTextAlign(32);
  prelim->Draw();
}

TGraphAsymmErrors* PlotWithSystErrorBand_mine(TH1F*& nom,
                                              std::vector<TH1F*>& ups,
                                              std::vector<TH1F*>& dns,
                                              int col, int errCol,
                                              float headroom, bool newaxis)
{
  if(col == -1){
    col = kRed;
    errCol = kRed - 9;
  }
  else if(errCol == -1) errCol = col-9; // hopefully a lighter version
  
  nom->SetLineColor(col);
  nom->GetXaxis()->CenterTitle();
  nom->GetYaxis()->CenterTitle();
  //if(newaxis) nom->Draw("hist ]["); // Set the axes up
  
  double yMax = nom->GetBinContent(nom->GetMaximumBin());

  TGraphAsymmErrors* g = new TGraphAsymmErrors;
  
  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
    const double y = nom->GetBinContent(binIdx);
    g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
    
    const double w = nom->GetXaxis()->GetBinWidth(binIdx);
    
    double errUp = 0;
    double errDn = 0;
      
    for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
      double hi = ups[systIdx]->GetBinContent(binIdx)-y;
      double lo = y-dns[systIdx]->GetBinContent(binIdx);
        
        
      if(lo <= 0 && hi <= 0) std::swap(lo, hi);
        
      errUp += hi*hi;
      errDn += lo*lo;
        
      // TODO: what happens if they're both high or both low?
    } // end for systIdx
      
    g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
  } // end for i
  
  g->SetFillColor(errCol);
  
  for(TH1* up: ups) delete up;
  for(TH1* dn: dns) delete dn;
  return g;
}

void MakePlots(TH2F* data, 
               std::vector<TH2F*> mc,
               TH1F* sig_weights, TH1F* cc_weights,
               TH1F* nc_weights, 
               std::vector<TH2F*>shifts_up,
               std::vector<TH2F*>shifts_down,
               std::string dataname,
               std::string weightname)
{
  //PID Axis Plots
  for(int x = 1; x <= data->GetYaxis()->GetNbins(); x++){

    std::stringstream low;      
    low<<std::fixed<<std::setprecision(2)<<
      data->GetYaxis()->GetBinLowEdge(x);
    std::stringstream high;
    high<<std::fixed<<std::setprecision(2)<<
      data->GetYaxis()->GetBinUpEdge(x);
    std::string caption = low.str() + " #leq KE_{#pi^{0}} < " + high.str();

    TH1F* dataproj = (TH1F*)data->ProjectionX("dataproj",x,x);
    dataproj->Rebin(2);    
    std::vector<TH1F*> shift_up_proj;
    std::vector<TH1F*> shift_down_proj;
    std::vector<TH1F*> shift_up_proj_ratio;
    std::vector<TH1F*> shift_down_proj_ratio;

    

    for(uint j = 0; j < shifts_up.size(); j++){
      char name_up[50];
      char name_down[50];      
      sprintf(name_up, "%s_%s_%i", "proj", "up", j);
      sprintf(name_down, "%s_%s_%i", "proj", "down", j);
      shift_up_proj.push_back((TH1F*)shifts_up[j]->ProjectionX(name_up,x,x));
      shift_down_proj.push_back((TH1F*)shifts_down[j]->ProjectionX(name_down,x,x));
      sprintf(name_up, "%s_%s_%i", "proj", "up_ratio", j);
      sprintf(name_down, "%s_%s_%i", "proj", "down_ratio", j);
      shift_up_proj_ratio.push_back((TH1F*)shifts_up[j]->ProjectionX(name_up,x,x));
      shift_down_proj_ratio.push_back((TH1F*)shifts_down[j]->ProjectionX(name_down,x,x));

      shift_up_proj[j]->Rebin(2);
      shift_down_proj[j]->Rebin(2);
      shift_up_proj_ratio[j]->Rebin(2);
      shift_down_proj_ratio[j]->Rebin(2);
    }

    std::vector<TH1F*> hMC;
    std::vector<TH1F*> hMC_wgt;
    for(uint i = 0; i < mc.size(); i++){
      char name[50];
      char name_wgt[50];

      sprintf(name, "%s_%i", "nominal", i);
      sprintf(name_wgt, "%s_%i", "weight", i);      

      hMC.push_back((TH1F*)mc[i]->ProjectionX(name, x,x));
      hMC_wgt.push_back((TH1F*)mc[i]->ProjectionX(name_wgt, x,x));
      hMC[i]->Rebin(2);
      hMC_wgt[i]->Rebin(2);
    }
    
    hMC_wgt[1]->Scale(sig_weights->GetBinContent(x));
    hMC_wgt[6]->Scale(cc_weights->GetBinContent(x));
    hMC_wgt[7]->Scale(nc_weights->GetBinContent(x));    

    TH1F *hTotalMC_wgt = (TH1F*)hMC_wgt[1]->Clone("hTotalMC_wgt");
    hTotalMC_wgt->Add(hMC_wgt[6]);
    hTotalMC_wgt->Add(hMC_wgt[7]);
    hTotalMC_wgt->SetLineColor(kBlue);

 
    TH1F* hBkgNom = (TH1F*)hMC[6]->Clone("hBkgNom");
    hBkgNom->Add(hMC[7]);
    TH1F* hBkgScl = (TH1F*)hMC_wgt[6]->Clone("hBkgScl");
    hBkgScl->Add(hMC_wgt[7]);
  
    TH1F* hSigNom = (TH1F*)hMC[1]->Clone("hSigNom");
    TH1F* hSigScl = (TH1F*)hMC_wgt[1]->Clone("hSigScl");
   
    hSigNom->SetLineColor(30); // sig before scaling
    hBkgNom->SetLineColor(kOrange-3);  // bkg before scaling

    hSigScl->SetLineColor(3); // sig after scaling
    hBkgScl->SetLineColor(kOrange+8);  // bkg after scaling
   

    TH1F* hDataRatio = (TH1F*) hMC[0]->Clone("data_ratio");
    TH1F* hDataRatio_Wgt = (TH1F*)hTotalMC_wgt->Clone("data_ratio_wgt");

    hDataRatio_Wgt->SetTitle("");
    hDataRatio->GetYaxis()->SetRangeUser(0.1,1.9);
   
    hDataRatio->Divide(dataproj);
    hDataRatio_Wgt->Divide(dataproj);

    TH1F* hSignalRatio = (TH1F*)hSigScl->Clone("hSignalRatio");
    hSignalRatio->Divide(hSigNom);

    TH1F* hBkgRatio = (TH1F*)hBkgScl->Clone("hBkgRatio");
    hBkgRatio->Divide(hBkgNom);

    for(uint ibin = 0; ibin < shift_up_proj_ratio.size(); ibin++){
      shift_up_proj_ratio[ibin]->Divide(dataproj);
      shift_down_proj_ratio[ibin]->Divide(dataproj);
    }


    float chisquared = 0;

    for(int i = 1; i <= dataproj->GetXaxis()->GetNbins(); i++){
      float d_i = dataproj->GetBinContent(i);
      //  float m_i = hTotalMC_wgt->GetBinContent(i);
      float m_i = hTotalMC_wgt->GetBinContent(i);
      float error = dataproj->GetBinError(i)*dataproj->GetBinError(i) +
        hTotalMC_wgt->GetBinError(i)*hTotalMC_wgt->GetBinError(i);

      if(error == 0) continue;

      chisquared += (d_i - m_i) * (d_i - m_i)/error;

    }

    std::cout << "******ChiSq in bin :*****" << x << "\t" <<  chisquared << "\t" << "******And ChiSq/DoF is******* : " <<  chisquared/(dataproj->GetXaxis()->GetNbins() -3) << std::endl;
   

    std::stringstream low_a;    
    low_a<<std::fixed<<std::setprecision(2)<<
      chisquared/(dataproj->GetXaxis()->GetNbins() -3) ;
    std::stringstream high_a;
    high_a<<std::fixed<<std::setprecision(2)<<
      dataproj->GetXaxis()->GetNbins() - 3;
    //std::string caption = low.str() + " #leq E_{#nu} < " + high.str();

    std::string caption_a = "#chi^{2}/dof = " +
      low_a.str() + " (" +
      high_a.str() + ") ";

    hDataRatio->SetLineColor(kRed);
    hDataRatio_Wgt->SetLineColor(kBlue);

    dataproj->GetYaxis()->SetTitle("Events/8.09#times10^{20}POT");
    dataproj->SetMarkerStyle(20);
    //
    //TLegend *leg = new TLegend(0.3,0.4,0.89,0.89);
    TLegend *leg = new TLegend(0.4,0.65,0.89,0.89);

    leg->SetNColumns(2);
    
    leg->AddEntry(hMC[0], "Nominal MC", "PL");
    leg->AddEntry(hTotalMC_wgt, "Scaled MC", "PL");
    leg->AddEntry(hSigNom, "Nominal Sig", "PL");
    leg->AddEntry(hBkgNom, "Nominal Bkg", "PL");
    leg->AddEntry(hSigScl, "Scaled Sig", "PL");
    leg->AddEntry(hBkgScl, "Scaled Bkg", "PL");
    leg->AddEntry(dataproj, "Fake Data", "PL");


    TCanvas *cPID1 = new TCanvas("cPID1","cPID1",700,700);
    gStyle->SetPadBorderMode(0);
    gStyle->SetFrameBorderMode(0);
    Float_t small = 1e-5;
    cPID1->Divide(1,2,small,small);
    cPID1->cd(1);
    gPad->SetBottomMargin(small);
    CenterTitles(dataproj);

 
    //     dataproj->GetYaxis()->SetRangeUser(0,2000);

    dataproj->GetXaxis()->SetRangeUser(-1,0.92);
    //dataproj->GetYaxis()->SetRangeUser(0,1.3*dataproj->GetMaximum());
    dataproj->SetTitle(caption.c_str());
    dataproj->SetTitleOffset(0.08);
    dataproj->GetXaxis()->SetTitle("");
    dataproj->Draw();
    hMC[0]->Draw("hist same");
    hTotalMC_wgt->Draw("hist same");
    hSigNom->Draw("hist same");
    hBkgNom->Draw("hist same");
    hSigScl->Draw("hist same");
    hBkgScl->Draw("hist same");
    TGraphAsymmErrors* asymm =
      PlotWithSystErrorBand_mine(hMC[0],shift_up_proj,
                                 shift_down_proj,
                                 -1,-1,1.3,false);
  
    asymm->Draw("e2 same");
    dataproj->Draw("same");
    hMC[0]->Draw("hist same");
    hTotalMC_wgt->Draw("hist same");
    hSigNom->Draw("hist same");
    hBkgNom->Draw("hist same");
    hSigScl->Draw("hist same");
    hBkgScl->Draw("hist same");
    gPad->RedrawAxis();
    leg->Draw();
    // PrintChiSq(caption_a);
    Simulation();
    cPID1->cd(2);
    gPad->SetTopMargin(small);
    hDataRatio->GetYaxis()->SetTitle("MC/Fake Data");
    hDataRatio->GetXaxis()->SetRangeUser(-1,0.92);
    hDataRatio->GetYaxis()->SetRangeUser(0.1,1.9);
    hDataRatio->GetYaxis()->CenterTitle();
    hDataRatio->GetXaxis()->CenterTitle();
    hDataRatio->Draw("hist");
    TGraphAsymmErrors* asymm2 =
      PlotWithSystErrorBand_mine(hDataRatio,shift_up_proj_ratio,
                                 shift_down_proj_ratio,
                                 -1,-1,1.3,false);
    asymm2->Draw("e2 same");
    hDataRatio->Draw("hist same");
    hDataRatio_Wgt->Draw("hist same");
    gPad->RedrawAxis();
    char outname[50];
    sprintf(outname, "%s%s_bin_%i.pdf", "PlotsFinal/", "bdt", x);
    cPID1->SaveAs(outname);
  

    TLegend *legA = new TLegend(0.4,0.6,0.89,0.89);
    legA->AddEntry(hSignalRatio, "Sclaed Sig/Nom. Sig", "PL");
    legA->AddEntry(hBkgRatio, "Sclaed Bkg/Nom. Bkg", "PL");


  }//loop over analysis bins

  //Now For the Analysis Variable
  TH1F* data_y = (TH1F*)data->ProjectionY("data_y");

  data_y->SetMarkerStyle(20);
  data_y->GetXaxis()->SetRangeUser(0.2,1.5);

  std::vector<TH1F*> mc_y;
  std::vector<TH1F*> mc_y_wgt;
  for(uint i = 0; i < mc.size(); i++){
    char name[50];
    char name_wgt[50];
    sprintf(name, "%s_%i", "yproj", i);
    sprintf(name_wgt, "%s_%i", "yproj_wgt", i);
    mc_y.push_back((TH1F*)mc[i]->ProjectionY(name));
    mc_y_wgt.push_back((TH1F*)mc[i]->ProjectionY(name_wgt));   
  }

  std::vector<TH1F*> shift_up_y;
  std::vector<TH1F*> shift_down_y;
  std::vector<TH1F*> shift_up_y_ratio;
  std::vector<TH1F*> shift_down_y_ratio;
  for(uint i = 0; i < shifts_up.size(); i++){
    char name[50];
    char name_wgt[50];
    sprintf(name, "%s_%i", "syst", i);
    sprintf(name_wgt, "%s_%i", "syst_down", i);   
 
    shift_up_y.push_back((TH1F*)shifts_up[i]->ProjectionY(name));
    shift_down_y.push_back((TH1F*)shifts_up[i]->ProjectionY(name_wgt));

    sprintf(name, "%s_%i_ratio", "syst", i);
    sprintf(name_wgt, "%s_%i_ratio", "syst_down", i);  
  
    shift_up_y_ratio.push_back((TH1F*)shifts_up[i]->ProjectionY(name));
    shift_down_y_ratio.push_back((TH1F*)shifts_up[i]->ProjectionY(name_wgt));

  }

  mc_y_wgt[1]->Multiply(sig_weights);
  mc_y_wgt[6]->Multiply(cc_weights);
  mc_y_wgt[7]->Multiply(nc_weights);

  TH1F* hMCWeightedTot = (TH1F*)mc_y_wgt[1]->Clone("hMCWeightedTotal");
  hMCWeightedTot->Add(mc_y_wgt[6]);
  hMCWeightedTot->Add(mc_y_wgt[7]);
  hMCWeightedTot->SetLineColor(kBlue);
  
  TH1F* hDataRatio = (TH1F*)mc_y[0]->Clone("hDataRatio");
  TH1F* hDataRatio_Wgt = (TH1F*)hMCWeightedTot->Clone("hDataRatio_Wgt");
  
  hDataRatio->SetLineColor(kRed);
  hDataRatio_Wgt->SetLineColor(kBlue);
  hDataRatio->GetYaxis()->SetRangeUser(0.5,1.5);

  hDataRatio->Divide(data_y);
  hDataRatio_Wgt->Divide(data_y);

  for(uint i = 0; i < shifts_up.size(); i++){    
    shift_up_y_ratio[i]->Divide(data_y);
    shift_down_y_ratio[i]->Divide(data_y);
  }

  TLegend *leg = new TLegend(0.45,0.5,0.80,0.80);
  leg->AddEntry(data_y, "Fake Data", "PL");
  leg->AddEntry(mc_y[0], "Nominal MC", "PL");
  leg->AddEntry(hMCWeightedTot, "Scaled MC", "PL");

  float chisquared = 0;

  for(int i = data_y->GetXaxis()->FindBin(0.21);
      i <= data_y->GetXaxis()->FindBin(1.45); i++){
    float d_i = data_y->GetBinContent(i);
    float m_i = hMCWeightedTot->GetBinContent(i);
    float error = data_y->GetBinError(i)*data_y->GetBinError(i) +
      hMCWeightedTot->GetBinError(i)*hMCWeightedTot->GetBinError(i);

    if(error == 0) continue;

    chisquared += (d_i - m_i) * (d_i - m_i)/d_i;// error;
  }

  std::stringstream low_a;      
  low_a<<std::fixed<<std::setprecision(2)<<
    chisquared/(4) ;
  std::stringstream high_a;
  high_a<<std::fixed<<std::setprecision(2)<<
    4;
  //std::string caption = low.str() + " #leq E_{#nu} < " + high.str();

  std::string caption_a = "#chi^{2}/dof = " +
    low_a.str() + " (" +
    high_a.str() + ") ";

  std::cout << caption_a << std::endl;

  TCanvas *cPID2 = new TCanvas("cPID2","cPID2");
  cPID2->SetBottomMargin(0.15);
  gStyle->SetPadBorderMode(0);
  gStyle->SetFrameBorderMode(0);
  Float_t small = 1e-5;
  cPID2->Divide(1,2,small,small);
  cPID2->cd(1);
  // data_y->GetYaxis()->SetRangeUser(0,1.4*data_y->GetMaximum());
  data_y->GetYaxis()->SetTitle("Events/8.09#times10^{20}POT");
  data_y->GetYaxis()->SetRangeUser(0,50000);
  data_y->GetXaxis()->SetRangeUser(0.2,1.5);
  data_y->GetYaxis()->CenterTitle();
  data_y->GetXaxis()->SetTitle("");
  data_y->GetXaxis()->CenterTitle();
  data_y->Draw();
  TGraphAsymmErrors* asymm3 = 
    PlotWithSystErrorBand_mine(mc_y[0],shift_up_y,shift_down_y,
                               -1,-1,1.3,false);
  asymm3->Draw("e2 same");
  data_y->Draw("same");
  mc_y[0]->Draw("hist same");
  hMCWeightedTot->Draw("hist same");
  gPad->RedrawAxis();
  Simulation();
  leg->Draw();
  gPad->RedrawAxis();
  cPID2->cd(2);
  gPad->SetTopMargin(small);
  hDataRatio->GetYaxis()->SetTitle("Fake Data/MC");
  hDataRatio->GetYaxis()->SetRangeUser(0.4,1.6);
  hDataRatio->GetYaxis()->CenterTitle();
  hDataRatio->GetXaxis()->CenterTitle();  
  hDataRatio->Draw("hist");
  TGraphAsymmErrors* asymm4 = 
    PlotWithSystErrorBand_mine(hDataRatio,shift_up_y_ratio,
                               shift_down_y_ratio,
                               -1,-1,1.3,false);
  asymm4->Draw("e2 same");
  hDataRatio->Draw("hist same");
  hDataRatio_Wgt->Draw("hist same");
  gPad->RedrawAxis();
  cPID2->SaveAs("PlotsFinal/RecoKE.pdf");

}

void CompareToTrueSignal(TH2F* data, TH2F* data_sig, std::vector<TH2F*> mc,
                         TH1F* sig_weights, TH1F* cc_weights, 
                         TH1F* nc_weights,std::string dataname,
                         std::string weightname)
{
  TH1F* data_y_signal = (TH1F*)data_sig->ProjectionY("data_y_signal");

  std::vector<TH1F*> hMC_y;
  std::vector<TH1F*> hMC_y_wgt;
  for(uint i = 0; i < mc.size(); i++){
    char name[50];
    sprintf(name, "mc_%i", i);
    hMC_y.push_back((TH1F*)mc[i]->ProjectionY(name));
    char name_wgt[50];
    sprintf(name_wgt, "mc_wgt_%i", i);
    hMC_y_wgt.push_back((TH1F*)mc[i]->ProjectionY(name_wgt));
  
  }

  hMC_y_wgt[1]->Multiply(sig_weights);
  hMC_y_wgt[6]->Multiply(cc_weights);
  hMC_y_wgt[7]->Multiply(nc_weights);

  for(int i = 1; i < hMC_y_wgt[1]->GetXaxis()->GetNbins(); i++){
  
    float fSignal = hMC_y[1]->GetBinContent(i);
    float fSignalError = hMC_y[1]->GetBinError(i);
    float fTotal = hMC_y[0]->GetBinContent(i) - fSignal;
    float fWeight = sig_weights->GetBinContent(i);
    float fError  = sig_weights->GetBinError(i);

    float total_error = (fWeight * fSignal)* pow(pow((fSignalError/fSignal),2) + pow(fError/fWeight,2),0.5);

    hMC_y_wgt[1]->SetBinError(i,total_error); 
  }


  //Data True Signal-Like
  data_y_signal->SetLineColor(kOrange+7);
  data_y_signal->SetMarkerColor(kOrange+7);

  //MC Signal Template
  hMC_y[1]->SetLineColor(kRed);
  hMC_y[1]->SetMarkerColor(kRed); 

  //MC Signal Template
  hMC_y_wgt[1]->SetLineColor(kBlue);
  hMC_y_wgt[1]->SetMarkerColor(kBlue); 
  // hMC_y_wgt[1]->SetFillColor(kBlue);

  TLegend *leg = new TLegend(0.40,0.50,0.80,0.90);
  leg->AddEntry(data_y_signal, "Fake Data True Signal", "PL");
  leg->AddEntry(hMC_y[1], "Nominal Sig estimate", "PL");
  leg->AddEntry(hMC_y_wgt[1], "Scaled Sig estimate", "PL");  
  
  float chisquared = 0;
  float max_value = 0;

  for(int i = data_y_signal->GetXaxis()->FindBin(0.21);
      i <= data_y_signal->GetXaxis()->FindBin(1.45); i++){
    float d_i = data_y_signal->GetBinContent(i);
    float m_i = hMC_y_wgt[1]->GetBinContent(i);
    float e_d_i = data_y_signal->GetBinError(i);
    float e_m_i = hMC_y_wgt[1]->GetBinError(i);
    if(d_i > max_value) max_value = d_i;
    if(d_i == 0) continue;
    
    chisquared += (d_i - m_i) * (d_i - m_i)/d_i;// error;
  }
  
  std::stringstream low_a;      
  low_a<<std::fixed<<std::setprecision(2)<<
    chisquared/(4) ;
  std::stringstream high_a;
  high_a<<std::fixed<<std::setprecision(2)<<
    4;
  
  std::string caption_a = "#chi^{2}/dof = " +
    low_a.str() + " (" +
    high_a.str() + ") ";
  
  std::cout << caption_a << std::endl;
  
  TCanvas *c3 = new TCanvas("c3","c3");
  c3->SetBottomMargin(0.15);
  gStyle->SetPadBorderMode(0);
  gStyle->SetFrameBorderMode(0);
  Float_t small = 1e-5;
  c3->Divide(1,2,small,small);
  c3->cd(1);
  gPad->SetBottomMargin(small);
  data_y_signal->GetXaxis()->SetRangeUser(0.2,1.5);
  // data_y_signal->GetYaxis()->SetRangeUser(0, 50000);
  //data_y_signal->GetYaxis()->SetRangeUser(0, 2.0*max_value);
  data_y_signal->GetYaxis()->SetTitle("Signal events/8.09#times10^{20}POT");
  data_y_signal->GetYaxis()->CenterTitle();
  data_y_signal->SetMarkerStyle(20);
  data_y_signal->SetMarkerColor(kBlack);
  data_y_signal->SetLineColor(kBlack);
  data_y_signal->Draw();
  hMC_y_wgt[1]->Draw("hist same");
  hMC_y[1]->Draw("hist same");
  data_y_signal->Draw("same");
  gPad->RedrawAxis();
  Simulation();
  leg->Draw();
  // PrintChiSq2(caption_a);
  c3->cd(2);
  gPad->SetTopMargin(small);
  TH1F* hRatio_mc = (TH1F*)data_y_signal->Clone("hRatio_mc");
  TH1F* hRatio_wgt = (TH1F*)data_y_signal->Clone("hRatio_wgt");
  hRatio_mc->Divide(hMC_y[1]);
  hRatio_wgt->Divide(hMC_y_wgt[1]);
  hRatio_mc->SetLineColor(kRed);
  hRatio_wgt->SetLineColor(kBlue);
  hRatio_mc->GetYaxis()->SetRangeUser(0.8,1.2);
  hRatio_mc->GetYaxis()->SetTitle("Data/MC");
  hRatio_mc->GetYaxis()->CenterTitle();
  hRatio_mc->GetXaxis()->CenterTitle();
  hRatio_mc->Draw("hist");
  hRatio_wgt->Draw("hist same");
  char outname[50];
  sprintf(outname,"%s%s_%s_%s.pdf","PlotsFinal/", "signal_prediction",
          dataname.c_str(), weightname.c_str());
  c3->SaveAs("PlotsFinal/SignalPrediction.pdf");
  //  c3->Close();
  //delete c3;
}

void MakePlots()
{


  std::string outname_flux = 
    "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/PPFX_universes_bdte_energy_new_bdt_bins.root";
  TFile *flux_file = new TFile(outname_flux.c_str(), "read");

  std::string outname_genie = 
    "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/GENIE_universes_bdte_energy_new_bdt_bins.root";
  TFile *genie_file = new TFile(outname_genie.c_str(), "read");


  std::string outname_systs = 
    "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/Systs_samples_bdte_energy_new_bdt_bins.root";
  TFile *systs_file = new TFile(outname_systs.c_str(), "read");

  std::string outname_nominal = 
    "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/Nominal_MC_bdte_energy_new_bdt_bins.root";
  TFile *nominal_file = new TFile(outname_nominal.c_str(),"read");

  std::string outname_data =
    "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/Data_bdte_energy_new_bdt_bins.root";
  TFile *data_file = new TFile(outname_data.c_str(),"read");

 

  ////////////////////////////////////////////////////////////////////////
  //////////// PPFX Multiverse
  /*
    std::vector<std::unique_ptr<ana::Spectrum>> ppfx_spects;
    for(int i = 0; i <= 100; i++){
    if( i < 100){
    char name_ppfx[50];
    sprintf(name_ppfx, "%s_%i", "flux_universe", i);
    ppfx_spects.push_back(Spectrum::LoadFrom(flux_file->GetDirectory(name_ppfx)));
    }
    if( i == 100)
    ppfx_spects.push_back(Spectrum::LoadFrom(nominal_file->GetDirectory("nominal_mc")));
    }
  
    FluxMultiverseSyst fluxmulti(100,ppfx_spects,true);

    Spectrum* sPPFX_Nominal = fluxmulti.CV();
    Spectrum* sPPFX_Upper   = fluxmulti.UpperSigma();
    Spectrum* sPPFX_Lower   = fluxmulti.LowerSigma();

    TH2F* hPPFX_Nominal     = (TH2F*)sPPFX_Nominal->ToTH2(pot);
    TH2F* hPPFX_Upper       = (TH2F*)sPPFX_Upper->ToTH2(pot);
    TH2F* hPPFX_Lower       = (TH2F*)sPPFX_Lower->ToTH2(pot);
  
    TH1F* hPPFX_Upper_1d    = (TH1F*)sPPFX_Upper->ToTH1(pot);
    TH1F* hPPFX_Lower_1d    = (TH1F*)sPPFX_Lower->ToTH1(pot);
  
    std::vector<float> histvalue_up;
    std::vector<float> histvalue_low;
    for(int i = 1; i < hPPFX_Upper_1d->GetXaxis()->GetNbins(); i++){
    histvalue_up.push_back(hPPFX_Upper_1d->GetBinContent(i));
    histvalue_low.push_back(hPPFX_Lower_1d->GetBinContent(i));
    }

    double counter = 0;
    for(int i = 1; i <= hPPFX_Nominal->GetXaxis()->GetNbins(); i++){
    for(int j = 1; j <= hPPFX_Nominal->GetYaxis()->GetNbins(); j++){
    hPPFX_Upper->SetBinContent(i,j,histvalue_up[counter]);
    hPPFX_Lower->SetBinContent(i,j,histvalue_low[counter]);
    counter++;
    }
    } 
  */
  ////////////////////////////////////////////////////////////////////////
  //////////// Genie Multiverse
  std::vector<std::unique_ptr<ana::Spectrum>>genie_spects;
  for(int i = 0; i < 100; i++){
    char name[50];
    sprintf(name, "%s_%il", "genie_universe",i);
    genie_spects.push_back(Spectrum::LoadFrom(genie_file->GetDirectory(name)));
  }

  GenieMultiverseSpectra geniemulti = 
    GenieMultiverseSpectra(100,genie_spects,true);

  const Spectrum* sGenie_Nominal = geniemulti.Nominal();
  const Spectrum* sGenie_Upper   = geniemulti.UpperSigma();
  const Spectrum* sGenie_Lower   = geniemulti.LowerSigma();

  TH2F* hGenie_Nominal     = (TH2F*)sGenie_Nominal->ToTH2(pot);
  TH2F* hGenie_Upper       = (TH2F*)sGenie_Nominal->ToTH2(pot);
  TH2F* hGenie_Lower       = (TH2F*)sGenie_Nominal->ToTH2(pot);

  TH1F* hGenie_Upper_1d    = (TH1F*)sGenie_Upper->ToTH1(pot);
  TH1F* hGenie_Lower_1d    = (TH1F*)sGenie_Lower->ToTH1(pot);

  std::vector<float> histvalue_up_genie;
  std::vector<float> histvalue_low_genie;
  for(int i = 1; i < hGenie_Upper_1d->GetXaxis()->GetNbins(); i++){
    histvalue_up_genie.push_back(hGenie_Upper_1d->GetBinContent(i));
    histvalue_low_genie.push_back(hGenie_Lower_1d->GetBinContent(i));
  }

  double counter = 0;
  for(int i = 1; i <= hGenie_Nominal->GetXaxis()->GetNbins(); i++){
    for(int j = 1; j <= hGenie_Nominal->GetYaxis()->GetNbins(); j++){
      hGenie_Upper->SetBinContent(i,j,histvalue_up_genie[counter]);
      hGenie_Lower->SetBinContent(i,j,histvalue_low_genie[counter]);
      counter++;
    }
  }   

  ////////////////////////////////////////////////////////////////////////
  //////////// Other Systematics
  std::vector<ana::Spectrum> systs_spects;
  std::vector<ana::Spectrum> systs_spects_up;
  std::vector<ana::Spectrum> systs_spects_down;
 
  systs_spects_up.push_back(*Spectrum::LoadFrom(systs_file->GetDirectory("upcal_systs")));
  systs_spects_down.push_back(*Spectrum::LoadFrom(systs_file->GetDirectory("dwncal_systs")));
  systs_spects_up.push_back(*Spectrum::LoadFrom(systs_file->GetDirectory("lightup_systs")));
  systs_spects_down.push_back(*Spectrum::LoadFrom(systs_file->GetDirectory("lightdwn_systs")));

  systs_spects.push_back(*Spectrum::LoadFrom(systs_file->GetDirectory("calibshape_systs")));
  systs_spects.push_back(*Spectrum::LoadFrom(systs_file->GetDirectory("cherenkov_systs")));

  std::vector<TH2F*> systs_hists;
  std::vector<TH2F*> systs_hists_up;
  std::vector<TH2F*> systs_hists_down;

  for(uint i = 0; i < systs_spects.size(); i++){
    systs_hists.push_back((TH2F*)systs_spects[i].ToTH2(pot));
  }
  for(uint i = 0; i < systs_spects_up.size(); i++){
    systs_hists_up.push_back((TH2F*)systs_spects_up[i].ToTH2(pot));
    systs_hists_down.push_back((TH2F*)systs_spects_down[i].ToTH2(pot));
  }  

  systs_hists_up.push_back(hGenie_Upper);
  systs_hists_down.push_back(hGenie_Lower);  

  ////////////////////////////////////////////////////////////////////////
  //////////// Nominal MC
  std::vector<ana::Spectrum> nominal_spects;
  for(int i = 0; i < knumchns; i++){
    char name[50];
    sprintf(name, "%s_%s", "template", chns[i].name.c_str());
    nominal_spects.push_back(*Spectrum::LoadFrom(nominal_file->GetDirectory(name)));
  }

  std::vector<TH2F*> nominal_hists;
  for(uint i = 0; i < nominal_spects.size(); i++){
    nominal_hists.push_back((TH2F*)nominal_spects[i].ToTH2(pot));
  }
  
  ////////////////////////////////////////////////////////////////////////
  //////////// Data

  std::string dataname="data";
  std::string weightname = "nominal";
  //   std::string weightname = "maccres";

  char grabData[50];
  sprintf(grabData, "%s_TotMC", dataname.c_str()); // only for nominal one
  //sprintf(grabData, "%s_%s_TotMC", dataname.c_str(),weightname.c_str());   // fake data with sighigh and other
  
  Spectrum data_spect = *Spectrum::LoadFrom(data_file->GetDirectory(grabData));
  TH2F* data_hist = (TH2F*)data_spect.ToTH2(pot);

  char grabData_signal[50];
  sprintf(grabData_signal, "%s_Sig", dataname.c_str() );
  //  sprintf(grabData_signal, "%s_%s_Sig", dataname.c_str(),weightname.c_str());

  Spectrum data_signal_spect = *Spectrum::LoadFrom(data_file->GetDirectory(grabData_signal));
  TH2F* data_hist_signal = (TH2F*)data_signal_spect.ToTH2(pot);


  ////////////////////////////////////////////////////////////////////////
  //////////// Create Histograms to hold weights
  TH1F* sig_weights = (TH1F*)data_hist->ProjectionY("sig_weights");
  TH1F* cc_weights =  (TH1F*)data_hist->ProjectionY("cc_weights");
  TH1F* nc_weights =  (TH1F*)data_hist->ProjectionY("nc_weights");


  //////////////////////////////////////////////////////////////////////////
  ///////////////
  //Start Doing the Projections and Fitting
  ///////////////
  //////////////////////////////////////////////////////////////////////////
  for(int i = 1; i <= data_hist->GetYaxis()->GetNbins(); i++){

    //for(int i = data_hist->GetYaxis()->FindBin(0.0); i <= data_hist->GetYaxis()->FindBin(0.99); i++){
    std::stringstream low;      
    low<<std::fixed<<std::setprecision(2)<<
      data_hist->GetYaxis()->GetBinLowEdge(i);
    std::stringstream high;
    high<<std::fixed<<std::setprecision(2)<<
      data_hist->GetYaxis()->GetBinUpEdge(i);
    std::string caption = low.str() + "  #leq KE_{#pi^{0}} < " + high.str();


    /////////////////////////////////////////////////////////////////
    //Project Data and MC templates
    TH1F* hData1D = (TH1F*)data_hist->ProjectionX("hData1D",i,i);
    hData1D->Rebin(2);

    std::vector<TH1D*>hMC;
    for(int j = 0; j < (int)nominal_hists.size(); j++){
      char name[50];
      sprintf(name, "mc_bin_%i_sample_%i", i,j);
      hMC.push_back((TH1D*)nominal_hists[j]->ProjectionX(name,i,i));
      hMC[j]->Rebin(2);
    }

  

    /////////////////////////////////////////////////////////////////
    //Calculate Errors due to Each Systematic Sample
    std::vector<TH1D*>hSysts;
    for(int j = 0; j < (int)systs_hists.size(); j++){
      char name[50];
      sprintf(name, "syst_bin_%i_sample_%i", i,j);
      TH1D* Systs_holder = systs_hists[j]->ProjectionX(name,i,i);
      for(int bin = 1; bin <= Systs_holder->GetXaxis()->GetNbins(); bin++){
        double err = Systs_holder->GetBinContent(bin);
        double cv  = hMC[0]->GetBinContent(bin);
        Systs_holder->SetBinContent(bin,fabs(err-cv));
      }
      hSysts.push_back(Systs_holder);
      hSysts.back()->Rebin(2);
    }


    for(int j = 0; j < (int) systs_hists_up.size(); j++){
      char name[50];
      sprintf(name, "syst_up_bin_%i,smaple_%i", i, j);
      char name2[50];
      sprintf(name2, "syst_down_bin_%i,smaple_%i", i, j);
      TH1D* Systs_holder_up = systs_hists_up[j]->ProjectionX(name,i,i);
      TH1D* Systs_holder_down = systs_hists_down[j]->ProjectionX(name,i,i);

      for(int bin =1; bin <=Systs_holder_up->GetXaxis()->GetNbins(); bin++){
        double hi = Systs_holder_up->GetBinContent(bin);
        double lo = Systs_holder_down->GetBinContent(bin);
        double cv = hMC[0]->GetBinContent(bin);
        double value = 0.5 * (fabs(cv - hi) + fabs(cv - lo));
        Systs_holder_up->SetBinContent(bin, value);
      }

      hSysts.push_back(Systs_holder_up);
      hSysts.back()->Rebin(2);
    }

    /////////////////////////////////////////////////////////////////
    //Calculate Systematic Covariance Matrix Using hSysts and hMC[0]
    for(int j = 1; j <= hSysts[0]->GetXaxis()->GetNbins(); j++){
      for(int k = 1; k <= hSysts[0]->GetXaxis()->GetNbins(); k++){
        
        double multiple = 0;
        double x_j      = 0;
        double x_k      = 0;
        double counter  = 0;
        for(int iSyst = 0; iSyst < (int)hSysts.size(); iSyst++){
          double val_j = hMC[0]->GetBinContent(j) + 
            hSysts[iSyst]->GetBinContent(j); 
          double val_k = hMC[0]->GetBinContent(k) +
            hSysts[iSyst]->GetBinContent(k);
          multiple += val_j * val_k;
          x_j      += val_j;
          x_k      += val_k;
          counter++;
        }       

        double covariance = multiple/counter - ((x_j/counter)*(x_k/counter));

        covMatrix[j-1][k-1] = covariance;
        
      }
    }
    
 
    /////////////////////////////////////////////////////////////////
    //Create Matrix for Statistical Uncertainty
    for(int j = 1; j <= hSysts[0]->GetXaxis()->GetNbins(); j++){
      for(int k = 1; k <=hSysts[0]->GetXaxis()->GetNbins(); k++){
        if(j==k){
          stat_covMatrix[j-1][k-1] = 
            hData1D->GetBinError(j)*hData1D->GetBinError(j) + 
            hMC[0]->GetBinError(k)*hMC[0]->GetBinError(k);

        }
        else{
          stat_covMatrix[j-1][k-1] = 0.0;
        }
      }
    }

 

    /////////////////////////////////////////////////////////////////
    //Fill arrays for minuit fit
    for(int x = 1; x <= hData1D->GetXaxis()->GetNbins(); x++)
      {
        float data = hData1D->GetBinContent(x);
        float total  = hMC[0]->GetBinContent(x);
        float signal = hMC[1]->GetBinContent(x);
        float ccbkg  = hMC[6]->GetBinContent(x);
        float ncbkg  = hMC[7]->GetBinContent(x);
      
        fDA[x-1]      = data;
        fSig[x-1]     = signal;
        fCCBkg[x-1]    = ccbkg;
        fNCBkg[x-1]    = ncbkg;
          
      }

    TH2F* covariance_hist = new TH2F("covariance_hist",";NC#pi^{0}ID;NC#pi^{0}ID",
                                     xbins, 0 ,xbins, xbins, 0, xbins);
    TH2F* stat_covariance_hist = new TH2F("stat_covariance_hist",
                                          ";NC#pi^{0}ID;NC#pi^{0}ID",
                                          xbins, 0 ,xbins, xbins, 0, xbins);

    TH2F* correlation_hist = new TH2F("correlation_hist",
                                      ";NC#pi^{0}ID;NC#pi^{0}ID",
                                      xbins, 0 ,xbins, xbins, 0, xbins);

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

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

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


    //Covariance Matrix Inversion First
    TMatrixF Covij = TMatrix(xbins,xbins);
    for(int i =0; i < xbins; i++){
      for(int j = 0; j < xbins; j++){
        Covij[i][j] = covMatrix[i][j];
        covariance_hist->SetBinContent(i+1,j+1, covMatrix[i][j]);
        stat_covariance_hist->SetBinContent(i+1,j+1, stat_covMatrix[i][j]);

        correlation_hist->SetBinContent(i+1,j+1, 
                                        covMatrix[i][j]/(sqrt(covMatrix[i][i]) *
                                                         sqrt(covMatrix[j][j])));
      }
    } 

    stat_covariance_hist->Add(covariance_hist);

    TCanvas *c9 = new TCanvas("c9","c9");
    gStyle->SetTitleAlign(32);
    c9->SetRightMargin(0.15);
    covariance_hist->SetTitle(caption.c_str());
    covariance_hist->Draw("COLZ");
    Simulation();
    char cov_name[50];
    sprintf(cov_name, "%s%s_bin_%i.pdf", "PlotsFinal/", "error_matrix", i);
    c9->SaveAs(cov_name);
   

    TCanvas *cCor = new TCanvas("cCor","cCor");
    gStyle->SetTitleAlign(32);
    cCor->SetRightMargin(0.15);
    correlation_hist->GetZaxis()->SetTitle("Correlation");
    correlation_hist->SetTitle(caption.c_str());
    correlation_hist->Draw("COLZ");
    Simulation();
    char cor_name[50];
    sprintf(cor_name, "%s%s_bin_%i.pdf", "PlotsFinal/", "CorMatrix", i);
    cCor->SaveAs(cor_name);

  
    TMinuit *gMinuit = new TMinuit(3);
    gMinuit->SetPrintLevel(-1);
    gMinuit->SetFCN(fcn);
    Double_t arglist[4];
    Int_t ierflg = 0;
    arglist[0] = 1;
    gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
    
    // double vstart[3] = {1.2,0.8,0.9}; 
    double vstart[3] = {1.0,1.0,1.0}; 
    double verror[3] = {0,0,0};
    double vstep[3] = {0.1,0.1,0.1};

    gMinuit->mnparm(0, "Signal scaling", vstart[0], vstep[0], 0, 0, ierflg);
    gMinuit->mnparm(1, "CC-Bkg Scaling", vstart[1], vstep[1], 0, 0, ierflg);
    gMinuit->mnparm(2, "NC-Bkg Scaling", vstart[2], vstep[2], 0, 0, ierflg); 

    arglist[0] = 0; //number of iterations 0 == no limit
    arglist[1] = 0;

    gMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);

    double fParNewStart;
    double fParamNewErr;
    for(int ll = 0; ll < 3; ll++){
      gMinuit->GetParameter(ll,fParNewStart,fParamNewErr);
      vstart[ll] = fParNewStart;
    }


    gMinuit->mnexcm("MIGRAD", arglist, 3, ierflg);
    gMinuit->mnexcm("SEEk",  arglist, 3, ierflg);
    for(int ll = 0; ll < 3; ll++){
      gMinuit->GetParameter(ll,fParNewStart,fParamNewErr);
      vstart[ll] = fParNewStart;
      //vstep[ll] = 0.05;
    }


    TMinuit *hMinuit = new TMinuit(3);
    // hMinuit->SetPrintLevel(-1);
    hMinuit->SetFCN(fcn);
    arglist[0] = 1;
    hMinuit->mnexcm("SET ERR", arglist, 1, ierflg);

    hMinuit->mnparm(0, "Signal Scaling", vstart[0], vstep[0], 0, 0, ierflg);
    hMinuit->mnparm(1, "CC-Bkg Scaling", vstart[1], vstep[1], 0, 0, ierflg);
    hMinuit->mnparm(2, "NC-Bkg Scaling", vstart[2], vstep[2], 0, 0, ierflg); 


    arglist[0] = 2; //1 == standard minimization strategy 
                    //2 == try to improve minimium (slower)
     

    hMinuit->mnexcm("SET STR",arglist,1,ierflg);  
    arglist[0] = 0;
    hMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);
    hMinuit->mnexcm("MIGRAD", arglist, 5, ierflg);
    hMinuit->mnexcm("SEEk",  arglist, 4, ierflg);
    hMinuit->mnexcm("MINOS",  arglist, 4, ierflg);

    for(int ll = 0; ll < 3; ll++){
      hMinuit->GetParameter(ll,fParNewStart,fParamNewErr);
      vstart[ll] = fParNewStart;
      verror[ll] = fParamNewErr;
    }

    bool   troubleFitting = false;

    Double_t fmin = 0;
    Double_t fedm = 0;
    Double_t errdef = 0;
    Int_t       npari = 0;
    Int_t       nparx = 0;
    Int_t       istat = 0;
    hMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat);
    if(istat == 2)
      troubleFitting = true;
 
    if(!troubleFitting){
      hMinuit->Command("SHOw FCNvalue");
    }

    if(troubleFitting){
      //Fix Background Weights together
      TMinuit *jMinuit = new TMinuit(2);
      jMinuit->SetPrintLevel(-1);
      jMinuit->SetFCN(fcn2Var);
      //Double_t arglist[4];
      //Int_t ierflg = 0;
      arglist[0] = 1;
      jMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
    
      jMinuit->mnparm(0, "Signal Scaling", vstart[0], vstep[0], 0, 0, ierflg);
      jMinuit->mnparm(1, "TotBkg Scaling", vstart[1], vstep[1], 0, 0, ierflg);

      arglist[0] = 0; //number of iterations 0 == no limit
      arglist[1] = 0;
      
      jMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);
      jMinuit->mnexcm("MIGRAD", arglist, 3, ierflg);
      jMinuit->mnexcm("SEEk",  arglist, 3, ierflg);

      jMinuit->GetParameter(0,fParNewStart,fParamNewErr);
      vstart[0] = fParNewStart;
      verror[0] = fParamNewErr;
      // vstart[2] = fParNewStart;
      //verror[2] = fParamNewErr;
      jMinuit->GetParameter(1,fParNewStart,fParamNewErr);
      vstart[1] = fParNewStart;
      verror[1] = fParamNewErr; 
      vstart[2] = fParNewStart;
      verror[2] = fParamNewErr; 

      jMinuit->Command("SHOw FCNvalue");
    }


    //    int stat = gMinuit->GetStatus();
    //std::cout << "Error status = " << stat << std::endl;   

    double fParamVal;
    double fParamErr;
    hMinuit->GetParameter(0,fParamVal,fParamErr);
    double fParamVal_2;
    double fParamErr_2;
    hMinuit->GetParameter(1,fParamVal_2,fParamErr_2);
    double fParamVal_3;
    double fParamErr_3;
    hMinuit->GetParameter(2,fParamVal_3,fParamErr_3);

    std::cout << "*********Fit for bin " << i << ": " << std::endl;
    std::cout << "Signal Scaling Factor= " << vstart[0] << " +- " <<
      verror[0] << std::endl;
    sig_weights->SetBinContent(i, vstart[0]);
    sig_weights->SetBinError(i, verror[0]);

    std::cout << "CC-Bkg Scaling Factor= "  << vstart[1] << " +- " <<
      verror[1] << std::endl;
    cc_weights->SetBinContent(i, vstart[1]);
    cc_weights->SetBinError(i, verror[1]);
    std::cout << "NC-Bkg Factor= "  << vstart[2] << " +- " << 
      verror[2] << std::endl;
    nc_weights->SetBinContent(i, vstart[2]);
    nc_weights->SetBinError(i, verror[2]);
  
    delete gMinuit;  
    delete hMinuit;

  } //end of loop over analysis bins




  MakePlots(data_hist,nominal_hists,sig_weights,cc_weights,
            nc_weights, systs_hists_up,systs_hists_down,dataname,weightname);
  CompareToTrueSignal(data_hist,data_hist_signal, nominal_hists,
                      sig_weights, cc_weights, 
                      nc_weights,dataname,weightname);
 

  TCanvas *cWeights = new TCanvas("cWeights","cWeights");
  sig_weights->GetYaxis()->SetTitle("Signal Weight");
  sig_weights->GetXaxis()->SetTitleOffset(0.75);
  sig_weights->GetXaxis()->SetRangeUser(0.2,1.5);
  sig_weights->GetYaxis()->SetRangeUser(0,1.5);
  sig_weights->GetXaxis()->CenterTitle();
  sig_weights->GetYaxis()->CenterTitle();
  sig_weights->Draw("e1");
  Simulation();
 
  char out2[50];
  cWeights->SaveAs("PlotsFinal/SigWts.pdf");


  ///////////////////////////////////////////////////////////////////////////
  //CC-bkg Weights
  TCanvas *cWeights3 = new TCanvas("cWeights3","cWeights3");
  cc_weights->GetYaxis()->SetTitle("CCBkg Weight");
  cc_weights->GetXaxis()->SetRangeUser(0.2,1.5);
  cc_weights->GetYaxis()->SetRangeUser(0,1.5);
  cc_weights->GetXaxis()->SetTitleOffset(0.75);
  cc_weights->GetXaxis()->CenterTitle();
  cc_weights->GetYaxis()->CenterTitle();
  cc_weights->Draw("e1");
  Simulation();

  cWeights3->SaveAs("PlotsFinal/CCBkgWts.pdf");
  

  ///////////////////////////////////////////////////////////////////////////
  //NC-bkg Weights
  TCanvas *cWeights5 = new TCanvas("cWeights5","cWeights5");
  nc_weights->GetYaxis()->SetTitle("NC Weight");
  nc_weights->GetXaxis()->SetTitleOffset(0.75);
  nc_weights->GetXaxis()->CenterTitle();
  nc_weights->GetYaxis()->CenterTitle();
  nc_weights->GetXaxis()->SetRangeUser(0.2,1.5);
  nc_weights->GetYaxis()->SetRangeUser(0,1.5);
  nc_weights->Draw("e1");
  Simulation();
  cWeights5->SaveAs("PlotsFinal/NCWtd.pdf");
  
 

  // data_file->Close();
  // nominal_file->Close();
  //flux_file->Close();
  //genie_file->Close();
  //systs_file->Close();

}