caf_numu_validation.C
Go to the documentation of this file.
1 // Fills spectra for systematic error band
2 #ifdef __CINT__
3 void caf_numu_validation(std::string kInputFileName, std::string kOutputFileName, bool isFD)
4 {
5  std::cout << "Sorry, you must run in compiled mode" << std::endl;
6 }
7 #else
8 
9 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Cuts/NumuCuts.h"
11 #include "CAFAna/Core/Binning.h"
12 #include "CAFAna/Core/Spectrum.h"
13 #include "CAFAna/Core/Var.h"
15 #include "CAFAna/Systs/Systs.h"
16 #include "CAFAna/Systs/NumuSysts.h"
17 #include "CAFAna/Vars/NumuVars.h"
18 #include "CAFAna/Cuts/SpillCuts.h"
19 #include "CAFAna/Analysis/Plots.h"
20 #include "CAFAna/Analysis/Style.h"
21 #include "CAFAna/Vars/HistAxes.h"
22 //#include "Utilities/rootlogon.C"
23 
24 #include "TCanvas.h"
25 #include "TFile.h"
26 #include "TGraph.h"
27 #include "TH1.h"
28 #include "TH2.h"
29 #include "TMath.h"
30 #include "TGaxis.h"
31 #include "TMultiGraph.h"
32 #include "TLegend.h"
33 #include "TLatex.h"
34 #include "TStyle.h"
35 #include "THStack.h"
36 #include "TPaveText.h"
37 #include "TList.h"
38 #include "TGaxis.h"
39 #include "TAttLine.h"
40 #include "TAttMarker.h"
41 
42 #include <cmath>
43 #include <iostream>
44 #include <vector>
45 #include <list>
46 #include <sstream>
47 #include <string>
48 #include <sstream>
49 #include <fstream>
50 #include <iomanip>
51 
52 
53 using namespace ana;
54 
55 // Put a "NOvA Preliminary" tag in the corner
56 void Preliminary(){
57  TLatex* prelim = new TLatex(.9, .95, "NO#nuA Preliminary");
58  prelim->SetTextColor(kBlue);
59  prelim->SetNDC();
60  prelim->SetTextSize(2/30.);
61  prelim->SetTextAlign(32);
62  prelim->Draw();
63 }
64 
65 template<class T> class Tangible{
66 
67 public:
68  Tangible(const T& obj, const std::string& shortName,
69  const std::string& blurb ):
70  fObj(obj),
71  fShortName(shortName),
72  fBlurb(blurb)
73  {};
74 
75  T fObj;
78 
79 };
80 
81 const Cut kIsDytmanMEC({"mc.nnu", "mc.nu.mode"},
82  [](const caf::StandardRecord* sr)
83  {
84  if(sr->mc.nnu == 0) return false;
85  assert(sr->mc.nnu == 1);
86  return (sr->mc.nu[0].mode == caf::kMEC);
87  });
88 
91 
94 
95 const Var kNumuHadTrkE = SIMPLEVAR(energy.numusimp.hadtrkE);
96 const Var kQePId = SIMPLEVAR(sel.qepid.pid);
97 const Var kNonQeE = SIMPLEVAR(energy.numusimp.trknonqeE);
98 const Var kQepidOffE = SIMPLEVAR(sel.qepid.offE);
99 const Var kQepidEDiff = SIMPLEVAR(sel.qepid.ediff);
100 const Var kQepidEDiffZ = SIMPLEVAR(sel.qepid.ediffz);
101 const Var kQepidDedx = SIMPLEVAR(sel.qepid.dedx);
102 const Var kQeE = SIMPLEVAR(energy.numusimp.trkqeE);
103 const Var kQeAngleE = SIMPLEVAR(energy.numusimp.angleE);
104 
105 const Var kQeHadE({"energy.numusimp.trkqeE", "energy.mutrkE", "trk.nkalman",
106  "sel.remid.bestidx"},
107  [](const caf::StandardRecord *sr)
108  {
109  if (sr->trk.nkalman > 0 && sr->sel.remid.bestidx != 999)
110  return sr->energy.numusimp.trkqeE - sr->energy.mutrkE[sr->sel.remid.bestidx].E;
111  return -5.f;
112  });
113 
114 const Var kNumuTrackE({"energy.mutrkE", "trk.nkalman", "sel.remid.bestidx"},
115  [](const caf::StandardRecord *sr)
116  {
117  if (sr->trk.nkalman > 0 && sr->sel.remid.bestidx != 999)
118  return sr->energy.mutrkE[sr->sel.remid.bestidx].E;
119  return -5.f;
120  });
121 
122 const Var kNumuHadVisE({"energy.numusimp.hadcalE", "energy.numusimp.hadtrkE"},
123  [](const caf::StandardRecord *sr)
124  {
125  return sr->energy.numusimp.hadcalE + sr->energy.numusimp.hadtrkE;
126  });
127 const Var kSlcUnCalibNHit({"slc.nhit", "slc.ncalhit"},
128  [](const caf::StandardRecord *sr)
129  {
130  return sr->slc.nhit - sr->slc.ncalhit;
131  });
132 
133  const Var kSlcMeanTime({"slc.meantime"},
134  [](const caf::StandardRecord *sr)
135  {
136  return sr->slc.meantime / 1000.; // Nicer in \mu s
137  });
138 
139  const Var kSlcStartTime({"slc.starttime"},
140  [](const caf::StandardRecord *sr)
141  {
142  return sr->slc.starttime / 1000.; // Nicer in \mu s
143  });
144 
145  const Var kSlcEndTime({"slc.endtime"},
146  [](const caf::StandardRecord *sr)
147  {
148  return sr->slc.endtime / 1000.; // Nicer in \mu s
149  });
150 
151  const Var kSlcMinX({"slc.boxmin.fX"},
152  [](const caf::StandardRecord *sr)
153  {
154  return sr->slc.boxmin.X() / 100.; // Nicer in m
155  });
156 
157  const Var kSlcMinY({"slc.boxmin.fY"},
158  [](const caf::StandardRecord *sr)
159  {
160  return sr->slc.boxmin.Y() / 100.; // Nicer in m
161  });
162 
163  const Var kSlcMinZ({"slc.boxmin.fZ"},
164  [](const caf::StandardRecord *sr)
165  {
166  return sr->slc.boxmin.Z() / 100.; // Nicer in m
167  });
168 
169  const Var kSlcMaxX({"slc.boxmax.fX"},
170  [](const caf::StandardRecord *sr)
171  {
172  return sr->slc.boxmax.X() / 100.; // Nicer in m
173  });
174 
175  const Var kSlcMaxY({"slc.boxmax.fY"},
176  [](const caf::StandardRecord *sr)
177  {
178  return sr->slc.boxmax.Y() / 100.; // Nicer in m
179  });
180 
181  const Var kSlcMaxZ({"slc.boxmax.fZ"},
182  [](const caf::StandardRecord *sr)
183  {
184  return sr->slc.boxmax.Z() / 100.; // Nicer in m
185  });
186 
187  const Var kSlcExtentX({"slc.boxmax.fX", "slc.boxmin.fX"},
188  [](const caf::StandardRecord *sr)
189  {
190  return (sr->slc.boxmax.X() - sr->slc.boxmin.X())/ 100.; // Nicer in m
191  });
192 
193  const Var kSlcExtentY({"slc.boxmax.fY", "slc.boxmin.fY"},
194  [](const caf::StandardRecord *sr)
195  {
196  return (sr->slc.boxmax.Y() - sr->slc.boxmin.Y())/ 100.; // Nicer in m
197  });
198 
199  const Var kSlcExtentZ({"slc.boxmax.fZ", "slc.boxmin.fZ"},
200  [](const caf::StandardRecord *sr)
201  {
202  return (sr->slc.boxmax.Z() - sr->slc.boxmin.Z())/ 100.; // Nicer in m
203  });
204 
205  const Var kTrkNPlaneGap({"trk.kalman.nplanegap", "trk.nkalman", "sel.remid.bestidx"},
206  [](const caf::StandardRecord *sr)
207  {
208  if (sr->trk.nkalman > 0 && sr->sel.remid.bestidx != 999)
209  return sr->trk.kalman[sr->sel.remid.bestidx].nplanegap;
210  return (unsigned short) 500;
211  });
212  const Var kCosNumi({"hdr.det", "sel.remid.bestidx", "trk.kalman.dir.fX", "trk.kalman.dir.fY",
213  "trk.kalman.dir.fZ"},
214  [](const caf::StandardRecord *sr)
215  {
216  if (sr->trk.nkalman > 0 && sr->sel.remid.bestidx != 999){
217  if (sr->hdr.det == 1){
218  return sr->trk.kalman[sr->sel.remid.bestidx].dir.Dot(beamDirND);
219  }
220  if (sr->hdr.det == 2){
221  return sr->trk.kalman[sr->sel.remid.bestidx].dir.Dot(beamDirFD);
222  }
223  }
224  return -5.;
225  });
226  const Var kDirX({"trk.nkalman", "trk.kalman.dir.fX", "sel.remid.bestidx"},
227  [](const caf::StandardRecord *sr)
228  {
229  if(sr->trk.nkalman < 1) return -5.;
230  return sr->trk.kalman[sr->sel.remid.bestidx].dir.X();
231  });
232 
233  const Var kDirY({"trk.nkalman", "trk.kalman.dir.fY", "sel.remid.bestidx"},
234  [](const caf::StandardRecord *sr)
235  {
236  if(sr->trk.nkalman < 1) return -5.;
237  return sr->trk.kalman[sr->sel.remid.bestidx].dir.Y();
238  });
239  const Var kSlcCalEPerNHit({"slc.nhit", "slc.calE"},
240  [](const caf::StandardRecord *sr)
241  {
242  if (sr->slc.nhit > 0) return sr->slc.calE / (1.78 * sr->slc.nhit);
243  return -5.;
244  });
245  const Var kNonQeHadE({"energy.numusimp.trknonqeE", "energy.mutrkE", "trk.nkalman",
246  "sel.remid.bestidx"},
247  [](const caf::StandardRecord *sr)
248  {
249  if (sr->trk.nkalman > 0 && sr->sel.remid.bestidx != 999)
250  return sr->energy.numusimp.trknonqeE - sr->energy.mutrkE[sr->sel.remid.bestidx].E;
251  return -5.f;
252  });
253 
254 
255 class UsefulHist{
256 
257 public:
258 
260  fSel(sel),
261  fAxis(tanAxis),
262  fHist(loader, fAxis.fObj, fSel.fObj),
263  fName(tanAxis.fShortName + "-" + sel.fShortName){};
264 
265  Selection fSel;
266  TangibleAxis fAxis;
267  Spectrum fHist;
268  std::string fName;
269 
270 };
271 
272 void caf_numu_validation(std::string kInputFileName, std::string kOutputFileName, bool isFD)
273 {
274  std::cout << "\nrun : ---- Running CAF NuMU Validation.\n";
275  std::cout << "run : input file name: " << kInputFileName << std::endl;
276  std::cout << "run : output file name: " << kOutputFileName << std::endl;
277  std::cout << "run : FD binning/cuts: " << isFD << std::endl;
278 
279  gStyle->SetMarkerStyle(kFullCircle);
280  TGaxis::SetMaxDigits(3);
281 
282  SpectrumLoader loader(kInputFileName);
283 
284  TH1D* spillHist = new TH1D("reco-spills", ";Detector;Spills", 3, 0, 3);
285  SpillVar spillRun({"det"}, [](const caf::SRSpill* spill) {return spill->det;});
286 
288  loader.AddSpillHistogram(spillHist, spillRun, kStandardSpillCuts);
289 
290  const Binning kXYBins = Binning::Simple(55,-2.20,2.20);
291  const Binning kZBins = Binning::Simple(32, 0,16.00);
292  const Binning kFDXYBins = Binning::Simple(55,-8.0,8.0);
293  const Binning kFDZBins = Binning::Simple(55,0,60.00);
294  const Binning kEnergyBinning = Binning::Simple(50,0,5);
296 
297  std::vector<Selection> selections;
298 
299  if (isFD){
300  selections.emplace_back((!kIsDytmanMEC) && kNumuContainFD && kNumuQuality && kNumuCosmicRej, "NumuContainFD",
301  "Selected events pass numu containment, cosmic rejection and slice quality cuts. ");
302  selections.emplace_back((!kIsDytmanMEC) && kNumuContainFD && kNumuQuality && kNumuCosmicRej && kNumuNCRej, "NumuFD",
303  "Selected events pass numu containment, cosmic rejection, slice quality, and ReMId cuts. ");
304  }
305  else{
306  selections.emplace_back((!kIsDytmanMEC) && kNumuContainND && kNumuQuality, "NumuContainND",
307  "Selected events pass numu containment and slice quality cuts. ");
308  selections.emplace_back((!kIsDytmanMEC) && kNumuContainND && kNumuQuality && kNumuNCRej, "NumuND",
309  "Selected events pass numu containment, slice quality, and ReMId cuts. ");
310  }
311 
312  std::vector<Selection> selectionsForQePId;
313 
314  if (isFD){
315  selectionsForQePId.emplace_back((!kIsDytmanMEC) && kNumuContainFD && kNumuQuality && kNumuCosmicRej && kNumuNCRej
316  && SIMPLEVAR(trk.nkalman)==1, "NumuFD1trk",
317  "Selected events pass numu containment, cosmic rejection, slice quality, ReMId cuts and have 1 track. ");
318 
319  selectionsForQePId.emplace_back((!kIsDytmanMEC) && kNumuContainFD && kNumuQuality && kNumuCosmicRej && kNumuNCRej
320  && SIMPLEVAR(trk.nkalman)==2, "NumuFD2trk",
321  "Selected events pass numu containment, cosmic rejection, slice quality, ReMId cuts and have 2 tracks. ");
322  }
323  else{
324  selectionsForQePId.emplace_back((!kIsDytmanMEC) && kNumuContainND && kNumuQuality && kNumuNCRej
325  && SIMPLEVAR(trk.nkalman)==1, "NumuND1trk",
326  "Selected events pass numu containment, slice quality, ReMId cuts and have 1 track. ");
327 
328  selectionsForQePId.emplace_back((!kIsDytmanMEC) && kNumuContainND && kNumuQuality && kNumuNCRej
329  && SIMPLEVAR(trk.nkalman)==2, "NumuND2trk",
330  "Selected events pass numu containment, slice quality, ReMId cuts and have 2 tracks. ");
331  }
332 
333  std::vector<TangibleAxis> variables;
334 
335  variables.emplace_back(
336  HistAxis("Reconstructed Neutrino Energy [GeV]", kEnergyBinning, kCCE),
337  "reco-numuE", "CC energy estimator. ");
338 
339  variables.emplace_back(
340  HistAxis("Slice N_{Hit}", Binning::Simple(50, 0, 500), kNHit),
341  "reco-slcNHit", "Number of hits in slice. " );
342 
343  variables.emplace_back(
344  HistAxis("Slice Calorimetric Energy [GeV]", kEnergyBinning, kCaloE),
345  "reco-calE", "Calorimetric energy of slice. " );
346 
347  variables.emplace_back(
348  HistAxis("Reconstructed Hadronic Energy [GeV]", kHadronicEnergyBinning, kHadE),
349  "reco-hadE", "Hadronic enery, i.e. numu energy estimate minus muon track energy. " );
350 
351  variables.emplace_back(
352  HistAxis("Muon Track Energy [GeV]", Binning::Simple(50, 0, 5), kNumuTrackE),
353  "reco-numuTrackE", "Muon track Energy. " );
354 
355  variables.emplace_back(
356  HistAxis("Average Hadronic Energy Per Hit [GeV]", Binning::Simple(40, 0, 0.04), kHadEPerNHit),
357  "reco-hadEPerNHit", "Average energy per hit in hadronic cluster, i.e. had_E/had_n_hit. " );
358 
359  variables.emplace_back(
360  HistAxis("Average Track Energy Per Hit [GeV]", Binning::Simple(40, 0, 0.04), kTrkEPerNHit),
361  "reco-trkEPerNHit", "Average energy per hit on primary track, i.e. trk_E/trk_n_hit. " );
362 
363  variables.emplace_back(
364  HistAxis("Hadronic N_{Hit}", Binning::Simple(50, 0, 100), kHadNHit),
365  "reco-hadNHit", "Number of hits in hadronic cluster. ");
366 
367  variables.emplace_back(
368  HistAxis("Off-track Calorimetric Hadronic Energy [GeV]", kHadronicEnergyBinning, kNumuHadCalE),
369  "reco-hadCalE", "Sum of calibrated energy deposits in hadronic cluster. " );
370 
371  variables.emplace_back(
372  HistAxis("On-track Calorimetric Hadronic Energy [GeV]", Binning::Simple(50, 0, 0.5), kNumuHadTrkE),
373  "reco-hadTrkE", "Sum of calibrated hadronic energy deposits on the muon track. " );
374 
375  variables.emplace_back(
376  HistAxis("Calorimetric Hadronic Energy [GeV]", kHadronicEnergyBinning, kNumuHadVisE),
377  "reco-hadVisE", "Sum of calibrated hadronic energy, both on and off the muon track. " );
378 
379  if(isFD){
380  variables.emplace_back(
381  HistAxis("Slice Maximum Y [m]", kFDXYBins, kSlcMaxY),
382  "reco-maxy", "Maximum Y position of slice. " );
383  }
384 
385  else{
386  variables.emplace_back(
387  HistAxis("Slice Maximum Y [m]", kXYBins, kSlcMaxY),
388  "reco-maxy", "Maximum Y position of slice. " );
389  }
390 
391  variables.emplace_back(
392  HistAxis("Number of Tracks in Slice", Binning::Simple(14, 1, 15), kNKalman),
393  "reco-nkal", "Number of tracks in slice. " );
394 
395  variables.emplace_back(
396  HistAxis("Number of Hits in Slice", Binning::Simple(50, 0, 500), kNHit),
397  "reco-nhit", "Number of hits in slice. " );
398 
399  if(isFD){
400  variables.emplace_back(
401  HistAxis("Track Start X Position [m]", kFDXYBins, kTrkStartX),
402  "reco-trkStartX", "Track start x position. " );
403 
404  variables.emplace_back(
405  HistAxis("Track Start Y Position [m]", kFDXYBins, kTrkStartY),
406  "reco-trkStartY", "Track start y position. " );
407 
408  variables.emplace_back(
409  HistAxis("Track Start Z Position [m]", kFDZBins, kTrkStartZ),
410  "reco-trkStartZ", "Track start z position. " );
411 
412  variables.emplace_back(
413  HistAxis("Track End X Position [m]", kFDXYBins, kTrkEndX),
414  "reco-trkEndX", "Track stop x position. " );
415 
416  variables.emplace_back(
417  HistAxis("Track End Y Position [m]", kFDXYBins, kTrkEndY),
418  "reco-trkEndY", "Track stop y position. " );
419 
420  variables.emplace_back(
421  HistAxis("Track End Z Position [m]", kFDZBins, kTrkEndZ),
422  "reco-trkEndZ", "Track stop z position. " );
423 
424  variables.emplace_back(
425  HistAxis("Slice Duration [ns]", Binning::Simple(50,0,3000), kSliceDuration),
426  "reco-sliceDuration", "Slice duration. " );
427  }
428 
429  else{
430  variables.emplace_back(
431  HistAxis("Track Start X Position [m]", kXYBins, kTrkStartX),
432  "reco-trkStartX", "Track start x position. " );
433 
434  variables.emplace_back(
435  HistAxis("Track Start Y Position [m]", kXYBins, kTrkStartY),
436  "reco-trkStartY", "Track start y position. " );
437 
438  variables.emplace_back(
439  HistAxis("Track Start Z Position [m]", kZBins, kTrkStartZ),
440  "reco-trkStartZ", "Track start z position. " );
441 
442  variables.emplace_back(
443  HistAxis("Track End X Position [m]", kXYBins, kTrkEndX),
444  "reco-trkEndX", "Track stop x position. " );
445 
446  variables.emplace_back(
447  HistAxis("Track End Y Position [m]", kXYBins, kTrkEndY),
448  "reco-trkEndY", "Track stop y position. " );
449 
450  variables.emplace_back(
451  HistAxis("Track End Z Position [m]", kZBins, kTrkEndZ),
452  "reco-trkEndZ", "Track stop z position. " );
453 
454  variables.emplace_back(
455  HistAxis("Slice Duration [ns]", Binning::Simple(50,0,500), kSliceDuration),
456  "reco-sliceDuration", "Slice duration. " );
457  }
458 
459  variables.emplace_back(
460  HistAxis("Number of Hits on Primary Track", Binning::Simple(50,0,500), kTrkNhits),
461  "reco-trkNhits", "Number of hits on primary track. ");
462 
463  variables.emplace_back(
464  HistAxis("Primary Track Length [m]", Binning::Simple(50,0,16), kTrkLength),
465  "reco-trkLength", "Primary track length. " );
466 
467  variables.emplace_back(
468  HistAxis("ReMId Input: Scattering Log-likelihood", Binning::Simple(50,-0.5,0.5), kReMIdScatLLH),
469  "reco-scatLL", "ReMId scattering log log-likelihood for primary track. " );
470 
471  variables.emplace_back(
472  HistAxis("ReMId Input: dE/dx Log-likelihood", Binning::Simple(50,-3,1), kReMIdDEDxLLH),
473  "reco-dedxLL", "ReMId dE/dx log log-likelihood for primary track. " );
474 
475  variables.emplace_back(
476  HistAxis("ReMId Input: Non-hadronic Plane Fraction", Binning::Simple(50,0,1), kReMIdMeasFrac),
477  "reco-nonHadPlaneFrac", "ReMId Non-hadronic plane fraction. " );
478 
479  variables.emplace_back(
480  HistAxis("ReMId", kRemidBinning, kRemID),
481  "reco-remid", "ReMId kNN score. " );
482 
483  variables.emplace_back(
484  HistAxis("Number of Uncalibrated Hits in Slice", Binning::Simple(50, 0, 50), kSlcUnCalibNHit),
485  "reco-slcUnCalibNHit", "Uncalibrated slice nhit. " );
486 
487  variables.emplace_back(
488  HistAxis("Slice Mean Time [#mus]", Binning::Simple(60, -50, 550), kSlcMeanTime),
489  "reco-slcMeanTime", "Slice mean time TNS in microseconds. " );
490 
491  variables.emplace_back(
492  HistAxis("Slice Start Time [#mus]", Binning::Simple(60, -50, 550), kSlcStartTime),
493  "reco-slcStartTime", "Slice start time TNS in microseconds. " );
494 
495  variables.emplace_back(
496  HistAxis("Slice End Time [#mus]", Binning::Simple(60, -50, 550), kSlcEndTime),
497  "reco-slcEndTime", "Slice end time TNS in microseconds. " );
498 
499  if(isFD){
500  variables.emplace_back(
501  HistAxis("Slice Minimum X [m]", kFDXYBins, kSlcMinX),
502  "reco-slcMinX", "Slice Min X [m]. " );
503 
504  variables.emplace_back(
505  HistAxis("Slice Minimum Y [m]", kFDXYBins, kSlcMinY),
506  "reco-slcMinY", "Slice Min Y [m]. " );
507 
508  variables.emplace_back(
509  HistAxis("Slice Minimum Z [m]", kFDZBins, kSlcMinZ),
510  "reco-slcMinZ", "Slice Min Z [m]. " );
511 
512  variables.emplace_back(
513  HistAxis("Slice Maximum X [m]", kFDXYBins, kSlcMaxX),
514  "reco-slcMaxX", "Slice Max X [m]. " );
515 
516  variables.emplace_back(
517  HistAxis("Slice Maximum Y [m]", kFDXYBins, kSlcMaxY),
518  "reco-slcMaxY", "Slice Max Y [m]. " );
519 
520  variables.emplace_back(
521  HistAxis("Slice Maximum Z [m]", kFDZBins, kSlcMaxZ),
522  "reco-slcMaxZ", "Slice Max Z [m]. " );
523 
524  variables.emplace_back(
525  HistAxis("Slice Extent X [m]", Binning::Simple(55, 0, 16.0), kSlcExtentX),
526  "reco-slcExtentX", "Slice Extent X [m]. " );
527 
528  variables.emplace_back(
529  HistAxis("Slice Extent Y [m]", Binning::Simple(55, 0, 16.0), kSlcExtentY),
530  "reco-slcExtentY", "Slice Extent Y [m]. " );
531 
532  variables.emplace_back(
533  HistAxis("Slice Extent Z [m]", Binning::Simple(55, 0, 60.0), kSlcExtentZ),
534  "reco-slcExtentZ", "Slice Extent Z [m]. " );
535  }
536 
537  else{
538  variables.emplace_back(
539  HistAxis("Slice Minimum X [m]", kXYBins, kSlcMinX),
540  "reco-slcMinX", "Slice Min X [m]. " );
541 
542  variables.emplace_back(
543  HistAxis("Slice Minimum Y [m]", kXYBins, kSlcMinY),
544  "reco-slcMinY", "Slice Min Y [m]. " );
545 
546  variables.emplace_back(
547  HistAxis("Slice Minimum Z [m]", kZBins, kSlcMinZ),
548  "reco-slcMinZ", "Slice Min Z [m]. " );
549 
550  variables.emplace_back(
551  HistAxis("Slice Maximum X [m]", kXYBins, kSlcMaxX),
552  "reco-slcMaxX", "Slice Max X [m]. " );
553 
554  variables.emplace_back(
555  HistAxis("Slice Maximum Y [m]", kXYBins, kSlcMaxY),
556  "reco-slcMaxY", "Slice Max Y [m]. " );
557 
558  variables.emplace_back(
559  HistAxis("Slice Maximum Z [m]", kZBins, kSlcMaxZ),
560  "reco-slcMaxZ", "Slice Max Z [m]. " );
561 
562  variables.emplace_back(
563  HistAxis("Slice Extent X [m]", Binning::Simple(55, 0, 4.4), kSlcExtentX),
564  "reco-slcExtentX", "Slice Extent X [m]. " );
565 
566  variables.emplace_back(
567  HistAxis("Slice Extent Y [m]", Binning::Simple(55, 0, 4.4), kSlcExtentY),
568  "reco-slcExtentY", "Slice Extent Y [m]. " );
569 
570  variables.emplace_back(
571  HistAxis("Slice Extent Z [m]", Binning::Simple(18, 0, 16), kSlcExtentZ),
572  "reco-slcExtentZ", "Slice Extent Z [m]. " );
573  }
574 
575  variables.emplace_back(
576  HistAxis("Kalman Track Cos #theta_{X}", Binning::Simple(100, -1, 1), kDirX),
577  "reco-dirX", "X-direction of muon track. " );
578 
579  variables.emplace_back(
580  HistAxis("Kalman Track Cos #theta_{Y}", Binning::Simple(100, -1, 1), kDirY),
581  "reco-dirY", "Y-direction of muon track. " );
582 
583  variables.emplace_back(
584  HistAxis("Kalman Track Cos #theta_{Z}", Binning::Simple(50, 0, 1), kDirZ),
585  "reco-dirZ", "Z-direction of muon track. " );
586 
587  variables.emplace_back(
588  HistAxis("Kalman Track Cos #theta_{NuMI}", Binning::Simple(50, 0, 1), kCosNumi),
589  "reco-cosNumi", "Beam direction of muon track. " );
590 
591  variables.emplace_back(
592  HistAxis("Number of Missing Planes in Primary Track", Binning::Simple(50, 0, 50),
593  kTrkNPlaneGap),
594  "reco-trkNPlaneGap", "Track N Plane Gap. " );
595 
596  variables.emplace_back(
597  HistAxis("Visible Slice Energy Per Hit [GeV]", Binning::Simple(40, 0, 0.04),
599  "reco-slcCalEPerNHit", "Slice Energy Per Slice NHit. " );
600 
601  /////////////////////////////////////////////////////
602  std::vector<TangibleAxis> variablesForQePId;
603 
604  variablesForQePId.emplace_back(
605  HistAxis("QePId", Binning::Simple(51, -0.01, 1.01), kQePId),
606  "reco-qepid", "QePId kNN score. " );
607 
608  variablesForQePId.emplace_back(
609  HistAxis("Non QE Neutrino Energy [GeV]", Binning::Simple(50, 0, 5), kNonQeE),
610  "reco-nonQeE", "Non QE Energy. " );
611 
612  variablesForQePId.emplace_back(
613  HistAxis("Hadronic Non-QE Energy [GeV]", Binning::Simple(50, 0, 5), kNonQeHadE),
614  "reco-nonQeHadE", "Hadronic Non QE Energy. " );
615 
616  variablesForQePId.emplace_back(
617  HistAxis("QePId Input: Off-track Energy Ratio", Binning::Simple(100, 0, 1), kQepidOffE),
618  "reco-QepidOffE", "Off-track Energy Ratio, QePId. " );
619 
620  variablesForQePId.emplace_back(
621  HistAxis("QePId Input: Fractional Energy Difference", Binning::Simple(100, -1, 1),
622  kQepidEDiff),
623  "reco-QepidEDiff", "Fractional Energy Difference, QePId. " );
624 
625  variablesForQePId.emplace_back(
626  HistAxis("QePId Input: Fractional Energy Difference Z-test", Binning::Simple(100, -5, 5),
627  kQepidEDiffZ),
628  "reco-QepidEDiffZ", "Fractional Energy Difference Z-test, QePId. " );
629 
630  variablesForQePId.emplace_back(
631  HistAxis("QePId Input: dE/dx Ratio", Binning::Simple(100, 0, 5), kQepidDedx),
632  "reco-QepidDedx", "dE/dx Ratio, applied only to 2 trk events, QePId. " );
633 
634  variablesForQePId.emplace_back(
635  HistAxis("QE Neutrino Energy [GeV]", Binning::Simple(50, 0, 5), kQeE),
636  "reco-qeE", "QE energy. " );
637 
638  variablesForQePId.emplace_back(
639  HistAxis("Hadronic QE Energy [GeV]", Binning::Simple(50, 0, 5), kQeHadE),
640  "reco-qeHadE", "Hadronic energy for QE things. " );
641 
642  variablesForQePId.emplace_back(
643  HistAxis("QE Angle Neutrino Energy [GeV]", Binning::Simple(50, 0, 5), kQeAngleE),
644  "reco-qeAngleE", "QE angle formula energy. " );
645 
646  variablesForQePId.emplace_back(
647  HistAxis("Primary Track Energy [GeV]", Binning::Simple(50, 0, 5), kNumuTrackE),
648  "reco-numuTrackE", "Primary track Energy. " );
649 
650  variablesForQePId.emplace_back(
651  HistAxis("Visible Hadronic Energy [GeV]", Binning::Simple(50, 0, 2), kNumuHadVisE),
652  "reco-HadVisE", "Hadronic Energy . " );
653 
654 
655  std::vector<UsefulHist> hists;
656  hists.reserve(selections.size() * variables.size() + selectionsForQePId.size()
657  * variablesForQePId.size());
658 
659  for(const auto& sel:selections){
660  for(const auto& variable:variables){
661  hists.emplace_back(sel, variable, loader);
662  }
663  }
664 
665  for(const auto& sel:selectionsForQePId){
666  for(const auto& variable:variablesForQePId){
667  hists.emplace_back(sel, variable, loader);
668  }
669  }
670 
671  std::cout << "\nrun : --- run loader.\n";
672  loader.Go();
673  std::cout << "\nrun : --- done.\n";
674  double pot = hists[0].fHist.POT();
675 
676 
677  std::cout << "\nrun : --- save output.\n";
678  // TFile inputfile(kInputFileName.c_str(), "READ");
679  TFile outputfile(kOutputFileName.c_str(), "RECREATE");
680 
681  for(const auto& hist:hists){
682  TH1 *temp = hist.fHist.ToTH1(pot);
683  std::string tempName = hist.fName;
684  temp->Write(tempName.c_str());
685  }
686 
687  std::cout << "\nrun : --- save PoT.\n";
688 
689  TH1F *pothist = new TH1F("meta-TotalPOT", "TotalPOT", 1, 0, 1);
690  TH1F *evthist = new TH1F("meta-TotalEvents", "TotalEvents", 1, 0, 1);
691  pothist->SetBinContent(1, pot);
692  evthist->SetBinContent(1, spillHist->GetEntries());
693  pothist->Write("meta-TotalPOT");
694  evthist->Write("meta-TotalEvents");
695 
696  std::cout << "\nrun : --- close output files.\n";
697  // inputfile.Close();
698  outputfile.Write();
699  outputfile.Close();
700 
701  std::cout << "\nrun : --- done.\n";
702 }
703 
704 #endif
const Var kHadE
Definition: NumuVars.h:23
const Var kQepidEDiff
Det_t det
Detector, ND = 1, FD = 2, NDOS = 3.
Definition: SRHeader.h:28
Near Detector underground.
Definition: SREnums.h:10
const Var kHadNHit([](const caf::SRProxy *sr){unsigned int nought=0;if(sr->trk.kalman.ntracks< 1) return nought;return sr->slc.nhit-sr->trk.kalman.tracks[0].nhit;})
Definition: NumuVars.h:61
const Var kNKalman
Definition: NumuVars.cxx:540
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Far Detector at Ash River.
Definition: SREnums.h:11
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 kReMIdScatLLH
Definition: NumuVars.cxx:555
SRHeader hdr
Header branch: run, subrun, etc.
const Var kSlcMeanTime([](const caf::SRProxy *sr){return sr->slc.meantime/1000.;})
Definition: NumuVars.h:83
const Var kNumuHadTrkE
Definition: NumuVars.cxx:539
const Binning kRemidBinning
Binning for plotting remid attractively.
Definition: Binning.cxx:80
const Var kQepidDedx
void caf_numu_validation(std::string kInputFileName, std::string kOutputFileName, bool isFD)
const Var kSlcMaxZ([](const caf::SRProxy *sr){return sr->slc.boxmax.Z()/100.;})
Definition: NumuVars.h:93
const Var kSlcCalEPerNHit([](const caf::SRProxy *sr){if(sr->slc.nhit > 0) return sr->slc.calE/(1.78 *sr->slc.nhit);return-5.;})
Definition: NumuVars.h:101
const Var kSlcMaxY([](const caf::SRProxy *sr){return sr->slc.boxmax.Y()/100.;})
Definition: NumuVars.h:92
SRVector3D boxmax
Maximum coordinates box containing all the hits [cm].
Definition: SRSlice.h:47
float Y() const
Definition: SRVector3D.h:33
float starttime
start time [ns]
Definition: SRSlice.h:39
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
void Preliminary()
Put NOvA Preliminary on plots.
Tangible< Cut > Selection
Definition: DataMCPair.h:52
const Cut kNumuCosmicRej([](const caf::SRProxy *sr){return(sr->sel.cosrej.anglekal > 0.5 && sr->sel.cosrej.numucontpid2019 > 0.535 && sr->slc.nhit< 400);})
Definition: NumuCuts.h:32
const Var kSliceDuration([](const caf::SRProxy *sr){return(sr->slc.endtime-sr->slc.starttime);})
Definition: NumuVars.h:35
const Cut kNumuContainFD([](const caf::SRProxy *sr){ std::pair< int, int > planes=calcFirstLastLivePlane(sr->slc.firstplane, std::bitset< 14 >(sr->hdr.dibmask));int planestofront=sr->slc.firstplane-planes.first;int planestoback=planes.second-sr->slc.lastplane;return( sr->slc.ncellsfromedge > 1 &&planestofront > 1 &&planestoback > 1 &&sr->sel.contain.kalfwdcell > 10 &&sr->sel.contain.kalbakcell > 10 &&sr->sel.contain.cosfwdcell > 0 &&sr->sel.contain.cosbakcell > 0);})
Definition: NumuCuts.h:20
const Var kSlcExtentY([](const caf::SRProxy *sr){return(sr->slc.boxmax.Y()-sr->slc.boxmin.Y())/100.;})
Definition: NumuVars.h:96
void SetSpillCut(const SpillCut &cut)
const Var kTrkEPerNHit([](const caf::SRProxy *sr){return(sr->trk.kalman.tracks[0].calE/sr->trk.kalman.tracks[0].nhit);})
TString hists[nhists]
Definition: bdt_com.C:3
const Var kReMIdMeasFrac
Definition: NumuVars.cxx:557
const Var kDirZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].dir.Z();})
Definition: NumuVars.h:39
const Var kTrkLength([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].len/100;})
Definition: NumuVars.h:65
const Var kSlcMinX([](const caf::SRProxy *sr){return sr->slc.boxmin.X()/100.;})
Definition: NumuVars.h:87
const Var kHadEPerNHit([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 0.0f;int nHit=sr->slc.nhit-sr->trk.kalman.tracks[0].nhit;if(nHit<=0) return 0.0f;float hadE=sr->energy.numu.hadcalE;return hadE/nHit;})
Definition: NumuVars.h:63
const Var kNumuHadVisE([](const caf::SRProxy *sr){return kNumuHadCalE(sr)+kNumuHadTrkE(sr);})
Definition: NumuVars.h:124
const Var kSlcMaxX([](const caf::SRProxy *sr){return sr->slc.boxmax.X()/100.;})
Definition: NumuVars.h:91
const Cut kNumuContainND([](const caf::SRProxy *sr){return( sr->trk.kalman.ntracks > sr->trk.kalman.idxremid &&sr->slc.ncellsfromedge > 1 &&sr->slc.firstplane > 1 &&sr->slc.lastplane< 212 &&sr->trk.kalman.tracks[0].start.Z()< 1150 &&( sr->trk.kalman.tracks[0].stop.Z()< 1275 ||sr->sel.contain.kalyposattrans< 55) &&( sr->energy.numu.ndhadcalcatE +sr->energy.numu.ndhadcaltranE)< 0.03 &&sr->sel.contain.kalfwdcellnd > 4 &&sr->sel.contain.kalbakcellnd > 8);})
Definition: NumuCuts.h:22
const TVector3 beamDirFD
const Var kDirY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].dir.Y();})
Definition: NumuVars.h:38
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const Var kNumuHadCalE
Definition: NumuVars.cxx:538
const Var kSlcMinY([](const caf::SRProxy *sr){return sr->slc.boxmin.Y()/100.;})
Definition: NumuVars.h:88
const Cut kIsDytmanMEC
Definition: TruthCuts.cxx:187
const Var kQepidEDiffZ
float X() const
Definition: SRVector3D.h:32
Track finder for cosmic rays.
float endtime
end time [ns]
Definition: SRSlice.h:40
SRKalman kalman
Tracks produced by KalmanTrack.
Definition: SRTrackBranch.h:24
const Var kQepidOffE
SRVector3D boxmin
Minimum coordinates box containing all the hits [cm].
Definition: SRSlice.h:46
const Var kQePId
const Var kRemID
PID
Definition: Vars.cxx:81
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.
const Var kCaloE
Summed calorimetric energy of all hits in slice, uncorrected.
Definition: Vars.cxx:52
Tangible(const T &obj, const std::string &shortName, const std::string &blurb)
const Binning kHadronicEnergyBinning
const Var kNHit
Definition: Vars.cxx:71
UsefulHist(Selection sel, TangibleAxis tanAxis, SpectrumLoader &loader, const Cut &bkg=kNoCut)
const Var kTrkStartZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].start.Z()/100;})
Definition: NumuVars.h:53
SRRemid remid
Output from RecoMuonID (ReMId) package.
Definition: SRIDBranch.h:39
Tangible< HistAxis > TangibleAxis
Definition: DataMCPair.h:53
float calE
Calorimetric energy of the cluster [GeV].
Definition: SRSlice.h:38
const Var kCCE
Definition: NumuVars.h:21
const Binning kEnergyBinning
double energy
Definition: plottest35.C:25
const Var kNonQeHadE([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999) return float(sr->energy.numu.recotrkcchadE);return-5.f;})
Definition: NumuVars.h:77
#define pot
const Var kSlcStartTime([](const caf::SRProxy *sr){return sr->slc.starttime/1000.;})
Definition: NumuVars.h:84
const Binning kXYBins
Definition: VarsAndCuts.h:95
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
const Cut kNumuNCRej([](const caf::SRProxy *sr){return(sr->sel.remid.pid >0.75);})
Definition: NumuCuts.h:24
loader
Definition: demo0.py:10
short nnu
Number of neutrinos in nu vector (0 or 1)
Definition: SRTruthBranch.h:37
const Var kQeE
Definition: NumuVars.cxx:562
float Z() const
Definition: SRVector3D.h:34
unsigned int nhit
number of hits
Definition: SRSlice.h:22
TLatex * prelim
Definition: Xsec_final.C:133
const Var kQeHadE([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999) return sr->energy.numu.trkqeE-sr->energy.numu.recomuonE;return-5.f;})
Definition: NumuVars.h:79
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
const Var kTrkEndZ([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.Z()/100;})
Definition: NumuVars.h:57
OStream cout
Definition: OStream.cxx:6
const TVector3 beamDirND
std::string fBlurb
Definition: DataMCPair.h:46
unsigned int ncalhit
number of hits with calibration
Definition: SRSlice.h:23
const Var kQeAngleE
Definition: NumuVars.cxx:563
The StandardRecord is the primary top-level object in the Common Analysis File trees.
const Var kTrkNPlaneGap([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999) return int(sr->trk.kalman.tracks[0].nplanegap);return 500;})
Definition: NumuVars.h:99
const Var kNumuTrackE({"energy.mutrkE","trk.nkalman","sel.remid.bestidx"}, [](const caf::StandardRecord *sr){if(sr->trk.nkalman > 0 &&sr->sel.remid.bestidx!=999) return sr->energy.mutrkE[sr->sel.remid.bestidx].E;return-5.f;})
const Var kTrkEndY([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.Y()/100;})
Definition: NumuVars.h:56
const Var kTrkEndX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-10.0f;return sr->trk.kalman.tracks[0].stop.X()/100;})
Definition: NumuVars.h:55
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kSlcEndTime([](const caf::SRProxy *sr){return sr->slc.endtime/1000.;})
Definition: NumuVars.h:85
const Var kDirX([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return-5.f;return sr->trk.kalman.tracks[0].dir.X();})
Definition: NumuVars.h:37
#define SIMPLEVAR(CAFNAME)
For Vars where literally all you need is a single CAF variable.
Definition: Var.h:88
SRIDBranch sel
Selector (PID) branch.
assert(nhit_max >=nhit_nbins)
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
SRSlice slc
Slice branch: nhit, extents, time, etc.
TVector3 NuMIBeamDirection(caf::Det_t det)
Average direction of NuMI neutrinos in a given detector This function is a copy of geo::GeometryBase:...
Definition: Utilities.cxx:329
SRTruthBranch mc
Truth branch for MC: energy, flavor, etc.
const Var kReMIdDEDxLLH
Definition: NumuVars.cxx:556
string tempName
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
const Var kNonQeE
Definition: NumuVars.cxx:561
double T
Definition: Xdiff_gwt.C:5
const Cut kNumuQuality
Definition: NumuCuts.h:18
const Var kSlcExtentX([](const caf::SRProxy *sr){return(sr->slc.boxmax.X()-sr->slc.boxmin.X())/100.;})
Definition: NumuVars.h:95
SRTrackBranch trk
Track branch: nhit, len, etc.
SREnergyBranch energy
Energy estimator branch.
const Var kSlcExtentZ([](const caf::SRProxy *sr){return(sr->slc.boxmax.Z()-sr->slc.boxmin.Z())/100.;})
Definition: NumuVars.h:97
const Binning kZBins
Definition: VarsAndCuts.h:96
const Var kSlcUnCalibNHit([](const caf::SRProxy *sr){return sr->slc.nhit-sr->slc.ncalhit;})
Definition: NumuVars.h:81
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
const Var kCosNumi([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks > 0 &&sr->trk.kalman.idxremid!=999){if(sr->hdr.det==1){return sr->trk.kalman.tracks[0].dir.Dot(beamDirND);}if(sr->hdr.det==2){return sr->trk.kalman.tracks[0].dir.Dot(beamDirFD);}}return-5.f;})
Definition: NumuVars.h:43
const Var kSlcMinZ([](const caf::SRProxy *sr){return sr->slc.boxmin.Z()/100.;})
Definition: NumuVars.h:89
const Var kTrkNhits([](const caf::SRProxy *sr){if(sr->trk.kalman.ntracks< 1) return 65535;return int(sr->trk.kalman.tracks[0].nhit);})
Definition: NumuVars.h:59
float meantime
mean time, weighted by charge [ns]
Definition: SRSlice.h:41
enum BeamMode string