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"
36 
37 #include "CAFAna/Cuts/Cuts.h"
38 #include "CAFAna/Cuts/TruthCuts.h"
39 
44 
45 #include "TCanvas.h"
46 #include "TFile.h"
47 #include "TH1.h"
48 #include "TH2.h"
49 #include "TStyle.h"
50 
51 using namespace ana;
52 // using namespace ana::xsec::numubarcc;
53 
54 void estimate_energy(string outfileName = "numubar_energy_estimate.root")
55 {
56  // ND respin Monte Carlo files
58 
59  SpectrumLoader * loader = new SpectrumLoader(nd_mc);
61 
62  struct Plot2D
63  {
65  std::string labelx;
66  Binning binsx;
67  Var varx;
68  std::string labely;
69  Binning binsy;
70  Var vary;
71  Cut cut;
72  };
73 
75  // const Var kHadE = kTrueNuE - kIncXsecMuonE;
76 
77  const Cut kSelectCut = ana::xsec::numubarcc::kIsNueorbarCC &&
81 
82  const int kNumPlots2D = 3;
83 
84  Plot2D plots2D[kNumPlots2D] = {
85 
86  //---------------------------------------------------------------------------------
87  // Calculate for Muon Energy
88  // 1st : Reco Muon Track Length vs True Muon Energy
89  // 2nd : Reco Muon Track Length vs True Muon E - True Muon E in Catcher
90  // 3rd : Reco Muon Track Length in Catcher vs True Muon Energy in Catcher
91  //---------------------------------------------------------------------------------
92 
93  {"MuonE_hist_act","Reco Muon Track Length (m)",
96 
97  {"MuonE_hist_actAndCat","Reco Muon Track Length (m)",
98  ana::xsec::numubarcc::kTrkLenActBinning, ana::xsec::numubarcc::kTrkLenAct, "True Muon E - True Muon E in Catcher (GeV)" ,
100 
101  {"MuonE_hist_cat","Reco Muon Track Length in Catcher (m)",
102  ana::xsec::numubarcc::kTrkLenCatBinning, ana::xsec::numubarcc::kTrkLenCat, "True Muon Energy in Catcher (GeV)" ,
103  ana::xsec::numubarcc::kMuonECatBinning, ana::xsec::numubarcc::kTrueCatcherE, kSelectCut && ana::xsec::numubarcc::kActiveAndCatcher},
104 
105  // {"HadE","Visible Hadronic Energy (GeV)",
106  // kVisHadEBinning, kVisHadE, "True Nu E - Reco Muon Energy (GeV)" ,
107  // kMuonEActBinning, kHadE, kCut},
108  };
109 
110  Spectrum* sSpect2D[kNumPlots2D];
111 
112  for(int i = 0; i < kNumPlots2D; i++)
113  {
114  Plot2D p = plots2D[i];
115  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);
116  } // loop over 2D plots
117 
118  loader->Go();
119 
120  TFile * fout = TFile::Open(outfileName.c_str(),"RECREATE");
121 
122  for(int i = 0; i < kNumPlots2D; i++)
123  {
124  Plot2D p = plots2D[i];
125  sSpect2D[i]->SaveTo(fout, p.name);
126  } // loop over 2D plots
127 
128 } //end of void
129 
130 
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 char * p
Definition: xmltok.h:285
const Cut kMuonIDCut([](const caf::SRProxy *sr){return muonid_classifier::kMuonID(sr) >=kMuonIDCutVal;})
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 kActive([](const caf::SRProxy *sr){int ibesttrk=muonid_classifier::kBestMuonTrack(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:40
const Binning kTrkLenCatBinning
const Var kTrkLenCat([](const caf::SRProxy *sr){int ibesttrk=muonid_classifier::kBestMuonTrack(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;})
const Binning kMuonECatBinning
virtual void Go() override
Load all the registered spectra.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:506
loader
Definition: demo0.py:10
const Var kTrkLenAct([](const caf::SRProxy *sr){int ibesttrk=muonid_classifier::kBestMuonTrack(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;})
std::vector< float > Spectrum
Definition: Constants.h:759
void estimate_energy(string outfileName="numubar_energy_estimate.root")
const std::string nominal_dataset
Dataset definitions.
const SystShifts kNoShift
Definition: SystShifts.cxx:22
const Binning kMuonEActBinning
const ana::Var std_wgt
const Cut cut
Definition: exporter_fd.C:30
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Binning kTrkLenActBinning
const Cut kActiveAndCatcher([](const caf::SRProxy *sr){int ibesttrk=muonid_classifier::kBestMuonTrack(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 SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
enum BeamMode string