caf_numu_reco_minus_true.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_reco_minus_true(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 
25 #include "CAFAna/Analysis/Plots.h"
26 #include "CAFAna/Analysis/Style.h"
27 
28 // #include "NumuEnergy/NumuEAlg.h"
29 
30 #include "TCanvas.h"
31 #include "TH1.h"
32 #include "TH2.h"
33 #include "TLine.h"
34 #include "TFile.h"
35 #include "TStyle.h"
36 #include "TLegend.h"
37 #include "TLine.h"
38 #include "TLatex.h"
39 #include "THStack.h"
40 
41 using namespace ana;
42 
43 // Put a "NOvA Preliminary" tag in the corner
44 void Preliminary(){
45  TLatex* prelim = new TLatex(.9, .95, "NO#nuA Preliminary");
46  prelim->SetTextColor(kBlue);
47  prelim->SetNDC();
48  prelim->SetTextSize(2/30.);
49  prelim->SetTextAlign(32);
50  prelim->Draw();
51 }
52 
53 const Cut kIsDytmanMEC({"mc.nnu", "mc.nu.mode"},
54  [](const caf::StandardRecord* sr)
55  {
56  if(sr->mc.nnu == 0) return false;
57  assert(sr->mc.nnu == 1);
58  return (sr->mc.nu[0].mode == caf::kMEC);
59  });
60 
61 template<class T> class Tangible{
62 
63 public:
64  Tangible(const T& obj, const std::string& shortName,
65  const std::string& blurb ):
66  fObj(obj),
67  fShortName(shortName),
68  fBlurb(blurb)
69  {};
70 
71  T fObj;
74 
75 };
76 
77 
80 
81 
82 class UsefulHist{
83 
84 public:
85 
86  // what's wrong with doing the following?
87  // UsefulHist(Selection selA, TangibleAxis tanAxis, SpectrumLoader& loader, const Cut& bkg=kNoCut):
88  // fSelA(selA),
89  // fAxis(tanAxis),
90  // fHist(loader, fAxis.fObj, fSelA.fObj),
91  // fName(tanAxis.fShortName + "-" + selA.fShortName ){};
92 
94  fSelA(selA),
95  fSelB(selB),
96  fSelC(selC),
97  fAxis(tanAxis),
98  fHist(loader, fAxis.fObj, fSelA.fObj && fSelB.fObj && fSelC.fObj && !kIsDytmanMEC),
99  //fName(tanAxis.fShortName + "-" + selA.fShortName + "_" + selB.fShortName + "_" + selC.fShortName ){};
100  fName(tanAxis.fShortName + "-" + selA.fShortName + "_" + selB.fShortName + "_" + selC.fShortName ){};
101 
102  Selection fSelA;
105  TangibleAxis fAxis;
106  Spectrum fHist;
107  std::string fName;
108 
109 };
110 
111 
112 void caf_numu_reco_minus_true(std::string kInputFileName, std::string kOutputFileName, bool isFD){
113 
114  std::cout << "\nrun : ---- Running CAF NuMU ND Validation.\n";
115  std::cout << "run : input file name: " << kInputFileName << std::endl;
116  std::cout << "run : output file name: " << kOutputFileName << std::endl;
117  std::cout << "Are we using FD binning/cuts? " << isFD << std::endl;
118 
119  // make loader with user input dataset
120  SpectrumLoader loader(kInputFileName);
121 
122  // make spill histogram tokeep record of spills. Going to use to count total number of events
123  TH1D* spillHist = new TH1D("reco-spills", ";Detector;Spills", 3, 0, 3);
124  SpillVar spillRun({"det"}, [](const caf::SRSpill* spill) {return spill->det;});
125 
126  // set spill cuts
128  loader.AddSpillHistogram(spillHist, spillRun, kStandardSpillCuts);
129 
130 
131  //////////////////// Cuts ////////////////////////
132  // do cuts for contained and interaction type and true numu CC
133 
134  // kIsNumuCC
135  // kIsQE
136  // !kIsQE
137  // kNumuFD and kNumuND. Use ifFD boolean to switch
138  // kTruthContainedND. Only apply if ifFD==0
139  // !kTruthContainedND. Only apply if ifFD==0
140 
141  std::vector<Selection> selections_base;
142  std::vector<Selection> selections_int;
143  std::vector<Selection> selections_contain;
144 
145  //base
146  selections_base.emplace_back(kIsNumuCC, "true_numucc",
147  "Selected events are truly numu cc");
148  if(isFD) selections_base.emplace_back(kNumuFD, "numuFD",
149  "Selected events pass the numu FD selection");
150  else selections_base.emplace_back(kNumuND, "numuND",
151  "Selected events pass the numu ND selection");
152  //int.
153  selections_int.emplace_back(kNoCut, "no_cut",
154  "No selection applied");
155  selections_int.emplace_back(kIsQE, "true_qe",
156  "Selected events are truly quasi-elastic");
157  selections_int.emplace_back(!kIsQE, "true_non_qe",
158  "Selected events are truly non quasi-elastic");
159  //contain
160  selections_contain.emplace_back(kNoCut, "no_cut",
161  "No selection applied");
162 
163  selections_contain.emplace_back(kTruthContainedND, "true_nd_contained",
164  "Selected events are truly contained in the ND");
165  selections_contain.emplace_back(!kTruthContainedND, "true_nd_uncontained",
166  "Selected events are truly un-contained in the ND");
167 
168  ////////////////////////////////////////////////////////////////////////////////
169 
170 
171  //////////////////// Vars ////////////////////////
172  const Var kTrueNHitSlc({"mc.nu", "mc.nnu", "mc.nu.nhitslc"},
173  [](const caf::StandardRecord* sr)
174  {return (sr->mc.nnu == 0) ? 0 : sr->mc.nu[0].nhitslc;});
175  const Var kTrueNHitTot({"mc.nu", "mc.nnu", "mc.nu.nhittot"},
176  [](const caf::StandardRecord* sr)
177  {return (sr->mc.nnu == 0) ? 0 : sr->mc.nu[0].nhittot;});
178  const Var kRecoNHitTot({"slc.nhit"},
179  [](const caf::StandardRecord* sr)
180  {return sr->slc.nhit;});
181  const Var kTruehadE({"mc.nu", "mc.nnu", "mc.nu.*"},
182  [](const caf::StandardRecord* sr)
183  {return (sr->mc.nnu == 0 || sr->mc.nu[0].prim.size() == 0) ? 0 :
184  (sr->mc.nu[0].visE - sr->mc.nu[0].prim[0].visE);
185  });
186  const Var kTruehadESlc({"mc.nu", "mc.nnu", "mc.nu.*"},
187  [](const caf::StandardRecord* sr)
188  {return (sr->mc.nnu == 0 || sr->mc.nu[0].prim.size() == 0) ? 0 :
189  (sr->mc.nu[0].visEinslc - sr->mc.nu[0].prim[0].visEinslc);
190  });
191  const Var kRecohadE({"energy.numusimp.hadcalE"},
192  [](const caf::StandardRecord* sr)
193  {return (sr -> energy.numusimp.hadcalE);
194  });
195  const Var kTrueMuE({"mc.nu", "mc.nnu", "mc.nu.*"},
196  [](const caf::StandardRecord* sr)
197  {return (sr->mc.nu[0].prim.size() == 0) ? 0 :
198  (sr->mc.nu[0].prim[0].p.E());
199  });
200  const Var kTrueMuESlc({"mc.nu", "mc.nnu", "mc.nu.*"},
201  [](const caf::StandardRecord* sr)
202  {return (sr->mc.nu[0].prim.size() == 0) ? 0 :
203  sr->mc.nu[0].prim[0].p.E();
204  });
205  const Var kRecoMuE({"energy.mutrkE.E","sel.remid.bestidx"},
206  [](const caf::StandardRecord* sr)
207  {
208  if(sr->sel.remid.bestidx < 0.f) return 0.f;
209  if(sr->sel.remid.bestidx == 999) return 0.f;
210  return (sr -> energy.mutrkE[sr->sel.remid.bestidx].E);
211  });
212  const Var kTrueNuEnergy({"mc.nu", "mc.nnu", "mc.nu.E"},
213  [](const caf::StandardRecord* sr)
214  {return (sr->mc.nnu == 0) ? 0 :
215  sr->mc.nu[0].E;
216  });
217  const Var kRecoNuEnergy({"energy.numusimp.trkccE"},
218  [](const caf::StandardRecord* sr)
219  {return (sr->mc.nnu == 0) ? 0 :
220  sr->energy.numusimp.trkccE;
221  });
222 
223 
224  // reco - true / true
225  const Var kRMTOT_NHitSlc({"mc.nu", "mc.nnu", "mc.nu.nhitslc", "slc.nhit"},
226  [](const caf::StandardRecord* sr)
227  {
228  if(sr->mc.nnu == 0) return -10.;
229 
230  double nhitReco = (double)sr->slc.nhit;
231  double nhitTrue = (double)sr->mc.nu[0].nhitslc;
232  return ((nhitReco - nhitTrue) / nhitTrue);
233 
234  });
235 
236 
237  const Var kRMTOT_Energy({"mc.nu", "mc.nnu", "mc.nu.E", "energy.numusimp.trkccE"},
238  [](const caf::StandardRecord* sr)
239  {return (sr->mc.nnu == 0) ? -10 :
240  ( (sr->energy.numusimp.trkccE - sr->mc.nu[0].E) / (sr->mc.nu[0].E) );
241  });
242  //ehad
243  const Var kRMTOT_hadE({"mc.nu", "mc.nnu", "mc.nu.*", "energy.numusimp.hadcalE"},
244  [](const caf::StandardRecord* sr)
245  {
246  if(sr->mc.nnu == 0) return -10.f;
247  if(sr->mc.nu[0].prim.size() == 0) return -10.f;
248  float trueHadE = (sr -> mc.nu[0].E) - (sr -> mc.nu[0].prim[0].p.E());
249  float muE=sr->energy.mutrkE[sr->sel.remid.bestidx].E;
250  float recoHadE=sr->energy.numusimp.trkccE-muE;
251 
252  if(trueHadE == 0) return -10.f;
253  return ( (recoHadE - trueHadE) / trueHadE );
254  });
255  // E_mu
256  const Var kRMTOT_Emu({"mc.nu", "mc.nnu", "mc.nu.*", "energy.mutrkE.E", "sel.remid.bestidx"},
257  [](const caf::StandardRecord* sr)
258  {
259  if(sr->mc.nnu == 0) return -10.f;
260  if(sr->mc.nu[0].prim.size() == 0) return -10.f;
261  if(sr->sel.remid.bestidx < 0.f) return -10.f;
262  if(sr->sel.remid.bestidx == 999) return -10.f;
263  float trueMuE = (sr -> mc.nu[0].prim[0].p.E());
264  float recoMuE = (sr -> energy.mutrkE[sr->sel.remid.bestidx].E);
265  if(trueMuE == 0) return -10.f;
266  return ( (recoMuE - trueMuE) / trueMuE );
267  });
268 
269  ////////////////////////////////////////////////////////////////////////////////
270 
271  // Binning
272  const Binning kXYBins = Binning::Simple(55,-2.20,2.20);
273  const Binning kZBins = Binning::Simple(32, 0,16.00);
274  const Binning kFDXYBins = Binning::Simple(55,-8.0,8.0);
275  const Binning kFDZBins = Binning::Simple(55,0,60.00);
276  const Binning kEnergyBinning = Binning::Simple(50,0,5);
278 
279  //HistAxis
280  const HistAxis axRMTOT_NHitSlc("nhit in slc (Reco. - true)/true", Binning::Simple(100,-2,2), kRMTOT_NHitSlc);
281  const HistAxis axRMTOT_Energy("Neutrino energy (Reco. - true)/true", Binning::Simple(100,-2,2), kRMTOT_Energy);
282  const HistAxis axRMTOT_HadE("Hadronic energy (Reco. - true)/true", Binning::Simple(100,-2,2), kRMTOT_hadE);
283  const HistAxis axRMTOT_Emu("Muon energy (Reco. - true)/true", Binning::Simple(100,-2,2), kRMTOT_Emu);
284 
285  const HistAxis axTrueNHit("NHits", Binning::Simple(100,0,500), kTrueNHitTot);
286  const HistAxis axTrueNHitSlc("NHits", Binning::Simple(100,0,500), kTrueNHitSlc);
287  const HistAxis axRecoNHit("NHits", Binning::Simple(100,0,500), kRecoNHitTot);
288  const HistAxis axTrueHadE("Had E. (GeV)", kHadronicEnergyBinning, kTruehadE);
289  const HistAxis axTrueHadESlc("Had E. (GeV)", kHadronicEnergyBinning, kTruehadESlc);
290  const HistAxis axRecoHadE("Had E. (GeV)", kHadronicEnergyBinning, kRecohadE);
291  const HistAxis axTrueMuE("Muon E. (GeV)", kEnergyBinning, kTrueMuE);
292  const HistAxis axTrueMuESlc("Muon E. (GeV)", kEnergyBinning, kTrueMuESlc);
293  const HistAxis axRecoMuE("Muon E. (GeV)", kEnergyBinning, kRecoMuE);
294  const HistAxis axTrueNuEnergy("Neutrino E. (GeV)", kEnergyBinning, kTrueNuEnergy);
295  const HistAxis axRecoNuEnergy("Neutrino E. (GeV)", kEnergyBinning, kRecoNuEnergy);
296 
297  ////////////////////////////////////////////////////////////////////////////////
298 
299  std::vector<TangibleAxis> variables;
300  variables.emplace_back(axRMTOT_NHitSlc, "res-nhit_slc_reco_minus_true_over_true",
301  "Number of hits in slc. (reco-true)/true");
302  variables.emplace_back(axRMTOT_Energy, "res-neutrino_energy_reco_minus_true_over_true",
303  "Neutrino energy (reco-true)/true");
304  variables.emplace_back(axRMTOT_HadE, "res-hadronic_energy_reco_minus_true_over_true",
305  "Hadronic energy (reco-true)/true");
306  variables.emplace_back(axRMTOT_Emu, "res-muon_energy_reco_minus_true_over_true",
307  "Muon energy (reco-true)/true");
308 
309  std::vector<TangibleAxis> dist_variables;
310  // include these when the variables are corrected
311  // dist_variables.emplace_back(axTrueNHit, "true-nhit_truth", "NHits");
312  // dist_variables.emplace_back(axTrueNHitSlc, "true-nhit_slc_truth", "NHits");
313  // dist_variables.emplace_back(axRecoNHit, "true-nhit_reco", "NHits");
314  // dist_variables.emplace_back(axTrueHadE, "true-hade_truth", "HadE (GeV)");
315  // dist_variables.emplace_back(axTrueHadESlc, "true-hade_slc_truth", "HadE (GeV)");
316  // dist_variables.emplace_back(axRecoHadE, "true-hade_reco", "HadE (GeV)");
317  // dist_variables.emplace_back(axTrueMuE, "true-mue_truth", "Mu E (GeV)");
318  // dist_variables.emplace_back(axTrueMuESlc, "true-mue_slc_truth", "Mu E (GeV)");
319  // dist_variables.emplace_back(axRecoMuE, "true-mue_reco", "Mu E (GeV)");
320  // dist_variables.emplace_back(axTrueNuEnergy, "true-nu_energy_truth", "Nu E (GeV)");
321  // dist_variables.emplace_back(axRecoNuEnergy, "true-nu_energy_reco", "Nu E (GeV)");
322 
323  ////////////////////////////////////////////////////////////////////////////////
324 
325  std::vector<UsefulHist> hists;
326  hists.reserve(selections_base.size() * selections_int.size() * selections_contain.size() * variables.size() + selections_base.size() * selections_int.size() * selections_contain.size() * dist_variables.size());
327 
328 
329  if(!isFD){
330  for(const auto& sel_base:selections_base){
331  for(const auto& sel_int:selections_int){
332  for(const auto& sel_contain:selections_contain){
333  for(const auto& variable:variables){
334  hists.emplace_back(sel_base, sel_int, sel_contain, variable, loader);
335  }
336  // include these when the variables are corrected
337  // for(const auto& dist:dist_variables){
338  // hists.emplace_back(sel_base, sel_int, sel_contain, dist, loader);
339  // }
340  }
341  }
342  }
343  }
344 
345  if(isFD){
346  for(const auto& sel_base:selections_base){
347  for(const auto& sel_int:selections_int){
348  for(const auto& variable:variables){
349  hists.emplace_back(sel_base, sel_int, selections_contain[0], variable, loader);
350  }
351  // for(const auto& dist:dist_variables){
352  // hists.emplace_back(sel_base, sel_int, selections_contain[0], dist, loader);
353  // }
354  }
355  }
356  }
357 
358  std::cout << "\nrun : --- run loader.\n";
359  loader.Go();
360  std::cout << "\nrun : --- done.\n";
361  double pot = hists[0].fHist.POT();
362 
363  std::cout << "\nrun : --- save output.\n";
364  // TFile inputfile(kInputFileName.c_str(), "READ");
365  TFile outputfile(kOutputFileName.c_str(), "RECREATE");
366 
367  for(const auto& hist:hists){
368  TH1 *temp = hist.fHist.ToTH1(pot);
369  std::string tempName = hist.fName;
370  temp->Write(tempName.c_str());
371  }
372 
373  std::cout << "\nrun : --- save PoT.\n";
374 
375  TH1F *pothist = new TH1F("meta-TotalPOT", "TotalPOT", 1, 0, 1);
376  TH1F *evthist = new TH1F("meta-TotalEvents", "TotalEvents", 1, 0, 1);
377  pothist->SetBinContent(1, pot);
378  evthist->SetBinContent(1, spillHist->GetEntries());
379  pothist->Write("meta-TotalPOT");
380  evthist->Write("meta-TotalEvents");
381 
382  std::cout << "\nrun : --- close output files.\n";
383  outputfile.Close();
384  // inputfile.Close(); // will produce error if running over SAM dataset
385  std::cout << "\nrun : --- done.\n";
386 
387  // dictionary to translate histogram names
388  /*
389 
390  */
391 
392  return;
393 }
394 
395 #endif
const Cut kIsQE
Definition: TruthCuts.cxx:104
const Cut kTruthContainedND([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;if(fabs(sr->mc.nu[0].vtx.X()) > 200|| fabs(sr->mc.nu[0].vtx.Y()) > 200|| sr->mc.nu[0].vtx.Z()< 0|| sr->mc.nu[0].vtx.Z() > 1650) return false;return true;})
Definition: Cuts.h:88
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
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
std::string fShortName
Definition: DataMCPair.h:47
Tangible< Cut > Selection
Definition: DataMCPair.h:52
void caf_numu_reco_minus_true(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 Var kTrueMuE
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
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)
const Binning kHadronicEnergyBinning
SRRemid remid
Output from RecoMuonID (ReMId) package.
Definition: SRIDBranch.h:39
const Var kTrueNuEnergy
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
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
unsigned int nhit
number of hits
Definition: SRSlice.h:22
TLatex * prelim
Definition: Xsec_final.C:133
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.
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
UsefulHist(Selection selA, Selection selB, Selection selC, TangibleAxis tanAxis, SpectrumLoader &loader, const Cut &bkg=kNoCut)
SRIDBranch sel
Selector (PID) branch.
assert(nhit_max >=nhit_nbins)
void Preliminary()
Put NOvA Preliminary on plots.
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
SRSlice slc
Slice branch: nhit, extents, time, etc.
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
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
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
std::vector< SRNeutrino > nu
implemented as a vector to maintain mc.nu structure, i.e. not a pointer, but with 0 or 1 entries...
Definition: SRTruthBranch.h:25
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
enum BeamMode string