estimate_energy.C

Go to the documentation of this file.

//---------------------------------------------------------
// Neutrino Energy Estimator for Inclusive numu CC Analysis
// 
// This script is mainly focused to generate 2D histogram  
// which are used in the estimation of muon and hadronic 
// energy.
//
// This script need to be run thrice.
// 1st : Calculate for Muon Energy
// -- Reco Muon Track Length in Active  vs True Muon Energy
// -- Reco Muon Track Length in Active  vs True Muon E - True Muon E in Catcher
// -- Reco Muon Track Length in Catcher vs True Muon Energy in Catcher
//
// 2nd : After Reconstruction of Muon Energy
// -- visible hadronic energy vs True Nu E - Reco Muon Energy
//
// 3rd : Neutrino Energy
// -- Sum of the reco muon energy and hadronic energy
//
//
// Author : Biswaranjan Behera
// email  : bbehera@fnal.gov
//
//---------------------------------------------------------

#include "CAFAna/Core/Loaders.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Var.h"
#include "CAFAna/Core/Cut.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/TruthCuts.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "CAFAna/Vars/PPFXWeights.h"
#include "3FlavorAna/Cuts/NumuCuts.h"

#include "CAFAna/Cuts/Cuts.h"
#include "CAFAna/Cuts/TruthCuts.h"

#include "NDAna/numubarcc_inc/NumubarCCIncVars.h"
#include "NDAna/numubarcc_inc/NumubarCCIncCuts.h"
#include "NDAna/numubarcc_inc/NumubarCCIncBins.h"
#include "NDAna/numubarcc_inc/NumubarCCIncSystDefs.h"

#include "TCanvas.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TStyle.h"

using namespace ana;
// using namespace ana::xsec::numubarcc;

void estimate_energy(string outfileName = "numubar_energy_estimate.root")
{
  // ND respin Monte Carlo files
  const std::string nd_mc = ana::xsec::numubarcc::nominal_dataset;
  
  SpectrumLoader * loader = new SpectrumLoader(nd_mc);
  loader->SetSpillCut(kStandardSpillCuts);
  
  struct Plot2D
  {
    std::string name;
    std::string labelx;
    Binning binsx;
    Var varx;
    std::string labely;
    Binning binsy;
    Var vary;
    Cut cut;
  };
  
  const Var kTrueMuEAct = ana::xsec::numubarcc::kTrueMuE - ana::xsec::numubarcc::kTrueCatcherE;
  // const Var kHadE       = kTrueNuE - kIncXsecMuonE;
 
  const Cut kSelectCut        = ana::xsec::numubarcc::kIsNueorbarCC && 
                          ana::xsec::numubarcc::kTrueMuon && 
                          ana::xsec::numubarcc::kPreselect && 
                          ana::xsec::numubarcc::kMuonIDCut;
  
  const int kNumPlots2D = 3;
  
  Plot2D plots2D[kNumPlots2D] = {

    //--------------------------------------------------------------------------------- 
    // Calculate for Muon Energy
    // 1st : Reco Muon Track Length            vs True Muon Energy
    // 2nd : Reco Muon Track Length            vs True Muon E - True Muon E in Catcher
    // 3rd : Reco Muon Track Length in Catcher vs True Muon Energy in Catcher
    //--------------------------------------------------------------------------------- 
    
    {"MuonE_hist_act","Reco Muon Track Length (m)",
     ana::xsec::numubarcc::kTrkLenActBinning, ana::xsec::numubarcc::kTrkLenAct, "True Muon Energy (GeV)" ,
     ana::xsec::numubarcc::kMuonEActBinning, ana::xsec::numubarcc::kTrueMuE, kSelectCut && ana::xsec::numubarcc::kActive},

    {"MuonE_hist_actAndCat","Reco Muon Track Length (m)",
     ana::xsec::numubarcc::kTrkLenActBinning, ana::xsec::numubarcc::kTrkLenAct, "True Muon E - True Muon E in Catcher (GeV)" ,
     ana::xsec::numubarcc::kMuonEActBinning, kTrueMuEAct, kSelectCut && ana::xsec::numubarcc::kActiveAndCatcher},

    {"MuonE_hist_cat","Reco Muon Track Length in Catcher (m)",
     ana::xsec::numubarcc::kTrkLenCatBinning, ana::xsec::numubarcc::kTrkLenCat, "True Muon Energy in Catcher (GeV)" ,
     ana::xsec::numubarcc::kMuonECatBinning, ana::xsec::numubarcc::kTrueCatcherE, kSelectCut && ana::xsec::numubarcc::kActiveAndCatcher},

    // {"HadE","Visible Hadronic Energy (GeV)",
    //  kVisHadEBinning, kVisHadE, "True Nu E - Reco Muon Energy (GeV)" ,
    //  kMuonEActBinning, kHadE, kCut},
  };
  
  Spectrum* sSpect2D[kNumPlots2D];
  
  for(int i = 0; i < kNumPlots2D; i++)
    {
      Plot2D p = plots2D[i];
      sSpect2D[i] = new Spectrum(p.labelx, p.labely, *loader, p.binsx, p.varx, p.binsy, p.vary, p.cut, kNoShift, ana::xsec::numubarcc::std_wgt);
    } // loop over 2D plots
    
  loader->Go();
    
  TFile * fout = TFile::Open(outfileName.c_str(),"RECREATE");
  
  for(int i = 0; i < kNumPlots2D; i++)
    {
      Plot2D p = plots2D[i];
      sSpect2D[i]->SaveTo(fout, p.name);
    } // loop over 2D plots

} //end of void