caf_numu_nuenergy_vs_xy.C
Go to the documentation of this file.
1 // Uses truth information. Only run on MC
2 
3 /*
4 */
5 #ifdef __CINT__
6 void caf_numu_nuenergy_vs_xy(std::string kInputFileName, std::string kOutputFileName, bool isFD)
7 {
8  std::cout << "Sorry, you must run in compiled mode" << std::endl;
9 }
10 #else
11 
12 #include "CAFAna/Core/Binning.h"
13 #include "CAFAna/Core/Spectrum.h"
16 
17 #include "CAFAna/Cuts/Cuts.h"
18 #include "CAFAna/Cuts/NumuCuts.h"
19 #include "CAFAna/Cuts/SpillCuts.h"
20 #include "CAFAna/Cuts/TruthCuts.h"
21 
22 #include "CAFAna/Vars/Vars.h"
23 #include "CAFAna/Vars/NumuVars.h"
24 #include "CAFAna/Vars/HistAxes.h"
25 
26 #include "CAFAna/Analysis/Plots.h"
27 #include "CAFAna/Analysis/Style.h"
28 
29 // #include "NumuEnergy/NumuEAlg.h"
30 
31 #include "TCanvas.h"
32 #include "TH1.h"
33 #include "TH2.h"
34 #include "TH3.h"
35 #include "TProfile2D.h"
36 #include "TLine.h"
37 #include "TFile.h"
38 #include "TStyle.h"
39 #include "TLegend.h"
40 #include "TLine.h"
41 #include "TLatex.h"
42 #include "THStack.h"
43 
44 using namespace ana;
45 
46 // Put a "NOvA Preliminary" tag in the corner
47 void Preliminary(){
48  TLatex* prelim = new TLatex(.9, .95, "NO#nuA Preliminary");
49  prelim->SetTextColor(kBlue);
50  prelim->SetNDC();
51  prelim->SetTextSize(2/30.);
52  prelim->SetTextAlign(32);
53  prelim->Draw();
54 }
55 
56 // const Cut kIsDytmanMEC({"mc.nnu", "mc.nu.mode"},
57 // [](const caf::StandardRecord* sr)
58 // {
59 // if(sr->mc.nnu == 0) return false;
60 // assert(sr->mc.nnu == 1);
61 // return (sr->mc.nu[0].mode == caf::kMEC);
62 // });
63 
64 template<class T> class Tangible{
65 
66 public:
67  Tangible(const T& obj, const std::string& shortName,
68  const std::string& blurb ):
69  fObj(obj),
70  fShortName(shortName),
71  fBlurb(blurb)
72  {};
73 
74  T fObj;
77 
78 };
79 
80 
83 
84 class UsefulHist{
85 
86 public:
87 
88  UsefulHist(Selection selA, TangibleAxis tanXAxis, TangibleAxis tanYAxis, TangibleAxis tanZAxis,
89  SpectrumLoader& loader, const Cut& cutMEC = !kIsDytmanMEC, const Cut& bkg=kNoCut):
90  fSelA(selA),
91  fXAxis(tanXAxis),
92  fYAxis(tanYAxis),
93  fZAxis(tanZAxis),
94  fHist(loader, fXAxis.fObj, fYAxis.fObj, fZAxis.fObj, fSelA.fObj && cutMEC),
95  fName("reco-" + fXAxis.fShortName + "_vs_" + fYAxis.fShortName + "_vs_" + fZAxis.fShortName + "-" + selA.fShortName ){};
96 
97  Selection fSelA;
101  Spectrum fHist;
102  std::string fName;
103 
104 };
105 
106 
107 void caf_numu_nuenergy_vs_xy(std::string kInputFileName, std::string kOutputFileName, bool isFD){
108 
109  std::cout << "\nrun : ---- Running CAF NuMU ND Validation.\n";
110  std::cout << "run : input file name: " << kInputFileName << std::endl;
111  std::cout << "run : output file name: " << kOutputFileName << std::endl;
112  std::cout << "Are we using FD binning/cuts? " << isFD << std::endl;
113 
114  // make loader with user input dataset
115  SpectrumLoader loader(kInputFileName);
116 
117  // make spill histogram tokeep record of spills. Going to use to count total number of events
118  TH1D* spillHist = new TH1D("reco-spills", ";Detector;Spills", 3, 0, 3);
119  SpillVar spillRun({"det"}, [](const caf::SRSpill* spill) {return spill->det;});
120 
121  // set spill cuts
123  loader.AddSpillHistogram(spillHist, spillRun, kStandardSpillCuts);
124 
125 
126  //////////////////// Cuts ////////////////////////
127  std::vector<Selection> selections_base;
128 
129  if(isFD) selections_base.emplace_back(kNumuFD, "numuFD",
130  "Selected events pass the numu FD selection");
131  else selections_base.emplace_back(kNumuND, "numuND",
132  "Selected events pass the numu ND selection");
133  ////////////////////////////////////////////////////////////////////////////////
134 
135  // Binning
136  const Binning kXYBins = Binning::Simple(55,-2.20,2.20); //ND
137  const Binning kZBins = Binning::Simple(32, 0,16.00); //ND
138  const Binning kFDXYBins = Binning::Simple(55,-8.0,8.0);
139  const Binning kFDZBins = Binning::Simple(55,0,60.00);
140  const Binning kEnergyBinning = Binning::Simple(50,0,5);
141 
142 
143  // variables
144  const Var kRecoMuE({"energy.mutrkE.E","sel.remid.bestidx","energy.mutrkE"},
145  [](const caf::StandardRecord* sr)
146  {
147  if(sr->sel.remid.bestidx < 0.f) return 0.f;
148  if(sr->sel.remid.bestidx == 999) return 0.f;
149  if(sr->energy.mutrkE.size() == 0) return 0.f;
150  return (sr -> energy.mutrkE[sr->sel.remid.bestidx].E);
151  });
152  const HistAxis axRecoMuE("Muon E. (GeV)", kEnergyBinning, kRecoMuE);
153 
154 
155  std::vector<TangibleAxis> variables;
156  variables.emplace_back(kNumuCCAxis, "numuE",
157  "CC energy estimator (GeV) ");
158 
159 
160  TangibleAxis tanTrkStartX (HistAxis("Track Start X Position (m)", kXYBins, kTrkStartX), "trkStartX", "Track start x position. ");
161  TangibleAxis tanTrkStartY (HistAxis("Track Start Y Position (m)", kXYBins, kTrkStartY), "trkStartY", "Track start y position. ");
162 
163 
164  std::vector<UsefulHist> hists;
165  hists.reserve(selections_base.size() * variables.size());
166 
167  for(const auto& sel_base:selections_base){
168  for(const auto& variable:variables){
169  hists.emplace_back(sel_base, tanTrkStartX, tanTrkStartY, variable, loader);
170  }
171  }
172 
173 
174  std::cout << "\nrun : --- run loader.\n";
175  loader.Go();
176  std::cout << "\nrun : --- done.\n";
177  double pot = hists[0].fHist.POT();
178 
179  std::cout << "\nrun : --- save output.\n";
180  TFile outputfile(kOutputFileName.c_str(), "RECREATE");
181 
182  for(const auto& hist:hists){
183  TH3 *temp = hist.fHist.ToTH3(pot);
184  std::string tempName = hist.fName;
185  TProfile2D* tempProj = (TProfile2D*)temp -> Project3DProfile("yx");
186  //temp->Write(tempName.c_str());
187  tempProj->Write(tempName.c_str());
188  }
189 
190  std::cout << "\nrun : --- save PoT.\n";
191 
192  TH1F *pothist = new TH1F("meta-TotalPOT", "TotalPOT", 1, 0, 1);
193  TH1F *evthist = new TH1F("meta-TotalEvents", "TotalEvents", 1, 0, 1);
194  pothist->SetBinContent(1, pot);
195  evthist->SetBinContent(1, spillHist->GetEntries());
196  pothist->Write("meta-TotalPOT");
197  evthist->Write("meta-TotalEvents");
198 
199  std::cout << "\nrun : --- close output files.\n";
200  outputfile.Close();
201  std::cout << "\nrun : --- done.\n";
202 
203 
204 
205  // dictionary to translate histogram names
206  /*
207 
208  */
209 
210  return;
211 }
212 
213 #endif
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
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 Var kTrkStartY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Y()/100;})
Definition: NumuVars.h:52
std::string fShortName
Definition: DataMCPair.h:47
Tangible< Cut > Selection
Definition: DataMCPair.h:52
void caf_numu_nuenergy_vs_xy(std::string kInputFileName, std::string kOutputFileName, bool isFD)
const Cut kNumuFD
Definition: NumuCuts.h:53
const Var kRecoMuE
void SetSpillCut(const SpillCut &cut)
TString hists[nhists]
Definition: bdt_com.C:3
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
void Preliminary()
Put NOvA Preliminary on plots.
const HistAxis kNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNumuEnergyBinning, kCCE)
Definition: HistAxes.h:8
virtual void AddSpillHistogram(TH1 *h, const SpillVar &var, const SpillCut &cut, const SpillVar &wei=kSpillUnweighted)
Uses include counting the total POT or spills in a run.
Tangible(const T &obj, const std::string &shortName, const std::string &blurb)
UsefulHist(Selection selA, TangibleAxis tanXAxis, TangibleAxis tanYAxis, TangibleAxis tanZAxis, SpectrumLoader &loader, const Cut &cutMEC=!kIsDytmanMEC, const Cut &bkg=kNoCut)
SRRemid remid
Output from RecoMuonID (ReMId) package.
Definition: SRIDBranch.h:39
Tangible< HistAxis > TangibleAxis
Definition: DataMCPair.h:53
const Cut kNumuND
Definition: NumuCuts.h:55
const Binning kEnergyBinning
double energy
Definition: plottest35.C:25
#define pot
const Binning kXYBins
Definition: VarsAndCuts.h:95
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
loader
Definition: demo0.py:10
TLatex * prelim
Definition: Xsec_final.C:133
const Var kTrkStartX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.X()/100;})
Definition: NumuVars.h:51
OStream cout
Definition: OStream.cxx:6
std::string fBlurb
Definition: DataMCPair.h:46
The StandardRecord is the primary top-level object in the Common Analysis File trees.
TangibleAxis fXAxis
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
SRIDBranch sel
Selector (PID) branch.
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
string tempName
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
double T
Definition: Xdiff_gwt.C:5
Template for Var and SpillVar.
SREnergyBranch energy
Energy estimator branch.
const Binning kZBins
Definition: VarsAndCuts.h:96
enum BeamMode kBlue
TangibleAxis fYAxis
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
enum BeamMode string