estimate_energy.C
Go to the documentation of this file.
1 //---------------------------------------------------------
2 // Neutrino Energy Estimator for Inclusive numu CC Analysis
3 //
4 // This script is mainly focused to generate 2D histogram
5 // which are used in the estimation of muon and hadronic
6 // energy.
7 //
8 // This script need to be run thrice.
9 // 1st : Calculate for Muon Energy
10 // -- Reco Muon Track Length in Active vs True Muon Energy
11 // -- Reco Muon Track Length in Active vs True Muon E - True Muon E in Catcher
12 // -- Reco Muon Track Length in Catcher vs True Muon Energy in Catcher
13 //
14 // 2nd : After Reconstruction of Muon Energy
15 // -- visible hadronic energy vs True Nu E - Reco Muon Energy
16 //
17 // 3rd : Neutrino Energy
18 // -- Sum of the reco muon energy and hadronic energy
19 //
20 //
21 // Author : Biswaranjan Behera
22 // email : bbehera@fnal.gov
23 //
24 //---------------------------------------------------------
25 
26 #include "CAFAna/Core/Loaders.h"
27 #include "CAFAna/Core/Spectrum.h"
29 #include "CAFAna/Core/Var.h"
30 #include "CAFAna/Core/Cut.h"
31 #include "CAFAna/Cuts/SpillCuts.h"
32 #include "CAFAna/Cuts/TruthCuts.h"
35 
40 
41 #include "TCanvas.h"
42 #include "TFile.h"
43 #include "TH1.h"
44 #include "TH2.h"
45 #include "TStyle.h"
46 
47 using namespace ana;
48 // using namespace ana::xsec::numubarcc;
49 
50 void estimate_energy(string outfileName = "numubar_energy_estimate.root")
51 {
52  // ND respin Monte Carlo files
53  const std::string nd_mc = "nd_rhc_prod5_decaf_minus_muonid_training_minus_fakedata";
54 
55  SpectrumLoader * loader = new SpectrumLoader(nd_mc);
57 
58  struct Plot2D
59  {
61  std::string labelx;
62  Binning binsx;
63  Var varx;
64  std::string labely;
65  Binning binsy;
66  Var vary;
67  Cut cut;
68  };
69 
70  // const Var kTrueMuEAct = ana::xsec::numubarcc::kTrueMuE - ana::xsec::numubarcc::kTrueCatcherE;
71  // const Var kHadE = kTrueNuE - kIncXsecMuonE;
72 
79 
80  const int kNumPlots2D = 3;
81 
82  Plot2D plots2D[kNumPlots2D] = {
83 
84  //---------------------------------------------------------------------------------
85  // Calculate for Muon Energy
86  // 1st : Reco Muon Track Length vs True Muon Energy
87  // 2nd : Reco Muon Track Length vs True Muon E - True Muon E in Catcher
88  // 3rd : Reco Muon Track Length in Catcher vs True Muon Energy in Catcher
89  //---------------------------------------------------------------------------------
90 
91  {"MuonE_hist_act","Reco Muon Track Length (m)",
94 
95  {"MuonE_hist_actAndCat","Reco Muon Track Length (m)",
96  ana::xsec::numubarcc::kTrkLenActBinning, ana::xsec::numubarcc::kTrkLenAct, "True Muon E - True Muon E in Catcher (GeV)" ,
98 
99  {"MuonE_hist_cat","Reco Muon Track Length in Catcher (m)",
100  ana::xsec::numubarcc::kTrkLenCatBinning, ana::xsec::numubarcc::kTrkLenCat, "True Muon Energy in Catcher (GeV)" ,
101  ana::xsec::numubarcc::kMuonECatBinning, ana::xsec::numubarcc::kTrueCatcherE, kSelectCut && ana::xsec::numubarcc::kActiveAndCatcher},
102 
103  // {"HadE","Visible Hadronic Energy (GeV)",
104  // kVisHadEBinning, kVisHadE, "True Nu E - Reco Muon Energy (GeV)" ,
105  // kMuonEActBinning, kHadE, kCut},
106  };
107 
108  Spectrum* sSpect2D[kNumPlots2D];
109 
110  for(int i = 0; i < kNumPlots2D; i++)
111  {
112  Plot2D p = plots2D[i];
113  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);
114  } // loop over 2D plots
115 
116  loader->Go();
117 
118  TFile * fout = TFile::Open(outfileName.c_str(),"RECREATE");
119 
120  for(int i = 0; i < kNumPlots2D; i++)
121  {
122  Plot2D p = plots2D[i];
123  sSpect2D[i]->SaveTo(fout, p.name);
124  } // loop over 2D plots
125 } //end of void
126 
127 
const Cut kRecoVtxNumuFiducialCut([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;return VtxInBounds(&sr->vtx.elastic.vtx, numucc_fiducial_min, numucc_fiducial_max);})
const XML_Char * name
Definition: expat.h:151
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Cut kActive([](const caf::SRProxy *sr){int ibesttrk=kBestMuonIDIndex(sr);if(sr->trk.kalman.ntracks< 1||ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return false;const caf::SRKalmanTrackProxy &bestmuon=sr->trk.kalman.tracks[ibesttrk];return(bestmuon.leninact > 0 && bestmuon.lenincat< 0);})
const char * p
Definition: xmltok.h:285
std::vector< double > Spectrum
Definition: Constants.h:746
void Plot2D(TH1 *total, std::map< std::string, TH1 * > systs, std::vector< std::string > to_plot, std::string basename)
void SetSpillCut(const SpillCut &cut)
const Cut kActiveAndCatcher([](const caf::SRProxy *sr){int ibesttrk=kBestMuonIDIndex(sr);if(sr->trk.kalman.ntracks< 1||ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return false;const caf::SRKalmanTrackProxy &bestmuon=sr->trk.kalman.tracks[ibesttrk];return(bestmuon.leninact > 0 && bestmuon.lenincat > 0);})
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
const Binning kTrkLenCatBinning
const Binning kMuonECatBinning
virtual void Go() override
Load all the registered spectra.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:534
const Var kTrkLenCat([](const caf::SRProxy *sr){int ibesttrk=kBestMuonIDIndex(sr);if(sr->trk.kalman.ntracks< 1||ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f;const caf::SRKalmanTrackProxy &bestmuon=sr->trk.kalman.tracks[ibesttrk]; if(bestmuon.leninact > 0 && bestmuon.lenincat<=0) return 0.f; if(bestmuon.leninact > 0 && bestmuon.lenincat > 0) return float(bestmuon.lenincat/100.); if(bestmuon.leninact<=0 && bestmuon.lenincat > 0) return float((bestmuon.leninact/100.) +(bestmuon.lenincat/100.));return-1000.f;})
loader
Definition: demo0.py:10
void estimate_energy(string outfileName="numubar_energy_estimate.root")
const SystShifts kNoShift
Definition: SystShifts.cxx:22
const Binning kMuonEActBinning
const ana::Var std_wgt
const Cut kTrueKalmanMuon([](const caf::SRProxy *sr){ int ibesttrk=kBestMuonIDIndex(sr);if(ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return false;return std::abs(sr->trk.kalman.tracks[ibesttrk].truth.pdg)==13;})
const Cut cut
Definition: exporter_fd.C:30
const Cut kSelectCut
const Var kTrkLenAct([](const caf::SRProxy *sr){int ibesttrk=kBestMuonIDIndex(sr);if(sr->trk.kalman.ntracks< 1||ibesttrk< 0||ibesttrk >=int(sr->trk.kalman.ntracks)) return-1000.f;const caf::SRKalmanTrackProxy &bestmuon=sr->trk.kalman.tracks[ibesttrk]; if(bestmuon.leninact > 0 && bestmuon.lenincat<=0) return float((bestmuon.leninact/100.) +(bestmuon.lenincat/100.)); if(bestmuon.leninact > 0 && bestmuon.lenincat > 0) return float(bestmuon.leninact/100.); if(bestmuon.leninact<=0 && bestmuon.lenincat > 0) return 0.f;return-1000.f;})
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Binning kTrkLenActBinning
const Var kTrueActiveE([](const caf::SRProxy *sr){return kTrueMuonE(sr)-kTrueCatcherE(sr);})
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const Cut kMuonIDCut([](const caf::SRProxy *sr){return kMuonID(sr) >=kMuonIDCutVal;})
enum BeamMode string