NueCCIncCrossSectionFunctions.h

Go to the documentation of this file.

#pragma once

#include "NDAna/nuecc_inc/NueCCIncTemplateFit.h"
#include "RooUnfoldResponse.h"
#include "RooUnfoldBayes.h"

void PrintChiSq(TString str)
{
  TLatex* prelim = new TLatex(.22, .85, str);
  prelim->SetTextColor(kBlack);
  prelim->SetNDC();
  prelim->SetTextSize(0.05);
  prelim->Draw();
}

void Plot2D(TH2F* hist,std::string outName, bool setLog = false)//char* outName)
{
  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str(),"c1");
  c1->SetLeftMargin(0.20);
  c1->SetRightMargin(0.20);
  if(setLog)c1->SetLogz();
  hist->GetZaxis()->SetTitleOffset(1.10);
  hist->GetXaxis()->CenterTitle();
  hist->GetYaxis()->CenterTitle();
  hist->GetZaxis()->CenterTitle();
  hist->Draw("COLZ");
  char out[50];
  sprintf(out, "Plots/%s", outName.c_str());
  c1->SaveAs(out);
  return;
}
/*
std::vector<TH3F*> ApplyFitResults(std::vector<TH3F*> mc, 
                                   std::vector<TH2F*> wei)
{
  std::vector<TH3F*> results;
  
  for(uint i = 0; i < mc.size(); i++){
    results.push_back((TH3F*)mc[i]->Clone(ana::UniqueName().c_str()));
  }

  if(wei.size() != 4){
  std::cout << "Something failed" << std::endl;
  return results;
  }

  for(int xbin = 0; xbin <= mc[0]->GetXaxis()->GetNbins(); xbin++){
    for(int ybin = 0; ybin <= mc[0]->GetYaxis()->GetNbins(); ybin++){
      
      float nuewei = wei[0]->GetBinContent(xbin,ybin);
      float nueerr = wei[0]->GetBinError(xbin,ybin);
      float numuwei = wei[1]->GetBinContent(xbin,ybin);
      float numuerr = wei[1]->GetBinError(xbin,ybin);
      float ncwei = wei[2]->GetBinContent(xbin,ybin);
      float ncerr = wei[2]->GetBinError(xbin,ybin);
      float correl = wei[3]->GetBinContent(xbin,ybin);

      for(int zbin = 0; zbin <= mc[0]->GetZaxis()->GetNbins(); zbin++){

        float n_sig = results[1]->GetBinContent(xbin,ybin,zbin);
        float n_nuebar = results[2]->GetBinContent(xbin,ybin,zbin);
        float n_numu = results[3]->GetBinContent(xbin,ybin,zbin);
        float n_numubar = results[4]->GetBinContent(xbin,ybin,zbin);
        float n_nc = results[5]->GetBinContent(xbin,ybin,zbin);
        float n_other = results[6]->GetBinContent(xbin,ybin,zbin);
        float n_nuenonfid = results[7]->GetBinContent(xbin,ybin,zbin);

        float err_sig = results[1]->GetBinError(xbin,ybin,zbin);
        float err_nuebar = results[2]->GetBinError(xbin,ybin,zbin);
        float err_numu = results[3]->GetBinError(xbin,ybin,zbin);
        float err_numubar = results[4]->GetBinError(xbin,ybin,zbin);
        float err_nc = results[5]->GetBinError(xbin,ybin,zbin);
        float err_other = results[6]->GetBinError(xbin,ybin,zbin);
        float err_nuenonfid = results[7]->GetBinError(xbin,ybin,zbin);



        //Propagate Only Systematic Uncertainties
        float val = nuewei*(n_sig+n_nuebar+n_nuenonfid)+
          numuwei*(n_numu+n_numubar)+ncwei*(n_nc)+n_other;
        float err_da = sqrt( pow(n_sig + n_nuebar + n_nuenonfid,2)*
                             pow(nueerr,2) +
                          pow(n_nc+n_numu+n_numubar,2)*pow(ncerr,2) +
                          2*(n_sig+n_nuebar+n_nuenonfid)*
                          (n_nc+n_numu+n_numubar)*correl*nueerr*ncerr);
        results[0]->SetBinContent(xbin,ybin,zbin,val);
        results[0]->SetBinError(xbin,ybin,zbin,err_da);

        float err_background = sqrt(pow(n_nuebar+n_nuenonfid,2)*pow(nueerr,2) +
                                    pow(n_nc+n_numu+n_numubar,2)*pow(ncerr,2) +
                                    2*(n_nuebar+n_nuenonfid)*
                                    (n_nc+n_numu+n_numubar)*
                                    correl*nueerr*ncerr);


        //err = sqrt(pow(nueerr,2)*pow(n_sig,2) + n_sig*pow(nuewei,2));
        //float err = sqrt(pow(err_da,2)*pow(err_background,2));
        float err = sqrt(pow(nueerr,2)*pow(n_sig,2));
        results[1]->SetBinContent(xbin,ybin,zbin,nuewei*n_sig);
        results[1]->SetBinError(xbin,ybin,zbin,err);

        //err = sqrt(pow(nueerr,2)*pow(n_nuebar,2) + n_nuebar*pow(nuewei,2));
        err = sqrt(pow(nueerr,2)*pow(n_nuebar,2));
        results[2]->SetBinContent(xbin,ybin,zbin,nuewei*n_nuebar);
        results[2]->SetBinError(xbin,ybin,zbin,err);

        //err = sqrt(pow(numuerr,2)*pow(n_numu,2) + n_numu*pow(numuwei,2));
        err = sqrt(pow(numuerr,2)*pow(n_numu,2));
        results[3]->SetBinContent(xbin,ybin,zbin,numuwei*n_numu);
        results[3]->SetBinError(xbin,ybin,zbin,err);

        //err = sqrt(pow(numuerr,2)*pow(n_numubar,2) 
        //           + n_numubar*pow(numuwei,2));
        err = sqrt(pow(numuerr,2)*pow(n_numubar,2));
        results[4]->SetBinContent(xbin,ybin,zbin,numuwei*n_numubar);
        results[4]->SetBinError(xbin,ybin,zbin,err);

        //err = sqrt(pow(ncerr,2)*pow(n_nc,2) + n_nc*pow(ncwei,2));
        err = sqrt(pow(ncerr,2)*pow(n_nc,2));
        results[5]->SetBinContent(xbin,ybin,zbin,ncwei*n_nc);
        results[5]->SetBinError(xbin,ybin,zbin,err);

        //err = sqrt(pow(nueerr,2)*pow(n_nuenonfid,2) + 
        //         n_nuenonfid*pow(nuewei,2));
        err = sqrt(pow(nueerr,2)*pow(n_nuenonfid,2));
        results[7]->SetBinContent(xbin,ybin,zbin,nuewei*n_nuenonfid);
        results[7]->SetBinError(xbin,ybin,zbin,err);
      }
    }
  }


  //results[0]->Add(mc[0],-1);

  //for(uint i = 1; i < results.size(); i++) 
  //results[0]->Add(results[i],1);
  
  return results;
}
*/

TH3F* GetEfficiency(TH3F* num, TH3F* denom)
{
  TH3F* eff = (TH3F*)num->Clone("eff");
  eff->Divide(denom);

  return eff;
}

TH1F* DoUnfolding(TH1F* hMeasured, TH2F* hResponse, int iter)
{
  /*
  //Pretty Colors
  //-------------------------------------------------------------------------
  const Int_t NRGBs = 9;
  const Int_t NCont = 255;
  
  Double_t stops[NRGBs] = { 0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};
  Double_t red[9]   = { 251./255., 235./255., 223./255., 214./255., 196./255., 188./255., 157./255., 102./255.,  37./255.};
  Double_t green[9] = { 251./255., 185./255., 132./255., 91./255., 67./255., 37./255., 25./255., 29./255., 37./255.};
  Double_t blue[9]  = { 251./255., 187./255., 137./255., 98./255., 66./255., 45./255., 33./255., 32./255., 37./255.};
  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
  gStyle->SetNumberContours(NCont);
*/


  RooUnfoldResponse ResponseMatrix(0,0,hResponse, "ResponseMatrix");
  std::cout << "Using Overflow? : " << 
    ResponseMatrix.UseOverflowStatus() << std::endl;
  if(!ResponseMatrix.UseOverflowStatus()){
    ResponseMatrix.UseOverflow(1);
    std::cout << "Using Overflow, now? : " << 
      ResponseMatrix.UseOverflowStatus() << std::endl;
  }

  RooUnfoldBayes unfold(&ResponseMatrix, hMeasured, iter);
  TH1F* hunfolded = (TH1F*)unfold.Hreco();
  //TMatrixD covMatrix = unfold.Ereco();

  //Draw Response Matrix
  TCanvas *c1 = new TCanvas("c1","c1");
  hResponse->SetMarkerColor(kBlack);
  hResponse->Draw("COLZ TEXT");
  gStyle->SetPaintTextFormat("4.f");
  Simulation();
  c1->SaveAs("Plots/response_matrix.png");
  c1->Close();
  delete c1;

  /*
  //Draw covarivance matrix
  TCanvas *c2 = new TCanvas("c2","c2");
  covMatrix.Draw("COLZ TEXT");
  gStyle->SetPaintTextFormat("4.f");
  Simulation();
  c2->SaveAs("Plots/covariance_matrix.png");
  c2->Close();
  delete c2;
  */

  //Draw Unfolded Histogram
  TCanvas *c3 = new TCanvas("c3","c3");
  hunfolded->SetTitle("");
  hunfolded->Draw();
  Simulation();
  c3->SaveAs("Plots/unfolded_signal_estimate.png");
  c3->Close();
  delete c3;

  return hunfolded;
}

void CalculateXSec(TH3F* hNumerator,TH1F* hFlux, TH3F* hEff,
                   std::string dataName,std::string xsecResult, 
                   double numTargets)
{
  TFile* outf = new TFile("CrossSection_MCResults.root","update");
  
  if(xsecResult == "energy"){
    TH3F* hXSec = (TH3F*)hNumerator->Clone("hXSec");
    
    //3D Efficiency Correction
    hXSec->Divide(hEff);
    
    TH1F* hXSec1D = (TH1F*)hXSec->ProjectionZ("hXSec1D");
    TH1F* hTrueVar = (TH1F*)hXSec1D->Clone("clone");
    hFlux->Scale(1e-4);//nu/cm^2 from nu/m^2
    hXSec1D->Divide(hFlux);
    hXSec1D->Scale(1./numTargets);
    hXSec1D->SetTitle("");
    hXSec1D->GetXaxis()->SetTitle("Neutrino Energy, E_{#nu} (GeV)");
    hXSec1D->GetYaxis()->SetTitle("#sigma (cm^{2}/nucleon)");
    hXSec1D->GetXaxis()->CenterTitle();
    hXSec1D->GetYaxis()->CenterTitle();

    char name[50];
    sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
            "xsec");
    hXSec1D->SetName(name);
    hXSec1D->Write();
    sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
            "truevar_eff_corrected");
    hTrueVar->SetName(name);
    hTrueVar->Write();
    sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
            "efficiency"); 
    hEff->SetName(name);
    hEff->Write();

    outf->Close();    
  }
  if(xsecResult == "electron"){
    TH3F* hXSec = (TH3F*)hNumerator->Clone("hXSec");
    //3D Efficiency Correction
    hXSec->Divide(hEff);
    
    TH1F* hXSec1D = (TH1F*)hXSec->ProjectionY("hXSec1D");
    TH1F* hTrueVar = (TH1F*)hXSec1D->Clone("clone");
    hFlux->Scale(1e-4);//nu/cm^2 from nu/m^2
    double fluxtot = hFlux->Integral();
    hXSec1D->Scale(1./fluxtot);
    hXSec1D->Scale(1./numTargets);
    hXSec1D->Scale(1,"width");
    hXSec1D->SetTitle("");
    hXSec1D->GetXaxis()->SetTitle("Electron Energy, E_{e} (GeV)");
    hXSec1D->GetYaxis()->SetTitle("#sigma (cm^{2}/nucleon GeV)");
    hXSec1D->GetXaxis()->CenterTitle();
    hXSec1D->GetYaxis()->CenterTitle();

    char name[50];
    sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
            "xsec");
    hXSec1D->SetName(name);

    hXSec1D->Write();
    sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
            "truevar_eff_corrected");
    hTrueVar->SetName(name);
    hTrueVar->Write();
    sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
            "efficiency"); 
    hEff->SetName(name);
    hEff->Write();

    outf->Close();    
  }
  if(xsecResult == "cos"){
    TH3F* hXSec = (TH3F*)hNumerator->Clone("hXSec");
    
    //3D Efficiency Correction
    hXSec->Divide(hEff);
    
    TH1F* hXSec1D = (TH1F*)hXSec->ProjectionX("hXSec1D");
    TH1F* hTrueVar = (TH1F*)hXSec1D->Clone("clone");

    hFlux->Scale(1e-4);//nu/cm^2 from nu/m^2
    double fluxtot = hFlux->Integral();
    hXSec1D->Scale(1./fluxtot);
    hXSec1D->Scale(1./numTargets);
    hXSec1D->Scale(1,"width");
    hXSec1D->SetTitle("");
    hXSec1D->GetXaxis()->SetTitle("cos #theta_{e}");
    hXSec1D->GetYaxis()->SetTitle("#sigma (cm^{2}/nucleon cos #theta)");
    hXSec1D->GetXaxis()->CenterTitle();
    hXSec1D->GetYaxis()->CenterTitle();

    char name[50];
    sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
            "xsec");
    hXSec1D->SetName(name);

    hXSec1D->Write();
    sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
            "truevar_eff_corrected");
    hTrueVar->SetName(name);
    hTrueVar->Write();
    sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
            "efficiency"); 
    hEff->SetName(name);
    hEff->Write();

    outf->Close();    
  }
}