makeEnergyEstimator.C

Go to the documentation of this file.

// I found that the energy estimator used by the nuebar oscillation analysis
// was not sufficient for this analysis due to training on a sample of 
// neutrinos with energy less than 4.5 GeV.
//
// This script marks the first piece of training a new estimator.
// Following docdb-26696, we need to normalize the true energy distribution
// to avoid biasing the energy estimation toward our peak beam energy.
// 
// Using the spectra obtained from getFlatWeightedFittingSpectra.C
// Fit the energy estimator

// dddoyle@colostate.edu

#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Analysis/Style.h"

#include "TH1D.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "TVectorD.h"
#include "TF2.h"
#include "TF1.h"

#include "NDAna/nuebarcc_inc/NuebarCCIncVars.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncCuts.h"
#include "NDAna/nuebarcc_inc/NuebarCCIncExtra.h"

#include <fstream>

using namespace ana;

void PrintPlot(TH1D* h, std::string name, std::string option);
void PrintPlot(TH2D* h, std::string name, std::string option);
void SetErrors(TH2D* hFlat_avg_energy, TH2* hFlat);
TH2D* MakeAverageTrueEnergySurface(TH2D* hFlat_energy, TH2D* hFlat);

Double_t PolN(Double_t *x, Double_t *par);

class PolyDef {
public:
  PolyDef(Int_t _nparam, TString _name, TF2* _func)
    : nparam(_nparam), name(_name), func(_func)
  {}
  PolyDef(TString _name, int _degree, bool _cross_terms)
    : nparam(9),
      name(_name),
      degree(_degree),
      cross_terms(_cross_terms)
  {
    func = new TF2(name, PolN, 0, 10, 0, 10, 9);
    FixTerms();
  }
  void FixTerms()
  {
    if(degree == 2) {
      func->FixParameter(4, 0);
      func->FixParameter(5, 0);
      func->FixParameter(7, 0);
      func->FixParameter(8, 0);
      if(!cross_terms) {
        func->FixParameter(6, 0);
      }
    }
    else {
      if(!cross_terms) {
        func->FixParameter(6, 0);
        func->FixParameter(7, 0);
        func->FixParameter(8, 0);
      }
    }
  }
  TF2* GetFunction() {return func;}
  TString GetName() {return name;}
  const Int_t nparam;
  TString name;
  TF2* func;
  bool cross_terms;
  int degree;
  
};

void PrintFitResults(TH2D* hist, std::ofstream& fout, PolyDef* pol);
void PrintErrorSurface(TH2D* hist, PolyDef* pol, std::string name);

void makeEnergyEstimator()
{
  std::string directory = "/nova/ana/users/ddoyle/NuebarEnergyEstimatorFiles/";
  TFile* fFlat = TFile::Open((directory + "flat_energy_spectra.root").c_str());

  Spectrum* sEM_vs_HAD_nom = Spectrum::LoadFrom(fFlat, "em_vs_had_nominal")).release(;
  Spectrum* sEM_vs_HAD_flat = Spectrum::LoadFrom(fFlat, "em_vs_had_flat")).release(;
  Spectrum* sEM_vs_HAD_flat_energy_sum = Spectrum::LoadFrom(fFlat, "em_vs_had_flat_energy_sum")).release(;

  Spectrum* sTrueE_nom = Spectrum::LoadFrom(fFlat, "trueE_nominal")).release(;
  Spectrum* sTrueE_flat = Spectrum::LoadFrom(fFlat, "trueE_flat")).release(;
  fFlat->Close();

  double POT = sTrueE_flat->POT();
  
  TH2D* hEM_vs_HAD_nom = (TH2D*) sEM_vs_HAD_nom->ToTH2(POT);
  TH2D* hEM_vs_HAD_flat = (TH2D*) sEM_vs_HAD_flat->ToTH2(POT);
  TH2D* hEM_vs_HAD_flat_energy_sum = (TH2D*) sEM_vs_HAD_flat_energy_sum->ToTH2(POT);
  TH1D* hTrueE_nom = (TH1D*)sTrueE_nom->ToTH1(POT);
  TH1D* hTrueE_flat = (TH1D*) sTrueE_flat->ToTH1(POT);  

  // Calculate average true nuetrino energy surface
  // Set errors on the average energy as avg_energy / sqrt(N) (error on the mean)
  TH2D* avg_energy_surface = MakeAverageTrueEnergySurface(hEM_vs_HAD_flat_energy_sum, hEM_vs_HAD_flat);

  // Make some plots
  PrintPlot(hEM_vs_HAD_nom, directory + "Plots/em_vs_had_nom.pdf", "colz");
  PrintPlot(hEM_vs_HAD_flat, directory + "Plots/em_vs_had_flat.pdf", "colz");
  PrintPlot(hEM_vs_HAD_flat_energy_sum, directory + "Plots/em_vs_had_flat_energy_sum.pdf", "colz");
  PrintPlot(avg_energy_surface, directory + "Plots/avg_energy_surface.pdf", "colz");
  PrintPlot(hTrueE_nom, directory + "Plots/trueE_nom.pdf", "hist");
  PrintPlot(hTrueE_flat, directory + "Plots/trueE_flat.pdf", "hist");
    
  // Do a fit to a polynomial
  // Try a few different function
  PolyDef* pol2  = new PolyDef("pol2", 2, false);
  PolyDef* pol2X = new PolyDef("pol2x", 2, true);
  PolyDef* pol3  = new PolyDef("pol3", 3, false);
  PolyDef* pol3X = new PolyDef("pol3x", 3, true);
 
  // Do the fits
  std::cout << "Starting fits...." << std::endl;
  std::cout << "-----------------------------------------\n";
  avg_energy_surface->Fit(pol2->name, "+");
  avg_energy_surface->Fit(pol2X->name, "+");
  avg_energy_surface->Fit(pol3->name, "+");
  avg_energy_surface->Fit(pol3X->name, "+");
  std::cout << "-----------------------------------------\n";
  
  // Print the fit parameters so they can be put into a header
  std::ofstream fout;
  fout.open("FitResults/all_fit_results.txt");
  PrintFitResults(avg_energy_surface, fout, pol2  );
  PrintFitResults(avg_energy_surface, fout, pol2X );
  PrintFitResults(avg_energy_surface, fout, pol3  );
  PrintFitResults(avg_energy_surface, fout, pol3X );
  fout.close();

  PrintErrorSurface(avg_energy_surface, pol2 , directory + "Plots/error_surface_pol2.pdf"  );
  PrintErrorSurface(avg_energy_surface, pol2X, directory + "Plots/error_surface_pol2x.pdf" );
  PrintErrorSurface(avg_energy_surface, pol3 , directory + "Plots/error_surface_pol3.pdf"  );
  PrintErrorSurface(avg_energy_surface, pol3X, directory + "Plots/error_surface_pol3x.pdf" );
}

void PrintErrorSurface(TH2D* hist, PolyDef* pol, std::string name)
{
  Double_t* pars = pol->GetFunction()->GetParameters();
  int dof = pol->GetFunction()->GetNDF();

  int NX = hist->GetNbinsX();
  int NY = hist->GetNbinsY();
  double x1 = hist->GetXaxis()->GetBinLowEdge(1);
  double x2 = hist->GetXaxis()->GetBinUpEdge(NX);
  double y1 = hist->GetYaxis()->GetBinLowEdge(1);
  double y2 = hist->GetYaxis()->GetBinUpEdge(NY);
  
  TH2D* ratio = new TH2D("ratio_" + pol->GetName(), 
                         "(Reco - True) /  True #bar{#nu} Energy " + pol->GetName(),
                         NX, x1, x2, NY, y1, y2);
  int nonempty = 0;
  for(int i = 1; i < NX + 1; i++) {
    for(int j = 1; j < NY + 1; j++) {
      if(hist->GetBinContent(i, j) > 0) {
        nonempty++;
        double x = hist->GetXaxis()->GetBinCenter(i);
        double y = hist->GetYaxis()->GetBinCenter(j);
        double pcoord[2] = {x, y};
        Double_t* coord = pcoord;       
        double fval = PolN(coord, pars);
        ratio->SetBinContent(i, j, (fval - hist->GetBinContent(i, j)) / hist->GetBinContent(i, j));
      }
    }
  }
  ratio->SetMaximum(0.5);
  ratio->SetMinimum(-0.5);
  ratio->GetXaxis()->SetTitle("EM Shower Energy [GeV]");
  ratio->GetYaxis()->SetTitle("Hadronic Energy [GeV]");
  
  SetPaletteBlueRedWhite();

  TDirectory* temp = gDirectory;
  TCanvas* c = new TCanvas();
  c->cd();
  ratio->Draw("colz");
  c->SaveAs(name.c_str());
  temp->cd();
  delete c;
}

void PrintFitResults(TH2D* hist, std::ofstream& fout, PolyDef* pol)
{
  TF1* fun = hist->GetFunction(pol->name);
  int dof = fun->GetNDF();
  // First display on screen
  std::cout << "-----------------------------------------\n";
  std::cout << "Fit Results: " << pol->name << std::endl;
  std::cout << "chi2 = " << fun->GetChisquare() << std::endl;
  std::cout << "chi2/dof = " << fun->GetChisquare()/dof << std::endl;
  for(int i = 0; i < pol->nparam; i++) {
    std::cout << "double p" << i << " = " << fun->GetParameter(i) << ";" << std::endl;
  }
  std::cout << "-----------------------------------------\n";
  
  // Now write it to file
  fout << "-----------------------------------------\n";
  fout << "Fit Results: " << pol->name << std::endl;
  fout << "chi2 = " << fun->GetChisquare() << std::endl;
  fout << "chi2/dof = " << fun->GetChisquare()/dof << std::endl;
  for(int i = 0; i < pol->nparam; i++) {
    fout << "double p" << i << " = " << fun->GetParameter(i) << ";" << std::endl;
  }
  fout << "-----------------------------------------\n";
  
  // now output each individual poly for parsing later
  std::ofstream polyout;
  polyout.open("FitResults/" + pol->name + "_fit_results.txt");
  for(int i = 0; i < pol->nparam; i++) {
    polyout << "p" << i << "=" << fun->GetParameter(i) << std::endl;
  }
  polyout << "chi2=" << fun->GetChisquare() << std::endl;
  polyout << "chi2/dof= " << fun->GetChisquare()/dof << std::endl;
  polyout.close();
  
}

Double_t PolN(Double_t *x, Double_t *par)
{
  Double_t result = par[0]*x[0] + par[1]*x[1]; // 1 order
  result += par[2]*x[0]*x[0] + par[3]*x[1]*x[1]; // 2 order 
  result += par[4]*x[0]*x[0]*x[0] + par[5]*x[1]*x[1]*x[1]; // 3 order
  result += par[6]*x[0]*x[1]; // 2 order cross terms
  result += par[7]*x[0]*x[0]*x[1] + par[8]*x[0]*x[1]*x[1]; // 3 order cross terms
  return result;
}

TH2D* MakeAverageTrueEnergySurface(TH2D* hFlat_energy, TH2D* hFlat)
{
  int NX = hFlat->GetNbinsX();
  int NY = hFlat->GetNbinsY();
  double x1 = hFlat->GetXaxis()->GetBinLowEdge(1);
  double x2 = hFlat->GetXaxis()->GetBinUpEdge(NX);
  double y1 = hFlat->GetYaxis()->GetBinLowEdge(1);
  double y2 = hFlat->GetYaxis()->GetBinUpEdge(NY);
  
  TH2D* ret = new TH2D("avg_energy", "Average True #bar{#nu} Energy [GeV]", NX, x1, x2, NY, y1, y2);
  ret->GetXaxis()->SetTitle("EM Shower Energy [GeV]");
  ret->GetYaxis()->SetTitle("Hadronic Energy [GeV]");
  for(int i = 1; i < NX + 1; i++) {
    for(int j = 1; j < NY + 1; j++) {
      double n_events = hFlat->GetBinContent(i, j);
      double avg_energy = hFlat_energy->GetBinContent(i, j) / n_events;
      double error = avg_energy / sqrt(n_events);
      if(avg_energy > 0) {
        ret->SetBinContent(i, j, avg_energy);
        ret->SetBinError(i, j, error);
      }
    }
  }  
  return ret;
}

void PrintPlot(TH1D* h, std::string name, std::string option)
{
  TDirectory* tmp = gDirectory;

  TCanvas* c = new TCanvas();
  c->cd();
  h->Draw(option.c_str());
  c->SaveAs(name.c_str());
  delete c;

  tmp->cd();
}

void PrintPlot(TH2D* h, std::string name, std::string option)
{
  TDirectory* tmp = gDirectory;

  TCanvas* c = new TCanvas();
  c->cd();
  h->Draw(option.c_str());
  c->SaveAs(name.c_str());
  delete c;

  tmp->cd();
}

void SetErrors(TH2D* hFlat_avg_energy, TH2* hFlat)
{
  int NX = hFlat->GetNbinsX();
  int NY = hFlat->GetNbinsY();
  for(int i = 1; i < NX + 1; i++) {
    for(int j = 1; j < NY + 1; j++) {
      double n_events = hFlat->GetBinContent(i, j);
      double avg_energy = hFlat_avg_energy->GetBinContent(i, j);
      double error = avg_energy / sqrt(n_events);
      hFlat_avg_energy->SetBinError(i, j, error);
    }
  }
}