NuebarCCIncTemplateFit.h

Go to the documentation of this file.

#pragma once

#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/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/XSec/Flux.h"
#include "CAFAna/XSec/TargetCount.h"
#include "CAFAna/Core/ReweightableSpectrum.h"
#include "CAFAna/Unfold/UnfoldIterative.h"
#include "CAFAna/Analysis/Plots.h"

//ROOT Includes
#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 "TFile.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TH1F.h"
#include "TH2F.h"
#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"
#include "TRandom.h"

//C++
#include <iostream>
#include <fstream>
#include <iomanip> 
#include <vector>
#include <cmath>
#include <math.h> 
#include <string>
#include <stdint.h>
#include <sstream>

using namespace ana;

const int pidBins = 21;
const bool shouldRebin = false;
double covMatrix[pidBins][pidBins];
double stat_covMatrix[pidBins][pidBins];
const float s_b_ratio = 0.5; //0.40; //ratio of signal/bkgd (in signal region)
                              // that needs to be exceeded 
                              // to perform fit in analysis bin
const float signalRegion = 0.85;

float fDA[pidBins];
float fSigLike[pidBins];
float fNumuLike[pidBins];
float fNC[pidBins];
float fOther[pidBins];

const std::vector<int> kNueCCColorDef = {
  632, //kRed, Total MC
  600,  //kBlue Signal
  616,  //kMagenta nuebar CC
  416-3, //kGreen - 3 NumuCC
  416+3, //kGreen + 3 Anti NumuCC
  800,   //kOrange NC
  400-3,  //kYellow - 3 Other
  416-3, //kGreen - 3 Total Background
};


std::string getCaption2D(TH3F* sample, int xbin, int ybin)
{
  std::stringstream low_X;
  low_X << std::fixed << std::setprecision(2) << 
    sample->GetXaxis()->GetBinLowEdge(xbin);
  std::stringstream hi_X;
  hi_X << std::fixed << std::setprecision(2) << 
    sample->GetXaxis()->GetBinUpEdge(xbin);

  std::stringstream low_Y;
  low_Y << std::fixed << std::setprecision(2) << 
    sample->GetYaxis()->GetBinLowEdge(ybin);
  std::stringstream hi_Y;
  hi_Y << std::fixed << std::setprecision(2) << 
    sample->GetYaxis()->GetBinUpEdge(ybin);

  std::string caption = low_X.str() + " #leq cos #theta_{e} < "+
    hi_X.str() + "   " + low_Y.str() + " #leq E_{e} < " + hi_Y.str();
  return caption;
}

void PrintCaption(TString str)
{

  TLatex* prelim = new TLatex(0.23, .80, str);
  prelim->SetTextColor(kBlack);
  prelim->SetNDC();
  prelim->SetTextSize(0.05);
  prelim->Draw();
}

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 numu    = 0;
  float nc      = 0;
  float other   = 0;

  //  Covariance Matrix Inversion First
  TMatrixF Vij_stat = TMatrix(pidBins,pidBins);
  TMatrixF Vij = TMatrix(pidBins,pidBins);
  for(int i =0; i < pidBins; i++){
    for(int j = 0; j < pidBins; 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,pidBins);
  TMatrixD FitVector2 = TMatrixD(pidBins,1);
  
  //Calculate Chi Sq
  for(int i =0; i < pidBins; i++){
    data    = fDA[i];
    signal  = fSigLike[i];
    numu    = fNumuLike[i];
    nc      = fNC[i];
    other   = fOther[i];  

    if(data < 30){
      FitVector[0][i] = 0;
      FitVector2[i][0] = 0;
      continue;
    }
    
    FitVector[0][i] = ( data - (par[1]*signal + par[0]*numu + 
                                par[2]*nc + other)); 
    
    FitVector2[i][0] = ( data -(par[1]*signal + par[0]*numu +
                                par[2]*nc +  other));       
  }
  
  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 numu    = 0;
  float nc      = 0;
  float other   = 0;

  //  Covariance Matrix Inversion First
  TMatrixF Vij_stat = TMatrix(pidBins,pidBins);
  TMatrixF Vij = TMatrix(pidBins,pidBins);
  for(int i =0; i < pidBins; i++){
    for(int j = 0; j < pidBins; 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();

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

  //Calculate Chi Sq
  for(int i =0; i < pidBins; i++){
    data    = fDA[i];
    signal  = fSigLike[i];
    numu    = fNumuLike[i];
    nc      = fNC[i];
    other   = fOther[i];  

    if(data < 30){
      FitVector[0][i] = 0;
      FitVector2[i][0] = 0;
      continue;

    }
    
    FitVector[0][i] = ( data - (par[1]*signal + par[0]*numu + 
                                par[0]*nc + other)); 
    
    FitVector2[i][0] = ( data -(par[1]*signal + par[0]*numu +
                                par[0]*nc + other));       
  }
  
  TMatrixD Chi = FitVector * VijInv;
  TMatrixD Chi2 = Chi * FitVector2;
  float secondterm = Chi2(0,0);
   
  f = secondterm; 
}





void PlotFitResults(TH1F* hData,
                    std::vector<TH1F*>hMCnom, std::vector<TH1F*>hMCwei,
                    TGraphAsymmErrors* g, std::string pidName,
                    std::string dataName, std::string caption,
                    int xbin,int ybin)
{

  TH1F* hMC_nom_signal = (TH1F*)hMCnom[1]->Clone(ana::UniqueName().c_str());
  hMC_nom_signal->Add(hMCnom[2],1);
  hMC_nom_signal->Add(hMCnom[7],1);
  TH1F* hMC_nom_bkgd = (TH1F*)hMCnom[3]->Clone(ana::UniqueName().c_str());
  hMC_nom_bkgd->Add(hMCnom[4],1);
  hMC_nom_bkgd->Add(hMCnom[5],1);
  hMC_nom_bkgd->Add(hMCnom[6],1);

  TH1F* hMC_wgt_signal = (TH1F*)hMCwei[1]->Clone(ana::UniqueName().c_str());
  hMC_wgt_signal->Add(hMCwei[2],1);
  hMC_wgt_signal->Add(hMCwei[7],1);
  TH1F* hMC_wgt_bkgd = (TH1F*)hMCwei[3]->Clone(ana::UniqueName().c_str());  
  hMC_wgt_bkgd->Add(hMCwei[2],1);
  hMC_wgt_bkgd->Add(hMCwei[4],1);
  hMC_wgt_bkgd->Add(hMCwei[5],1);
  hMC_wgt_bkgd->Add(hMCwei[6],1);

  hMC_nom_signal->SetLineColor(kBlue+3);
  hMC_nom_bkgd->SetLineColor(kGreen+3);

  TH1F* hMCWeight = (TH1F*)hMCwei[6]->Clone(ana::UniqueName().c_str());
  hMCWeight->Add(hMCwei[1]);
  hMCWeight->Add(hMCwei[2]);
  hMCWeight->Add(hMCwei[3]);
  hMCWeight->Add(hMCwei[4]);
  hMCWeight->Add(hMCwei[5]);
  hMCWeight->Add(hMCwei[7]);
  hMCWeight->SetLineColor(kRed);

  hMCnom[0]->SetLineColor(kMagenta-3);
  hData->SetMarkerStyle(20);
  hData->SetTitle(caption.c_str());
  hMCnom[0]->SetTitle(caption.c_str());

  if(hData->GetMaximum() > hMCnom[0]->GetMaximum())
    hMCnom[0]->GetYaxis()->SetRangeUser(0,hData->GetMaximum()*1.3);
  else(hMCnom[0]->GetYaxis()->SetRangeUser(0,hMCnom[0]->GetMaximum()*1.3));

  TCanvas* cCompare = new TCanvas(ana::UniqueName().c_str(),"cCompare");
  cCompare->SetLeftMargin(0.10);
  cCompare->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();
  cCompare->cd();
  pad1->Draw();
  pad1->cd();
  hMCWeight->SetTitle("");
  hMCnom[0]->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} POT");
  hMCnom[0]->GetXaxis()->SetLabelColor(kWhite);
  hMCnom[0]->GetYaxis()->SetLabelSize(0.03);
  hMCnom[0]->GetYaxis()->SetTitleSize(0.035);
  hMCnom[0]->SetTitle("");
  hMCnom[0]->Draw("hist");
  PrintCaption(caption.c_str());
  g->Draw("e2 same");
  hMCnom[0]->Draw("hist same");
  hMCWeight->Draw("hist same");
  hData->Draw("same");
  hMC_wgt_signal->Draw("hist same");
  hMC_nom_signal->Draw("hist same");
  hMC_wgt_bkgd->Draw("hist same");
  hMC_nom_bkgd->Draw("hist same");
  hData->Draw("same");
  TLegend* leg = ana::AutoPlaceLegend(0.3,0.3);
  leg->AddEntry(hData, "Fake Data", "p");
  leg->AddEntry(hMCWeight, "Scaled Total MC", "l");
  leg->AddEntry(hMC_wgt_signal, "Scaled Signal", "l");
  leg->AddEntry(hMC_wgt_bkgd,"Scaled Background", "l");
  leg->AddEntry(hMCnom[0],    "Total MC", "l");
  leg->AddEntry(hMC_nom_signal, "Signal", "l");
  leg->AddEntry(hMC_nom_bkgd,"Background", "l");
  leg->Draw();
  Simulation();
  gPad->RedrawAxis();
  cCompare->cd();
  pad2->Draw();
  pad2->cd();
  TH1F* hratio = (TH1F*)hMCWeight->Clone(ana::UniqueName().c_str());
  TH1F* hratio2 = (TH1F*)hMCnom[0]->Clone(ana::UniqueName().c_str());
  hratio->GetYaxis()->SetTitle("MC/Data");
  hratio->GetXaxis()->SetTitle("Electron Prong CVN Score");
  hratio->GetXaxis()->CenterTitle();
  hratio->SetMarkerColor(kBlack);
  hratio->SetLineColor(kBlack);
  hratio->SetMarkerStyle(20);
  hratio->Divide(hData);
  hratio2->Divide(hData);
  hratio->GetYaxis()->SetLabelSize(0.03);
  hratio->GetYaxis()->SetTitleSize(0.035); 
  hratio->GetYaxis()->SetRangeUser(0.4,1.5);
  hratio->Draw("same");
  hratio2->Draw("hist same");
  //hratio2->Draw("hist same");
  pad2->Draw("same");
  gPad->RedrawAxis();
  cCompare->cd();
  cCompare->Update();
  char outname[50];
  sprintf(outname, "%s%s_%s_postfit_ratio_bin_%i_%i.png", "Plots/", 
          pidName.c_str(),
          dataName.c_str(),xbin,ybin); 
  cCompare->SaveAs(outname);
  cCompare->Close();
}

TGraphAsymmErrors* PlotSystErrorBand(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();
  
  //  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);
        if(lo < 0 && hi > 0){ // if both are high, take largest deviation
          double max = std::max(abs(lo), abs(hi));
          lo = 0;
          hi = max;
        }
        if(lo > 0 && hi < 0) { // if both are low, take largest deviation
          double max = std::max(abs(lo), abs(hi));
          lo = max;
          hi = 0;
        }

        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);
  return g;
}


void CalculateCovarianceMatrix(TH1F* nom,std::vector<TH1F*> hSysts, 
                               int nPseudo,std::string pidName,
                               int xbin,int ybin)
{
  std::vector<TH1F*> vPseudoExp;
  
  TH2F* covariance_hist = 
    new TH2F(ana::UniqueName().c_str(),
             "Covariance Matrix;PID Bin;PID Bin;Covariance",
             nom->GetXaxis()->GetNbins(),
             0,nom->GetXaxis()->GetNbins(),
             nom->GetXaxis()->GetNbins(),
             0,nom->GetXaxis()->GetNbins());
                                  

  TRandom *r0 = new TRandom();
  for(int i = 0; i < nPseudo; i++){
    TH1F* hClone = (TH1F*)nom->Clone(ana::UniqueName().c_str());

    for(uint j = 0; j < hSysts.size(); j++){
      double wei = r0->Gaus(0,1);
      //std::cout << wei << std::endl;
      hClone->Add(hSysts[j],wei);
    }

    vPseudoExp.push_back(hClone);
  }

  //Mij = < (Val_i)*(Val_j)> - <Val_i> <Val_j>
  //covMatrix has to be filled
  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 i = 0; i < nPseudo; i++){
        double val_j = vPseudoExp[i]->GetBinContent(j); 
        double val_k = vPseudoExp[i]->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;
      covariance_hist->SetBinContent(j,k,covariance);
    }
  }  

  char out[50];
  sprintf(out, "Plots/covariance_%s_%i_%i.png", pidName.c_str(),xbin,ybin);
  TCanvas *c1 = new TCanvas("c1","c1");
  c1->SetLeftMargin(0.15);
  c1->SetRightMargin(0.15);
  covariance_hist->GetXaxis()->CenterTitle();
  covariance_hist->GetYaxis()->CenterTitle();
  covariance_hist->GetZaxis()->CenterTitle();
  covariance_hist->Draw("COLZ");
  c1->SaveAs(out);
  delete c1;
  return;
}

void CalculateTotalCovariance(TH1F* nom,std::vector<TH1F*> hSysts, 
                              std::string pidName,
                              int xbin,int ybin)
{
  std::vector<TH1F*> vPseudoExp;
  
  std::cout << "Hey Hey Right Here " << hSysts.size() << std::endl;

  TH2F* covariance_hist = 
    new TH2F(ana::UniqueName().c_str(),
             "Covariance Matrix;Systematic Number;Systematic Number;Covariance",
             (int)hSysts.size(),
             0,(int)hSysts.size(),
             (int)hSysts.size(),
             0,(int)hSysts.size());

  TH2F* correlation_hist = 
    new TH2F(ana::UniqueName().c_str(),
             "Correlation Matrix;;;Correlation",
             (int)hSysts.size(),
             0,(int)hSysts.size(),
             (int)hSysts.size(),
             0,(int)hSysts.size());

                                  
  std::vector<double> vAverages;
  for(int j = 0; j < (int)hSysts.size(); j++){
    double avg = 0;
    for(int i = 0; i < nom->GetXaxis()->GetNbins(); i++){
      avg += hSysts[j]->GetBinContent(i);
    }
    avg /= nom->GetXaxis()->GetNbins();
    vAverages.push_back(avg);
  }

  for(int i = 0; i < (int)vAverages.size(); i++){
    std::cout << vAverages[i] << std::endl;
  }

  const char *labels[6] = {"Cal.Shp.","Ckv","Calib","Light","XSec","PPFX"};


  //Sample Covariance = 1/N Sum( (x_i - <x>) (y_i - <y>))
  for(int j = 0; j < (int) hSysts.size(); j++){
    for(int k = 0; k < (int)hSysts.size(); k++){
      double summed_term = 0;
      double nterms = 0;
      for(int i = 0; i < nom->GetXaxis()->GetNbins(); i++){
        summed_term += (hSysts[j]->GetBinContent(i) - vAverages[j]) *
          (hSysts[k]->GetBinContent(i) - vAverages[k]);
        nterms++;
      } 
      double covariance = 1./nterms * summed_term;
      covariance_hist->SetBinContent(j+1,k+1,covariance);
      std::cout << j+1 << " " << k+1 << " " << covariance << std::endl;
    }
  } 

  for(int j = 1; j <= covariance_hist->GetXaxis()->GetNbins(); j++){
    for(int k = 1; k <= covariance_hist->GetXaxis()->GetNbins(); k++){
      double covariance = covariance_hist->GetBinContent(j,k);
      double correlation = 
        covariance/(sqrt(covariance_hist->GetBinContent(j,j)) *
                    sqrt(covariance_hist->GetBinContent(k,k)));
      correlation_hist->Fill(labels[j-1],labels[k-1],correlation);
    }
  } 


  char out[50];
  sprintf(out, "Plots/covariance_systs_%s_%i_%i.png", pidName.c_str(),xbin,ybin);
  TCanvas *c1 = new TCanvas("c1","c1");
  c1->SetLeftMargin(0.15);
  c1->SetRightMargin(0.15);
  covariance_hist->GetXaxis()->CenterTitle();
  covariance_hist->GetYaxis()->CenterTitle();
  covariance_hist->GetZaxis()->CenterTitle();
  covariance_hist->Draw("COLZ");
  c1->SaveAs(out);
  delete c1;


   correlation_hist->LabelsDeflate("X");
   correlation_hist->LabelsDeflate("Y");
   correlation_hist->LabelsOption("v");
  sprintf(out, "Plots/correlation_systs_%s_%i_%i.png", 
          pidName.c_str(),xbin,ybin);
  TCanvas *c2 = new TCanvas("c2","c2");
  c2->SetLeftMargin(0.15);
  c2->SetRightMargin(0.15);
  correlation_hist->GetXaxis()->CenterTitle();
  correlation_hist->GetYaxis()->CenterTitle();
  correlation_hist->GetZaxis()->CenterTitle();
  correlation_hist->Draw("COLZ text");
  c2->SaveAs(out);
  delete c2;
  return;

}

std::vector<TH3F*> BinByBinTemplateFit(TH3F* data,
                                       std::vector<TH3F*> templates,
                                       std::vector<TH3F*> systs_hists,
                                       std::vector<TH3F*> systs_hists_up,
                                       std::vector<TH3F*> systs_hists_down,
                                       std::vector<TH3F*> analysis_templates,
                                       std::string pidName, std::string varName,
                                       std::string dataName)
{
  std::vector<TH3F*> badentry;
  if(systs_hists_up.size() != systs_hists_down.size()){
    std::cerr << "Unequal up and down shifts. Exiting" << std::endl;
  }

  std::vector<TH3F*> hResults;
  for(uint i = 0; i < analysis_templates.size(); i++){
    hResults.push_back((TH3F*)analysis_templates[i]->
                      Clone(ana::UniqueName().c_str()));
  }


  for(int xbin = 0; xbin <= data->GetXaxis()->GetNbins(); xbin++){
    for(int ybin = 0; ybin <= data->GetYaxis()->GetNbins(); ybin++){
      
      std::string caption = getCaption2D(templates[0],xbin,ybin);
      

      //Variables to hold fit results
      std::pair<double,double> nue_wei  = std::make_pair (0,0);
      std::pair<double,double> numu_wei = std::make_pair (0,0);
      std::pair<double,double> nc_wei   = std::make_pair (0,0);
      Double_t covariance[3][3];

      //Project Data and MC to PID axis for template fit
      TH1F* hDataZ = (TH1F*)data->ProjectionZ(ana::UniqueName().c_str(),
                                              xbin,xbin,ybin,ybin);
      
      std::vector<TH1F*> hMCZ;
      for(uint numChns = 0; numChns < templates.size(); numChns++){
        hMCZ.push_back((TH1F*)templates[numChns]->
                       ProjectionZ(ana::UniqueName().c_str(),
                                   xbin,xbin,ybin,ybin));
        hMCZ[numChns]->SetLineColor(kNueCCColorDef[numChns]);
      }

      //Make Sure Fit should be performed in this bin (xbin,ybin)
      int pid_bin       = hDataZ->GetXaxis()->FindBin(signalRegion);
      float signal_like = hMCZ[1]->Integral(pid_bin,-1) +
        hMCZ[2]->Integral(pid_bin,-1) + hMCZ[7]->Integral(pid_bin,-1);
      float bkgd_amount = hMCZ[0]->Integral(pid_bin,-1) - signal_like;

      if(signal_like > 100 &&
         signal_like/bkgd_amount > s_b_ratio){
        
        //Calculate Uncertainties due to each systematic sample
        std::vector<TH1F*> hSystsZ;

        //For One-Sided Systematics
        for(int systNum = 0; systNum < (int) systs_hists.size(); systNum++){
          TH1F* hHolder = 
            (TH1F*)systs_hists[systNum]->ProjectionZ(ana::UniqueName().c_str(),
                                                     xbin,xbin,ybin,ybin);
          //Loop Over Bins, Calculate the absolute distance of systematic from
          //Nominal MC
          for(int bin = 1; bin <= hHolder->GetXaxis()->GetNbins(); bin++){
            double err = hHolder->GetBinContent(bin);
            double cv  = hMCZ[0]->GetBinContent(bin);
            hHolder->SetBinContent(bin,fabs(err-cv));
          }
          hSystsZ.push_back(hHolder);
        }

        //For Two-Side Systematics
        // Lets take a conservative approach, so furthest value away from
        // nominal is taken as the uncertainty for that bin. Should deal with
        // highly asymmetric systematics better
        for(int systNum = 0; systNum < (int)systs_hists_up.size(); systNum++){
          TH1F* hHolder     = (TH1F*)systs_hists_up[systNum]->
            ProjectionZ(ana::UniqueName().c_str(),xbin,xbin,ybin,ybin);
          TH1F* hHolder_dwn = (TH1F*)systs_hists_down[systNum]->
            ProjectionZ(ana::UniqueName().c_str(),xbin,xbin,ybin,ybin);

          //Loop over bins, convert two-sided uncertainty into a single value
          //which is the largest distance of each shift from nominal mc
          for(int bin = 1; bin <= hHolder->GetXaxis()->GetNbins(); bin++){
            double hi = hHolder->GetBinContent(bin);
            double lo = hHolder_dwn->GetBinContent(bin);
            double cv = hMCZ[0]->GetBinContent(bin);
            double val_hi = fabs(hi - cv);
            double val_lo = fabs(lo - cv);
            //      double value = 0.5*(fabs(cv-hi) + fabs(cv-lo));
            double value = (val_hi > val_lo)? val_hi : val_lo;
            hHolder->SetBinContent(bin,value);
          }
          hSystsZ.push_back(hHolder);
        }

        CalculateTotalCovariance(hMCZ[0],hSystsZ,pidName,xbin,ybin);
        //Calculate Systematic Covariance Matrix
        //Throw nPseudo experiments
        int nPseudo = 1000;
        CalculateCovarianceMatrix(hMCZ[0],hSystsZ,nPseudo,pidName,xbin,ybin);

        //Calculate Statistical Covariance Matrix
        //Diagonal with Data statistics for each pid bin
        for(int jj = 1; jj <= hSystsZ[0]->GetXaxis()->GetNbins(); jj++){
          for(int kk = 1; kk <= hSystsZ[0]->GetXaxis()->GetNbins(); kk++){
            if(jj==kk){
              stat_covMatrix[jj-1][kk-1] = 
                hDataZ->GetBinError(jj)*hDataZ->GetBinError(jj)
              +hMCZ[0]->GetBinError(jj)*hMCZ[0]->GetBinError(jj);
            }
            else{
              stat_covMatrix[jj-1][kk-1] = 0.0;
            }
          }
        }


        //Arrays are most(?) efficienct way to feed values to Minuit
        //Fill Arrays for MinuitFit
        for(int x = 1; x <= hDataZ->GetXaxis()->GetNbins(); x++){
          float data           = hDataZ->GetBinContent(x);
          //      float signal         = hMCZ[1]->GetBinContent(x);
          float nue         = hMCZ[1]->GetBinContent(x);
          float nuebar         = hMCZ[2]->GetBinContent(x);
          float numu           = hMCZ[3]->GetBinContent(x);
          float numubar        = hMCZ[4]->GetBinContent(x);
          float nc             = hMCZ[5]->GetBinContent(x);
          float other          = hMCZ[6]->GetBinContent(x);
          float nuenonfiducial = hMCZ[7]->GetBinContent(x);
          
          fDA[x-1]      = data;
          fSigLike[x-1]     = nue + nuebar + nuenonfiducial;
          fNumuLike[x-1]    = numu + numubar;
          fNC[x-1]      = nc;
          fOther[x-1]   = other;     
        }


        //Now perform the fit
        //Fit 1: Get Good Starting Fit Values, Quickly
        //Three parameter fit
        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] = {0.8,1.2,0.9};
        double verror[3] = {0,0,0};
        double vstep[3] = {0.01,0.01,0.01};

        //Define Fit Parameters
        gMinuit->mnparm(0, "Numu-Like Scaling", 
                        vstart[0], vstep[0], 0, 0, ierflg);
        gMinuit->mnparm(1, "Nue-Like Scaling", 
                        vstart[1], vstep[1], 0, 0, ierflg);
        gMinuit->mnparm(2, "NC Scaling", vstart[2], vstep[2], 0, 0, ierflg);
        
        arglist[0] = 0; //number of iterations 0 == no limit
        arglist[1] = 1; //1 == standard minimization strategy 
        
        //SIMPLEX should provide a good starting point for further fits
        gMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);

        //Fit 2: Improved Fit + Better Eror Estimation
        arglist[1] = 2; //2 == try to improve minimum (slower)
        gMinuit->mnexcm("SET STR", arglist,1,ierflg);

        gMinuit->mnexcm("SIMPLEX", arglist, 3, ierflg);
        gMinuit->mnexcm("MIGRAD", arglist, 5, ierflg);
        gMinuit->mnexcm("SEEk",  arglist, 4, ierflg);
        gMinuit->mnexcm("MINOS",  arglist, 4, ierflg);
        
        //Grab Final Results
        double fParNewStart;
        double fParNewErr;
        gMinuit->GetParameter(0,fParNewStart,fParNewErr);
        vstart[0] = fParNewStart;
        verror[0] = fParNewErr;
        gMinuit->GetParameter(1,fParNewStart,fParNewErr);
        vstart[1] = fParNewStart;
        verror[1] = fParNewErr;
        gMinuit->GetParameter(2,fParNewStart,fParNewErr);
        vstart[2] = fParNewStart;
        verror[2] = fParNewErr;

        //Before Ending Fits, Make Sure There Was No Issue With the Fit
        bool troubleFitting = false;
        Double_t fmin = 0,fedm = 0,ferrdef = 0;
        Int_t npari =0, nparx = 0, istat = 0;
        gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);

        Double_t corr_mat_test[3][3];
        gMinuit->mnemat(*corr_mat_test,3);

        //Turn this test error matix into a correlation matrix
        double corr_01 = corr_mat_test[0][1]/(corr_mat_test[0][0] *
                                              corr_mat_test[1][1]);
        double corr_02= corr_mat_test[0][2]/(corr_mat_test[0][0] *
                                             corr_mat_test[2][2]);
        double corr_12= corr_mat_test[1][2]/(corr_mat_test[1][1] *
                                             corr_mat_test[2][2]);
        if(istat == 2 || vstart[0] < 0 || vstart[1] < 0 || vstart[2] < 0 ||
           fabs(corr_01) > 0.96 || fabs(corr_02) > 0.96 || fabs(corr_12) > 0.96)
        troubleFitting=true;
        
        if(troubleFitting){
          //There Was trouble with the original fits,
          //most likely this was due to extremely high correlation (approx 1)
          //between two of the fitting parameters
          //Switch to a fit that adjusts all backgrounds
          TMinuit *jMinuit = new TMinuit(2);
          jMinuit->SetPrintLevel(-1);
          jMinuit->SetFCN(fcn2Var);
          arglist[0] = 1;
          jMinuit->mnexcm("SET ERR", arglist, 1, ierflg);

          jMinuit->mnparm(0, "Background", vstart[0], 
                          vstep[0], 0, 0, ierflg);
          jMinuit->mnparm(1, "Nue CC Scaling", vstart[1], 
                          vstep[1], 0, 0, ierflg);
        
        arglist[0] = 0; //number of iterations 0 == no limit
        arglist[1] = 1;
        
        jMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);
        jMinuit->mnexcm("MIGRAD", arglist, 3, ierflg);
        jMinuit->mnexcm("SEEk",  arglist, 3, ierflg);
        jMinuit->mnexcm("MINOS", arglist, 3, ierflg);


        jMinuit->GetParameter(0,fParNewStart,fParNewErr);
        vstart[0] = fParNewStart;
        verror[0] = fParNewErr;
        jMinuit->GetParameter(0,fParNewStart,fParNewErr);         
        vstart[2] = fParNewStart;
        verror[2] = fParNewErr;
        jMinuit->GetParameter(1,fParNewStart,fParNewErr);
        vstart[1] = fParNewStart;
        verror[1] = fParNewErr;
        //Get the covariance matrix
        Double_t emat[2][2];
        jMinuit->mnemat(*emat,2);

        nue_wei  = std::make_pair (vstart[1],verror[1]);
        numu_wei  = std::make_pair (vstart[0],verror[0]);
        nc_wei  = std::make_pair (vstart[2],verror[2]);
        covariance[0][0] = emat[0][0];
        covariance[0][1] = emat[0][1];
        covariance[0][2] = 0;
        covariance[1][0] = emat[1][0];
        covariance[1][1] = emat[1][1];
        covariance[1][2] = 0;   
        covariance[2][0] = emat[0][0];
        covariance[2][1] = emat[0][1];
        covariance[2][2] = emat[0][0];
        }
        else{
        gMinuit->GetParameter(0,fParNewStart,fParNewErr);
        vstart[0] = fParNewStart;
        verror[0] = fParNewErr;
        gMinuit->GetParameter(2,fParNewStart,fParNewErr);         
        vstart[2] = fParNewStart;
        verror[2] = fParNewErr;
        gMinuit->GetParameter(1,fParNewStart,fParNewErr);
        vstart[1] = fParNewStart;
        verror[1] = fParNewErr;

        //Get the covariance matrix
        Double_t emat[3][3];
        gMinuit->mnemat(*emat,3);

        nue_wei  = std::make_pair (vstart[1],verror[1]);
        numu_wei  = std::make_pair (vstart[0],verror[0]);
        nc_wei  = std::make_pair (vstart[2],verror[2]);
        covariance[0][0] = emat[0][0];
        covariance[0][1] = emat[0][1];
        covariance[0][2] = emat[0][2];
        covariance[1][0] = emat[1][0];
        covariance[1][1] = emat[1][1];
        covariance[1][2] = emat[1][2];  
        covariance[2][0] = emat[2][0];
        covariance[2][1] = emat[2][1];
        covariance[2][2] = emat[2][2];
        }
      }
      else{
        std::cout << "******** Skipping Fit For Bin " 
                  << caption << "\t\n\n\n\n" << std::endl;
        nue_wei = std::make_pair (1,0);
        numu_wei = std::make_pair (1,0);
        nc_wei   = std::make_pair (1,0);
        covariance[0][0] = 0;
        covariance[0][1] = 0;
        covariance[0][2] = 0;
        covariance[1][0] = 0;
        covariance[1][1] = 0;
        covariance[1][2] = 0;       
        covariance[2][0] = 0;
        covariance[2][1] = 0;
        covariance[2][2] = 0;       
      }
      
      std::cout << "*********Fit for bin " << caption << " : " << std::endl;
      std::cout << "NumuCC Scaling Factor= " << numu_wei.first << " +- " <<
        numu_wei.second << std::endl;      
      std::cout << "NueCC Scaling Factor= "  << nue_wei.first << " +- " <<
        nue_wei.second << std::endl;
      std::cout << "NC Scaling Factor= "  << nc_wei.first << " +- " << 
        nc_wei.second << std::endl;


      //Projections For Systematic Error Band Drawing
      std::vector<TH1F*> hSysts;
      for(int z = 0; z < (int)systs_hists.size(); z++){
        hSysts.push_back((TH1F*)systs_hists[z]->
                         ProjectionZ(ana::UniqueName().c_str(),
                                     xbin,xbin,ybin,ybin));
      }
      std::vector<TH1F*> hSystsUp;
      for(int z = 0; z < (int)systs_hists_up.size();z++){
        hSystsUp.push_back((TH1F*)systs_hists_up[z]->
                         ProjectionZ(ana::UniqueName().c_str(),
                                     xbin,xbin,
                                     ybin,ybin));       
      }
      std::vector<TH1F*> hSystsDown;
      for(int z = 0; z < (int)systs_hists_down.size();z++){
        hSystsDown.push_back((TH1F*)systs_hists_down[z]->
                         ProjectionZ(ana::UniqueName().c_str(),
                                     xbin,xbin,
                                     ybin,ybin));       
      }

      std::vector<TH1F*> shifts_up;
      std::vector<TH1F*> shifts_down;

      for(int z = 0; z < (int)hSysts.size();z++){
        shifts_up.push_back(hSysts[z]);
        shifts_down.push_back((TH1F*)hMCZ[0]->Clone(ana::UniqueName().c_str()));
      }

      for(int z = 0; z < (int)hSystsUp.size(); z++){
        shifts_up.push_back(hSystsUp[z]);
        shifts_down.push_back(hSystsDown[z]);
      }

      //Get Systematic Error Band
      TGraphAsymmErrors* g2 = PlotSystErrorBand(hMCZ[0],
                                                shifts_up,
                                                shifts_down,
                                                -1,-1,1.3,false);

      //Make Copies of Nominal MC template predictions
      std::vector<TH1F*> hMCnom;
      for(int z = 0; z < (int)hMCZ.size(); z++){
        hMCnom.push_back((TH1F*)hMCZ[z]->Clone(ana::UniqueName().c_str()));
      }
      
      //Scale MC templates by fitted weights
      hMCZ[1]->Scale(nue_wei.first);  //Nue
      hMCZ[2]->Scale(nue_wei.first);  //Nuebar
      hMCZ[3]->Scale(numu_wei.first); //Numu
      hMCZ[4]->Scale(numu_wei.first); //Numubar
      hMCZ[5]->Scale(nc_wei.first);   //NC
      hMCZ[6]->Scale(1);              //Other
      hMCZ[7]->Scale(nue_wei.first);  //Non-fiducial Nuebar


      //Plot Template Before after Fit
      PlotFitResults(hDataZ,hMCnom,hMCZ,g2, pidName,dataName,caption,xbin,ybin);
          
      //Apply Fit Results to Analysis Histogram
      for(int numChn = 0; numChn < (int)hResults.size(); numChn++){
        for(int zbin = 0; zbin < hResults[0]->GetZaxis()->GetNbins();zbin++){
          float nnue     = hResults[1]->GetBinContent(xbin,ybin,zbin);
          float nnuebar  = hResults[2]->GetBinContent(xbin,ybin,zbin);
          float nnumu    = hResults[3]->GetBinContent(xbin,ybin,zbin);
          float nnumubar = hResults[4]->GetBinContent(xbin,ybin,zbin);
          float nnc      = hResults[5]->GetBinContent(xbin,ybin,zbin);
          float nother   = hResults[6]->GetBinContent(xbin,ybin,zbin);
          float nnf      = hResults[7]->GetBinContent(xbin,ybin,zbin);
          if(numChn == 0){
            float binCont = nnue * nue_wei.first;
            binCont += nnuebar * nue_wei.first;
            binCont += nnumu * numu_wei.first;
            binCont += nnumubar * numu_wei.first;
            binCont += nnc * nc_wei.first;
            binCont += nother;
            binCont += nnf * nue_wei.first;
            hResults[0]->SetBinContent(xbin,ybin,zbin,binCont);

            float binErr2 =
              pow(covariance[0][0],2)*pow((nnumu+nnumubar),2) +
              pow(covariance[1][1],2)* pow((nnue+nnuebar+nnf),2) +
              pow(covariance[2][2],2)* pow((nnc),2) +
              2 * covariance[0][1] * (nnumu+nnumubar)*(nnue+nnuebar+nnf) +
              2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) +
              2 * covariance[1][2] * (nnue + nnuebar + nnf) * (nnc) +
              nnue * pow(nue_wei.first,2) + nnuebar * pow(nue_wei.first,2) +
              nnf * pow(nue_wei.first,2) + nnumu * pow(numu_wei.first,2) +
              nnumubar * pow(numu_wei.first,2) + nnc * pow(nc_wei.first,2) +
              nother;

            if(covariance[0][1] == 0 && covariance[0][2] == 0 &&
               covariance[1][2] == 0 && nue_wei.second == 0 &&
               numu_wei.second == 0 && nc_wei.second == 0){
              binErr2 = nnue + nnuebar + nnf + nnumu + nnumubar + 
                nnc + nother;
            }


            //if(binErr2 < 0) binErr2 = fabs(binErr2);

            hResults[0]->SetBinError(xbin,ybin,zbin,sqrt(binErr2));

            std::cout << zbin << "\t" << binCont << "\t" << 
              sqrt(binErr2) << std::endl;
            if(std::isnan(sqrt((binErr2))) || binErr2 > 1e6){
              std::cout << pow(covariance[0][0],2)*pow((nnumu+nnumubar),2) <<
                "\t" << pow(covariance[1][1],2)* pow((nnue+nnuebar+nnf),2) <<
                "\t" << pow(covariance[2][2],2)* pow((nnc),2) << "\t" <<
                 2 * covariance[0][1] * (nnumu+nnumubar)*(nnue+nnuebar+nnf) << 
                "\t" << 2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) <<
                "\t" << 2 * covariance[1][2] * (nnue + nnuebar + nnf) * (nnc) <<
                "\t" << nnue * pow(nue_wei.first,2) << "\t" <<
                nnuebar * pow(nue_wei.first,2) << "\t" <<
                nnf * pow(nue_wei.first,2) << "\t" <<
                nnumu * pow(numu_wei.first,2) << "\t" << 
                nnumubar * pow(numu_wei.first,2) << "\t" << 
                nnc * pow(nc_wei.first,2) << "\t" << nother << std::endl;
            }
          }
          if(numChn == 1){
            float binCont = nnue * nue_wei.first;
            float binErr  = sqrt( pow(nue_wei.second,2)*pow(nnue,2) + 
                                  nnue * pow(nue_wei.first,2));
            hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
            hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
          }
          if(numChn == 2){
            float binCont = nnuebar * nue_wei.first;
            float binErr  = sqrt( pow(nue_wei.second,2)*pow(nnuebar,2) +
                                  nnuebar * pow(nue_wei.first,2) +
                                  pow(nue_wei.second,2)*pow(nnue,2) + // Add uncertainty in ws to signal component
                                  nnue * pow(nue_wei.first,2));       //
            hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
            hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
          }
          if(numChn == 3){
            float binCont = nnumu * numu_wei.first;
            float binErr  = sqrt( pow(numu_wei.second,2)*pow(nnumu,2) +
                                  nnumu * pow(numu_wei.first,2));
            hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
            hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
          }
          if(numChn == 4){
            float binCont = nnumubar * numu_wei.first;
            float binErr  = sqrt( pow(numu_wei.second,2)*pow(nnumubar,2) +
                                  nnumubar * pow(numu_wei.first,2));
            hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
            hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
          }
          if(numChn == 5){
            float binCont = nnc * nc_wei.first;
            float binErr  = sqrt( pow(nc_wei.second,2)*pow(nnc,2) +
                                  nnc * pow(nc_wei.first,2));
            hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
            hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
          }
          if(numChn == 7){
            float binCont = nnf * nue_wei.first;
            float binErr  = sqrt( pow(nue_wei.second,2)*pow(nnf,2) +
                                  nnf * pow(nue_wei.first,2));
            hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
            hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
          }
        }//loop over zbins
      } // loop over interaction channels
    }//loop over ybins
  }//loop over xbins

  return hResults;
}