SlicerAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file SlicerAna_module.cc
3 // \brief Analysis module to validate slicer performance
4 // \version $Id: SlicerAna.cxx,v 1.4 2012-12-19 14:51:38 mbaird42 Exp $
5 // \author $Author: mbaird42 $
6 ////////////////////////////////////////////////////////////////////////
7 
8 // C/C++ includes
9 #include <cmath>
10 #include <iostream>
11 #include <map>
12 #include <vector>
13 
14 // ROOT includes
15 #include "TTree.h"
16 #include "TH1F.h"
17 #include "TH2F.h"
18 
19 // Framework includes
30 #include "fhiclcpp/ParameterSet.h"
32 
33 // NOvASoft includes
35 #include "RecoBase/CellHit.h"
36 #include "RecoBase/Cluster.h"
37 #include "Utilities/AssociationUtil.h"
41 
42 #include "Calibrator/Calibrator.h"
43 #include "MCCheater/BackTracker.h"
44 
45 // Required for CAFAna Cuts
47 #include "CAFAna/Cuts/Cuts.h"
48 #include "CAFAna/Cuts/TruthCuts.h"
49 #include "CAFAna/Cuts/TimingCuts.h"
50 #include "CAFAna/Vars/Vars.h"
51 
54 
55 
56 #include "SummaryData/POTSum.h"
57 
58 // SlicerAna header
59 namespace slicer {
60  class SlicerAna : public art::EDAnalyzer {
61  public:
62  explicit SlicerAna(fhicl::ParameterSet const& pset);
63  ~SlicerAna();
64 
65  void analyze(const art::Event& evt);
66  void reconfigure(const fhicl::ParameterSet& pset);
67  void beginJob();
68  void endJob();
69  virtual void endSubRun(const art::SubRun& sr);
70 
71  std::vector<cheat::NeutrinoEffPur> sortEffPur(const std::vector<cheat::NeutrinoEffPur>& nuEP);
72 
73  private:
74  void ClearVectors();
75 
76  std::string fSlicerLabel; ///< label for the process that made the slices
77  bool fIsMC; ///< is the file being processed MC?
78  std::string fCAFLabel; ///< label for the process that made the standard records
79  bool fApplyCAFCuts; ///< should we apply CAF cuts?
80  int fCutType; ///< what cuts should we apply (see CAFCutter.h)
81  float fTotalPOT; ///< what is the POT in this file?
82 
83  recovalid::CAFCutter fCAFCutter; ///< the CAFCutter helper class
84 
85  TTree* fOutTree;
86 
87 
88 
89  // Histograms created for production validation
90  TH1F *fNumNu;
91  TH1F *fNumMCT;
92  TH1F *fNumSlices;
94  TH1F *fNumNoise;
95  TH1F *fNumHits;
97 
98  TH1F *fComp;
99  TH1F *fPur;
100  TH2F *fCompPur;
101  TH1F *fCompNoise;
102  TH1F *fPurNoise;
103 
104  TH1F *fNSvsVE;
105  TH1F *fNGvsVE;
106  TH1F *fEffVsVisE;
107 
108 
109  TH2F *fPCXallhits;
110  TH2F *fPCYallhits;
111  TH1F *fTallhits;
112  TH2F *fPCXave;
113  TH2F *fPCYave;
114  TH1F *fTave;
115  TH2F *fPCXdiff;
116  TH2F *fPCYdiff;
117  TH1F *fTsd;
118  TH1F *fSlDur;
119 
125  TH1F *fTaveNoise;
126 
128 
129  TH1F *fTres;
130  TH1F *fTresNoise;
131 
132  TH1F *fPOT;
133 
134  //
135  // info to go into TTree
136  //
137 
138  // event level info
139  int run; ///< run number
140  int subrun; ///< subrun number
141  int event; ///< event number
142 
143  int NumNoise; ///< number of noise slices
144  int NumSlice; ///< number of non-noise slices
145 
146  int NumNuFLSHitALL; ///< number of nu counted by FLSHits
147  int NumNuCellHitALL; ///< number of nu counted by CellHits
148  int NumNuCellHit; ///< number of nu in non-noise slices counted by CellHits
149 
150  int NumMCTFLSHitALL; ///< total number of MCTruths counted by FLSHits
151  int NumMCTCellHitALL; ///< total number of MCTruths counted by CellHits
152  int NumMCTCellHit; ///< total number of MCTruths in non-noise slices counted by CellHits
153 
154  // slice level info
155  std::vector<unsigned int> IsNoiseVec; ///< 0 = non-noise slice, 1 = noise slice
156  std::vector<int> NumHitsXVec; ///< number of X hits per slice
157  std::vector<int> NumHitsYVec; ///< number of Y hits per slice
158  std::vector<double> slEVec; ///< slice.TotalGeV()
159  std::vector<float> slPEVec; ///< sum of PE for all hits in the slice
160 
161  // geometric info about all slices
162  std::vector<double> PaveVec; ///< average plane number for hits in a non-noise slice
163  std::vector<double> TaveVec; ///< average time for hits in a non-noise slice
164  std::vector<double> CXaveVec; ///< average cell number for X hits in a non-noise slice
165  std::vector<double> CYaveVec; ///< average cell number for Y hits in a non-noise slice
166 
167  std::vector<double> PsdVec; ///< st. dev. of plane number for hits in a non-noise slice
168  std::vector<double> TsdVec; ///< st. dev. time for hits in a non-noise slice
169  std::vector<double> CXsdVec; ///< st. dev. cell number for X hits in a non-noise slice
170  std::vector<double> CYsdVec; ///< st. dev. cell number for Y hits in a non-noise slice
171 
172  std::vector<double> PhiVec; ///< highest plane in the slice
173  std::vector<double> PloVec; ///< lowest plane in the slice
174  std::vector<double> ThiVec; ///< highest time in the slice
175  std::vector<double> TloVec; ///< lowest plane in the slice
176  std::vector<double> CXhiVec; ///< highest X cell in the slice
177  std::vector<double> CXloVec; ///< lowest X cell in the slice
178  std::vector<double> CYhiVec; ///< highest Y cell in the slice
179  std::vector<double> CYloVec; ///< lowest Y cell in the slice
180 
181  // info about slicer performance (matching to MCTruths...)
182  std::vector<float> nu1XVec; ///< x location for primary nu associated with slice
183  std::vector<float> nu1YVec; ///< y location for primary nu associated with slice
184  std::vector<float> nu1ZVec; ///< z location for primary nu associated with slice
185  std::vector<float> nu1EVec; ///< Energy for primary nu associated with slice
186 
187  std::vector<float> nu2XVec; ///< x location for second nu associated with slice
188  std::vector<float> nu2YVec; ///< y location for second nu associated with slice
189  std::vector<float> nu2ZVec; ///< z location for second nu associated with slice
190  std::vector<float> nu2EVec; ///< Energy for second nu associated with slice
191 
192  std::vector<float> nu3XVec; ///< x location for third nu associated with slice
193  std::vector<float> nu3YVec; ///< y location for third nu associated with slice
194  std::vector<float> nu3ZVec; ///< z location for third nu associated with slice
195  std::vector<float> nu3EVec; ///< Energy for third nu associated with slice
196 
197  std::vector<float> nu4XVec; ///< x location for fourth nu associated with slice
198  std::vector<float> nu4YVec; ///< y location for fourth nu associated with slice
199  std::vector<float> nu4ZVec; ///< z location for fourth nu associated with slice
200  std::vector<float> nu4EVec; ///< Energy for fourth nu associated with slice
201 
202  std::vector<float> nu5XVec; ///< x location for fifth nu associated with slice
203  std::vector<float> nu5YVec; ///< y location for fifth nu associated with slice
204  std::vector<float> nu5ZVec; ///< z location for fifth nu associated with slice
205  std::vector<float> nu5EVec; ///< Energy for fifth nu associated with slice
206 
207  std::vector<int> nu1pdgVec; ///< PDG code for primary nu associated with slice
208  std::vector<int> nu1CCNCVec; ///< CCNC for primary nu associated with slice
209  std::vector<int> nu1IntTypeVec; ///< Interaction type for primary neutrino associated with slice
210  std::vector<double> nu1VisEslVec; ///< visible energy for the best matched nu in the slice
211  std::vector<double> nu1VisEtotVec; ///< total visible energy for the best matched nu in the entire event
212 
213  std::vector<int> nu2pdgVec; ///< PDG code for second nu associated with slice
214  std::vector<int> nu2CCNCVec; ///< CCNC for second nu associated with slice
215  std::vector<int> nu2IntTypeVec; ///< Interaction type for second neutrino associated with slice
216  std::vector<double> nu2VisEslVec; ///< visible energy for the second best matched nu in the slice
217  std::vector<double> nu2VisEtotVec; ///< total visible energy for the second best matched nu in the entire event
218 
219  std::vector<int> nu3pdgVec; ///< PDG code for third nu associated with slice
220  std::vector<int> nu3CCNCVec; ///< CCNC for third nu associated with slice
221  std::vector<int> nu3IntTypeVec; ///< Interaction type for third neutrino associated with slice
222  std::vector<double> nu3VisEslVec; ///< visible energy for the third best matched nu in the slice
223  std::vector<double> nu3VisEtotVec; ///< total visible energy for the third best matched nu in the entire event
224 
225  std::vector<int> nu4pdgVec; ///< PDG code for fourth nu associated with slice
226  std::vector<int> nu4CCNCVec; ///< CCNC for fourth nu associated with slice
227  std::vector<int> nu4IntTypeVec; ///< Interaction type for fourth neutrino associated with slice
228  std::vector<double> nu4VisEslVec; ///< visible energy for the fourth best matched nu in the slice
229  std::vector<double> nu4VisEtotVec; ///< total visible energy for the fourth best matched nu in the entire event
230 
231  std::vector<int> nu5pdgVec; ///< PDG code for fifth nu associated with slice
232  std::vector<int> nu5CCNCVec; ///< CCNC for fifth nu associated with slice
233  std::vector<int> nu5IntTypeVec; ///< Interaction type for fifth neutrino associated with slice
234  std::vector<double> nu5VisEslVec; ///< visible energy for the fifth best matched nu in the slice
235  std::vector<double> nu5VisEtotVec; ///< total visible energy for the fifth best matched nu in the entire event
236 
237  std::vector<int> pdgPartMostVec; ///< the PDG code for the particle that contributed the most energy to this slice
238  std::vector<int> NumMatchesVec; ///< number of MCTruth matches per slice
239 
240  std::vector<double> nu1EffVec; ///< efficiency for primary nu associated with slice
241  std::vector<double> nu1PurVec; ///< purity for primary nu associated with slice
242  std::vector<double> nu2EffVec; ///< efficiency for second nu associated with slice
243  std::vector<double> nu2PurVec; ///< purity for second nu associated with slice
244  std::vector<double> nu3EffVec; ///< efficiency for third nu associated with slice
245  std::vector<double> nu3PurVec; ///< purity for third nu associated with slice
246  std::vector<double> nu4EffVec; ///< efficiency for fourth nu associated with slice
247  std::vector<double> nu4PurVec; ///< purity for fourth nu associated with slice
248  std::vector<double> nu5EffVec; ///< efficiency for fifth nu associated with slice
249  std::vector<double> nu5PurVec; ///< purity for fifth nu associated with slice
250 
251  std::vector<unsigned int> IsMichelPrimaryVec; ///< is this slice a Michel e? (the primary contribution by energy was from a michel e)
252  std::vector<unsigned int> IsMichelPresentVec; ///< did a Michel e contribute to this slice?
253  std::vector<unsigned int> IsLateSliceVec; ///< is the average time for this slice > 200 ns after the time of the neutrino interaction?
254 
255  std::vector<double> PercentNoiseVec; ///< Percentage of slice that is true noise (by hit)
256  std::vector<double> PercentNoiseByPEVec; ///< Percentage of slice that is true noise (weighted by PE)
257 
258  // MCTruth level info
259 
260  // the main F.O.M. number
261  unsigned int NgoodSL;
262 
263  };
264 }
265 
266 
267 
268 namespace slicer
269 {
270  //.......................................................................
272  : EDAnalyzer(pset)
273  {
274  this->reconfigure(pset);
275  }
276 
277  //......................................................................
279  {
280  //======================================================================
281  // Clean up any memory allocated by your module
282  //======================================================================
283  }
284 
285  //......................................................................
287  {
288  fSlicerLabel = pset.get<std::string> ("SlicerLabel");
289  fIsMC = pset.get<bool> ("IsMC");
290  fCAFLabel = pset.get<std::string> ("CAFLabel");
291  fApplyCAFCuts = pset.get<bool> ("ApplyCAFCuts");
292  fCutType = pset.get<int> ("CutType");
293  }
294 
295  //......................................................................
297  {
298  fTotalPOT = 0.0;
299  NgoodSL = 0;
300 
302 
303  // define the TTree and branch variables
304  fOutTree = tfs->make<TTree>("SlicerAna", "Slicer Analysis Tree");
305  fOutTree->Branch("run", &run);
306  fOutTree->Branch("subrun", &subrun);
307  fOutTree->Branch("event", &event);
308 
309  fOutTree->Branch("NumNoise", &NumNoise);
310  fOutTree->Branch("NumSlice", &NumSlice);
311  fOutTree->Branch("NumNuFLSHitALL", &NumNuFLSHitALL);
312  fOutTree->Branch("NumNuCellHitALL", &NumNuCellHitALL);
313  fOutTree->Branch("NumNuCellHit", &NumNuCellHit);
314  fOutTree->Branch("NumMCTFLSHitALL", &NumMCTFLSHitALL);
315  fOutTree->Branch("NumMCTCellHitALL", &NumMCTCellHitALL);
316  fOutTree->Branch("NumMCTCellHit", &NumMCTCellHit);
317 
318  fOutTree->Branch("IsNoise", &IsNoiseVec);
319  fOutTree->Branch("NumHitsX", &NumHitsXVec);
320  fOutTree->Branch("NumHitsY", &NumHitsYVec);
321  fOutTree->Branch("slE", &slEVec);
322  fOutTree->Branch("slPE", &slPEVec);
323 
324  fOutTree->Branch("Pave", &PaveVec);
325  fOutTree->Branch("Tave", &TaveVec);
326  fOutTree->Branch("CXave", &CXaveVec);
327  fOutTree->Branch("CYave", &CYaveVec);
328 
329  fOutTree->Branch("Psd", &PsdVec);
330  fOutTree->Branch("Tsd", &TsdVec);
331  fOutTree->Branch("CXsd", &CXsdVec);
332  fOutTree->Branch("CYsd", &CYsdVec);
333 
334  fOutTree->Branch("Phi", &PhiVec);
335  fOutTree->Branch("Thi", &ThiVec);
336  fOutTree->Branch("CXhi", &CXhiVec);
337  fOutTree->Branch("CYhi", &CYhiVec);
338 
339  fOutTree->Branch("Plo", &PloVec);
340  fOutTree->Branch("Tlo", &TloVec);
341  fOutTree->Branch("CXlo", &CXloVec);
342  fOutTree->Branch("CYlo", &CYloVec);
343 
344  fOutTree->Branch("nu1X", &nu1XVec);
345  fOutTree->Branch("nu1Y", &nu1YVec);
346  fOutTree->Branch("nu1Z", &nu1ZVec);
347  fOutTree->Branch("nu1E", &nu1EVec);
348  fOutTree->Branch("nu2X", &nu2XVec);
349  fOutTree->Branch("nu2Y", &nu2YVec);
350  fOutTree->Branch("nu2Z", &nu2ZVec);
351  fOutTree->Branch("nu2E", &nu2EVec);
352  fOutTree->Branch("nu3X", &nu3XVec);
353  fOutTree->Branch("nu3Y", &nu3YVec);
354  fOutTree->Branch("nu3Z", &nu3ZVec);
355  fOutTree->Branch("nu3E", &nu3EVec);
356  fOutTree->Branch("nu4X", &nu4XVec);
357  fOutTree->Branch("nu4Y", &nu4YVec);
358  fOutTree->Branch("nu4Z", &nu4ZVec);
359  fOutTree->Branch("nu4E", &nu4EVec);
360  fOutTree->Branch("nu5X", &nu5XVec);
361  fOutTree->Branch("nu5Y", &nu5YVec);
362  fOutTree->Branch("nu5Z", &nu5ZVec);
363  fOutTree->Branch("nu5E", &nu5EVec);
364 
365  fOutTree->Branch("nu1pdg", &nu1pdgVec);
366  fOutTree->Branch("nu1CCNC", &nu1CCNCVec);
367  fOutTree->Branch("nu1IntType", &nu1IntTypeVec);
368  fOutTree->Branch("nu1VisEsl", &nu1VisEslVec);
369  fOutTree->Branch("nu1VisEtot", &nu1VisEtotVec);
370  fOutTree->Branch("nu2pdg", &nu2pdgVec);
371  fOutTree->Branch("nu2CCNC", &nu2CCNCVec);
372  fOutTree->Branch("nu2IntType", &nu2IntTypeVec);
373  fOutTree->Branch("nu2VisEsl", &nu2VisEslVec);
374  fOutTree->Branch("nu2VisEtot", &nu2VisEtotVec);
375  fOutTree->Branch("nu3pdg", &nu3pdgVec);
376  fOutTree->Branch("nu3CCNC", &nu3CCNCVec);
377  fOutTree->Branch("nu3IntType", &nu3IntTypeVec);
378  fOutTree->Branch("nu3VisEsl", &nu3VisEslVec);
379  fOutTree->Branch("nu3VisEtot", &nu3VisEtotVec);
380  fOutTree->Branch("nu4pdg", &nu4pdgVec);
381  fOutTree->Branch("nu4CCNC", &nu4CCNCVec);
382  fOutTree->Branch("nu4IntType", &nu4IntTypeVec);
383  fOutTree->Branch("nu4VisEsl", &nu4VisEslVec);
384  fOutTree->Branch("nu4VisEtot", &nu4VisEtotVec);
385  fOutTree->Branch("nu5pdg", &nu5pdgVec);
386  fOutTree->Branch("nu5CCNC", &nu5CCNCVec);
387  fOutTree->Branch("nu5IntType", &nu5IntTypeVec);
388  fOutTree->Branch("nu5VisEsl", &nu5VisEslVec);
389  fOutTree->Branch("nu5VisEtot", &nu5VisEtotVec);
390 
391  fOutTree->Branch("pdgPartMost", &pdgPartMostVec);
392  fOutTree->Branch("NumMatches", &NumMatchesVec);
393 
394  fOutTree->Branch("nu1Eff", &nu1EffVec);
395  fOutTree->Branch("nu1Pur", &nu1PurVec);
396  fOutTree->Branch("nu2Eff", &nu2EffVec);
397  fOutTree->Branch("nu2Pur", &nu2PurVec);
398  fOutTree->Branch("nu3Eff", &nu3EffVec);
399  fOutTree->Branch("nu3Pur", &nu3PurVec);
400  fOutTree->Branch("nu4Eff", &nu4EffVec);
401  fOutTree->Branch("nu4Pur", &nu4PurVec);
402  fOutTree->Branch("nu5Eff", &nu5EffVec);
403  fOutTree->Branch("nu5Pur", &nu5PurVec);
404 
405  fOutTree->Branch("IsMichelPrimary", &IsMichelPrimaryVec);
406  fOutTree->Branch("IsMichelPresent", &IsMichelPresentVec);
407  fOutTree->Branch("IsLateSlice", &IsLateSliceVec);
408 
409  fOutTree->Branch("PercentNoise", &PercentNoiseVec);
410  fOutTree->Branch("PercentNoiseByPE", &PercentNoiseByPEVec);
411 
412 
413 
414  // book the histograms
415  fNumNu = tfs->make<TH1F>("NumNu",
416  "# of Nu;# of neutrinos in event;count",
417  101,-0.5,100.5);
418  fNumMCT = tfs->make<TH1F>("NumMCT",
419  "# of MCTruths;# of MCTruth objects in event;count",
420  201,-0.5,200.5);
421  fNumSlices = tfs->make<TH1F>("NumSlices",
422  "# of non-noise slices;# of non-noise slices per event;count",
423  101,-0.5,100.5);
424  fNumMCTSlices = tfs->make<TH2F>("NumMCTSlices",
425  "# of MCTruths that contributed to non-noise slices vs. # of non-noise slices;# of MCTruth objects;# of slices",
426  101,-0.5,100.5,101,-0.5,100.5);
427  fNumNoise = tfs->make<TH1F>("NumNoiseSlices",
428  "# of noise slices;# of noise slices per event;count",
429  21,-0.5,20.5);
430  fNumHits = tfs->make<TH1F>("NumSliceHits",
431  "# of hits per non-noise slice;# of hits per non-noise slice;count",
432  1000,0.0,5000.0);
433  fNumHitsNoise = tfs->make<TH1F>("NumSliceHitsNoise",
434  "# of hits in the noise slice;# of hits per noise slice;count",
435  8000,0.0,40000.0);
436 
437  fComp = tfs->make<TH1F>("SliceCompleteness",
438  "Completeness for most pure MCTruth in slice; slice completeness; count",
439  101,-0.005,1.005);
440  fPur = tfs->make<TH1F>("SlicePurity",
441  "Purity for most pure MCTruth in slice; slice purity; count",
442  101,-0.005,1.005);
443  fCompPur = tfs->make<TH2F>("SliceCompletenessVsPurity",
444  "Completeness/Purity for 'most pure' MCTruth by slice;Slice Purity;Slice Completeness",
445  101,-0.005,1.005,101,-0.005,1.005);
446  fCompNoise = tfs->make<TH1F>("NoiseSliceCompletenss",
447  "Completeness for all MCTruths in NOISE slice; completeness in noise slice; count",
448  101,-0.005,1.005);
449  fPurNoise = tfs->make<TH1F>("NoiseSlicePurity",
450  "Purity for all MCTruths in NOISE slice; purity in noise slice; count",
451  101,-0.005,1.005);
452 
453  fNSvsVE = new TH1F("SliceVisE","Slice VisE;slice visible energy (GeV)",80,0.0,20.0);
454  fNGvsVE = new TH1F("GoodSliceVisE","good slice VisE;good slice (> 90% efficiency, > 90% purity) visible energy (GeV)",80,0.0,20.0);
455  fNSvsVE->Reset();
456  fNGvsVE->Reset();
457  fEffVsVisE = tfs->make<TH1F>("EffVsVisE",
458  "Efficiency vs. visible energy of best matched MCTruth for all non-noise slices;Visible energy [GeV];slice efficiency",
459  80,0.0,20.0);
460 
461  fPCXallhits = tfs->make<TH2F>("SlicePlaneVsCellX",
462  "Plane vs. X cell for all hits in non-noise slices;plane;X cell",
463  1320,-20.0,1300.0,510,-10.0,500.0);
464  fPCYallhits = tfs->make<TH2F>("SlicePlaneVsCellY",
465  "Plane vs. Y cell for all hits in non-noise slices;plane;Y cell",
466  1320,-20.0,1300.0,510,-10.0,500.0);
467  fTallhits = tfs->make<TH1F>("SliceTime",
468  "Time for all hits in non-noise slices; time [ns]; count",
469  12000,-50000.0,550000.0);
470  fPCXave = tfs->make<TH2F>("SlicePlaneCellXAvg",
471  "avg hit plane vs. avg hit X cell;avg slice plane;avg slice X cell",
472  660,-20.0,1300.0,255,-10.0,500.0);
473  fPCYave = tfs->make<TH2F>("SlicePlaneCellYAvg",
474  "ave hit plane vs. ave hit Y cell;avg slice plane;avg slice Y cell",
475  660,-20.0,1300.0,255,-10.0,500.0);
476  fTave = tfs->make<TH1F>("SliceTAvg",
477  "ave hit time per slice; avg slice time [ns]; count",
478  10000,0.0,500000.0);
479  fPCXdiff = tfs->make<TH2F>("SlicePlaneCellXDiff",
480  "diff btwn hi and lo for plane vs X cell per slice;slice delta plane;slicedelta X cell",
481  500,0.0,1000.0,250,0.0,500.0);
482  fPCYdiff = tfs->make<TH2F>("SlicePlaneCellYDiff",
483  "diff btwn hi and lo for plane vs Y cell per slice;slice delta plane;slice delta Y cell",
484  500,0.0,1000.0,250,0.0,500.0);
485  fTsd = tfs->make<TH1F>("SliceTsd",
486  "st. dev. of hit time per slice; slice sigma time [ns]; count",
487  500000,0.0,500000.0);
488  fSlDur = tfs->make<TH1F>("SliceDur",
489  "Slice Duration (last hit time - first hit time);slice deltaT [ns];count",
490  200,0.0,2000.0);
491 
492  fPCXallhitsNoise = tfs->make<TH2F>("NoiseSlicePlaneCellX",
493  "Plane vs. X cell for all hits (noise slice);noise slice plane;noise slice X cell",
494  1320,-20.0,1300.0,510,-10.0,500.0);
495  fPCYallhitsNoise = tfs->make<TH2F>("NoiseSlicePlaneCellY",
496  "Plane vs. Y cell for all hits (noise slice);noise slice plane;noise slice Y cell",
497  1320,-20.0,1300.0,510,-10.0,500.0);
498  fTallhitsNoise = tfs->make<TH1F>("NoiseSliceT",
499  "Time for all hits (noise slice); noise slice hit time [ns]; count",
500  12000,-50000.0,550000.0);
501  fPCXaveNoise = tfs->make<TH2F>("NoiseSlicePlaneCellXAvg",
502  "avg hit plane vs. ave hit X cell (noise slice);avg noise slice plane;avg noise slice X cell",
503  660,-20.0,1300.0,255,-10.0,500.0);
504  fPCYaveNoise = tfs->make<TH2F>("NoiseSlicePlaneCellYAvg",
505  "avg hit plane vs. ave hit Y cell (noise slice);avg noise slice plane;avg noise slice Y cell",
506  660,-20.0,1300.0,255,-10.0,500.0);
507  fTaveNoise = tfs->make<TH1F>("NoiseSliceTAvg",
508  "avg hit time per slice (noise slice);avg noise slice time [ns]; count",
509  10000,0.0,500000.0);
510 
511  fEvtNoiseFraction = tfs->make<TH1F>("EvtNoiseFraction",
512  "Fraction of all hits in the noise slice; noise fraction; count",
513  10000,0.0,1.0001);
514 
515  fTres = tfs->make<TH1F>("SliceHitTimeResolution",
516  "Computed timing resolution for hits in non-noise slices;non-noise hit time resolution [ns];count",
517  500,0.0,500.0);
518  fTresNoise = tfs->make<TH1F>("NoiseHitTimeResolution",
519  "Computed timing resolution for hits in the noise slice;noise hit time resolution [ns];count",
520  500,0.0,500.0);
521  fPOT = tfs->make<TH1F>("fPOT","total POT",1,0,1);
522 
523  }
524 
525  //......................................................................
527  {
528  if(fIsMC) {
529  std::cout << "FOM: " << NgoodSL
530  << "\n========== End Job ==========\n\n";
531 
532  for(int i = 1; i <= fNSvsVE->GetNbinsX(); ++i) {
533  double ratio;
534  if(fNSvsVE->GetBinContent(i) != 0.0) {
535  ratio = (fNGvsVE->GetBinContent(i))/(fNSvsVE->GetBinContent(i));
536  fEffVsVisE->Fill(fNSvsVE->GetBinCenter(i),ratio);
537  }
538  }
539 
540  }
541  fPOT->Fill(0.5,fTotalPOT);
542  }
543 
544  //......................................................................
546  {
548  if(fIsMC) sr.getByLabel("generator",pot);
549  else sr.getByLabel("ifdbspillinfo",pot);
550 
551  if(!pot.failedToGet()) {
552  if(fIsMC) fTotalPOT += pot->totgoodpot;
553  else fTotalPOT += 1.0e12*(pot->totgoodpot);
554  }
555  }
556 
557  //......................................................................
559  {
560 
561  this->ClearVectors();
562  int NSlice = 0;
563  int NNoiseSlice = 0;
564  int NhitNoise = 0;
565  int NhitTotal = 0;
566 
567  // get the MCTruth info
569  if(fIsMC) evt.getByLabel("generator", mclist);
570 
571  // get the slices
573  evt.getByLabel(fSlicerLabel, slices);
574 
575  // Get the FMP for the Standard Records
576  art::FindManyP<caf::StandardRecord> recordFMP(slices, evt, fCAFLabel);
577 
578 
579  // make the needed services
582 
583  art::PtrVector<rb::Cluster> ARTslices;
584  for(unsigned int i = 0; i < slices->size(); ++i)
585  ARTslices.push_back(art::Ptr<rb::Cluster>(slices,i));
586 
587  std::vector< std::vector< cheat::NeutrinoEffPur > > ALLnuEffPur;
588  if(fIsMC)
589  ALLnuEffPur = BT->SlicesToMCTruthsTable(ARTslices);
590 
591  // sets to record info about neutrinos/MCTruths that interact in the detector
592  std::set< art::Ptr< simb::MCTruth > > NuListFLSHitALL;
593  std::set< art::Ptr< simb::MCTruth > > NuListCellHitALL;
594  std::set< art::Ptr< simb::MCTruth > > NuListCellHit;
595 
596  std::set< art::Ptr< simb::MCTruth > > MCTListFLSHitALL;
597  std::set< art::Ptr< simb::MCTruth > > MCTListCellHitALL;
598  std::set< art::Ptr< simb::MCTruth > > MCTListCellHit;
599 
600 
601 
602  run = evt.run();
603  subrun = evt.subRun();
604  event = evt.id().event();
605 
606 
607 
608  //
609  // Loop over all slices and gather info for each slice...
610  //
611 
612  // std::cout << "\n\n" << evt.id() << "\nsize of MClist = " << mclist->size() << "\n\n";
613 
614  for(unsigned int i = 0; i < slices->size(); ++i) {
615 
616  // Need to count noise slices before we make CAF cuts since the noise slice will never pass.
617  if( (*slices)[i].IsNoise() ) {
618  NNoiseSlice++;
619  }
620 
621  // Apply the CAF-level cuts
622  bool pass = true;
623  if( fApplyCAFCuts ) {
624  // get record associated with the slice
625  std::vector< art::Ptr<caf::StandardRecord> > records = recordFMP.at(i);
626  if (records.size() > 0 && recordFMP.isValid()) {
627  pass = fCAFCutter.passCuts(fCutType, records[0].get());
628  }
629  else { pass = false; }
630  }
631 
632  if(!pass) continue;
633 
634  // define default values for the slice info
635  unsigned int IsNoise = 0;
636  int NumHitsX = -1;
637  int NumHitsY = -1;
638  double slE = 0.0;
639  float slPE = 0.0;
640 
641  double Pave = 0.0;
642  double CXave = 0.0;
643  double CYave = 0.0;
644  double Tave = 0.0;
645 
646  double Psd = 0.0;
647  double CXsd = 0.0;
648  double CYsd = 0.0;
649  double Tsd = 0.0;
650 
651  double Phi = -1.0e6;
652  double Plo = 1.0e6;
653  double Thi = -1.0e6;
654  double Tlo = 1.0e6;
655  double CXhi = -1.0e6;
656  double CXlo = 1.0e6;
657  double CYhi = -1.0e6;
658  double CYlo = 1.0e6;
659 
660  float nu1X = 0.0, nu1Y = 0.0, nu1Z = 0.0, nu1E = -1.0;
661  int nu1pdg = 0, nu1CCNC = -1, nu1IntType = 0;
662  double nu1VisEsl = -1.0, nu1VisEtot = -1.0;
663 
664  float nu2X = 0.0, nu2Y = 0.0, nu2Z = 0.0, nu2E = -1.0;
665  int nu2pdg = 0, nu2CCNC = -1, nu2IntType = 0;
666  double nu2VisEsl = -1.0, nu2VisEtot = -1.0;
667 
668  float nu3X = 0.0, nu3Y = 0.0, nu3Z = 0.0, nu3E = -1.0;
669  int nu3pdg = 0, nu3CCNC = -1, nu3IntType = 0;
670  double nu3VisEsl = -1.0, nu3VisEtot = -1.0;
671 
672  float nu4X = 0.0, nu4Y = 0.0, nu4Z = 0.0, nu4E = -1.0;
673  int nu4pdg = 0, nu4CCNC = -1, nu4IntType = 0;
674  double nu4VisEsl = -1.0, nu4VisEtot = -1.0;
675 
676  float nu5X = 0.0, nu5Y = 0.0, nu5Z = 0.0, nu5E = -1.0;
677  int nu5pdg = 0, nu5CCNC = -1, nu5IntType = 0;
678  double nu5VisEsl = -1.0, nu5VisEtot = -1.0;
679 
680  int pdgPartMost = 0, NumMatches = 0;
681 
682  double nu1Eff = 0.0, nu1Pur = 0.0;
683  double nu2Eff = 0.0, nu2Pur = 0.0;
684  double nu3Eff = 0.0, nu3Pur = 0.0;
685  double nu4Eff = 0.0, nu4Pur = 0.0;
686  double nu5Eff = 0.0, nu5Pur = 0.0;
687 
688  unsigned int IsMichelPrimary = 0;
689  unsigned int IsMichelPresent = 0;
690  unsigned int IsLateSlice = 0;
691 
692  double PercentNoise = -1.0;
693  double PercentNoiseByPE = -1.0;
694 
695 
696 
697  if((*slices)[i].IsNoise()) IsNoise = 1;
698  NumHitsX = (*slices)[i].NXCell();
699  NumHitsY = (*slices)[i].NYCell();
700  // slE = (*slices)[i].TotalGeV();
701 
702  // count the non-noise slices (to verify that there is only 1 noise slice)
703  if(!(*slices)[i].IsNoise()) NSlice++;
704 
705 
706 
707  // Calculate average and st. dev. for plane #, cell #, and time for
708  // hits in the slice
709 
710  double nx = 0.0, ny = 0.0;
711 
712  // loop over all hits
713  for(unsigned int j = 0; j < (*slices)[i].NCell(); ++j) {
714 
715  const rb::CellHit* h = (*slices)[i].Cell(j).get();
716 
717  // Fill hit level histos
718  if((*slices)[i].IsNoise()) {
719  fTresNoise->Fill(cal->GetTimeRes(*h));
720  fTallhitsNoise->Fill(h->TNS());
721  if(h->View() == geo::kX) fPCXallhitsNoise->Fill(h->Plane(),h->Cell());
722  if(h->View() == geo::kY) fPCYallhitsNoise->Fill(h->Plane(),h->Cell());
723  }
724  else {
725  fTres->Fill(cal->GetTimeRes(*h));
726  fTallhits->Fill(h->TNS());
727  if(h->View() == geo::kX) fPCXallhits->Fill(h->Plane(),h->Cell());
728  if(h->View() == geo::kY) fPCYallhits->Fill(h->Plane(),h->Cell());
729  }
730 
731  // Problems occur if you do the time calculations in ns since the numbers
732  // are huge when you square them (even if you use doubles.) So I will
733  // convert to microseconds to do my sums for Tave and Tsd.
734  double Tus = (h->TNS())/1000.0;
735 
736  slPE += h->PE();
737 
738  Pave += h->Plane();
739  Tave += Tus;
740  Psd += (h->Plane())*(h->Plane());
741  Tsd += Tus*Tus;
742 
743  if(h->Plane() > Phi) Phi = h->Plane();
744  if(h->Plane() < Plo) Plo = h->Plane();
745  if(h->TNS() > Thi) Thi = h->TNS();
746  if(h->TNS() < Tlo) Tlo = h->TNS();
747 
748  if(h->View() == geo::kX) {
749  CXave += h->Cell();
750  CXsd += (h->Cell())*(h->Cell());
751  nx += 1.0;
752  if(h->Cell() > CXhi) CXhi = h->Cell();
753  if(h->Cell() < CXlo) CXlo = h->Cell();
754  }
755  else if(h->View() == geo::kY) {
756  CYave += h->Cell();
757  CYsd += (h->Cell())*(h->Cell());
758  ny += 1.0;
759  if(h->Cell() > CYhi) CYhi = h->Cell();
760  if(h->Cell() < CYlo) CYlo = h->Cell();
761  }
762  else {
763  // Trigger some kind of warning that an invalid view occured.
764  CXave = 1.0e9;
765  CYave = 1.0e9;
766  }
767  } // end loop over j
768 
769  if(nx == 0.0 && ny == 0.0) {
770  // assign a negative number to attract attention
771  Pave = -100.0;
772  Tave = -100.0;
773  Psd = -100.0;
774  Tsd = -100.0;
775  }
776  else {
777  Pave = Pave/(nx+ny);
778  Tave = Tave/(nx+ny);
779  Psd = sqrt(Psd/(nx+ny) - Pave*Pave);
780  Tsd = sqrt(Tsd/(nx+ny) - Tave*Tave);
781 
782  // convert back to ns
783  Tave = Tave*1000.0;
784  Tsd = Tsd*1000.0;
785  }
786 
787  if(nx != 0.0) {
788  CXave = CXave/nx;
789  CXsd = sqrt(CXsd/nx - CXave*CXave);
790  }
791  else {
792  CXave = -100.0;
793  CXsd = -100.0;
794  }
795 
796  if(ny != 0.0) {
797  CYave = CYave/ny;
798  CYsd = sqrt(CYsd/ny - CYave*CYave);
799  }
800  else {
801  CYave = -100.0;
802  CYsd = -100.0;
803  }
804 
805 
806 
807  // additional analysis for comparison to MCTruth
808  if(fIsMC) {
809 
810  // Get the list for this slice, keep only the neutrinos that contributed,
811  // and sort the remaining list by purity.
812  std::vector<cheat::NeutrinoEffPur> nuEffPur;
813  nuEffPur = this->sortEffPur(ALLnuEffPur[i]);
814 
815  if(nuEffPur.size() > 0) {
816 
817  // if this is a noise slice, count all Eff and Pur
818  if((*slices)[i].IsNoise()) {
819  for(unsigned int a = 0; a < nuEffPur.size(); ++a) {
820  fCompNoise->Fill(nuEffPur[a].efficiency);
821  fPurNoise ->Fill(nuEffPur[a].purity);
822  }
823  }
824 
825  simb::MCNeutrino nu = nuEffPur[0].neutrinoInt->GetNeutrino();
826  simb::MCParticle Nu = nu.Nu();
827 
828  /*
829  std::cout << i << ":\tsize of nuEffPur = " << nuEffPur.size() << "\t"
830  << "Nu pdg = " << Nu.PdgCode() << "\n";
831  */
832 
833  // count neutrinos/MCTruths that contributed cell hits
834  // keep track of the nu that contribute to ALL slices
835  for(unsigned int a = 0; a < nuEffPur.size(); ++a) {
836  MCTListCellHitALL.insert(nuEffPur[a].neutrinoInt);
837  if(nuEffPur[a].neutrinoInt->NeutrinoSet())
838  NuListCellHitALL.insert(nuEffPur[a].neutrinoInt);
839  }
840 
841  // keep track of only the nu that contribute to the non-noise slices
842  if( !((*slices)[i].IsNoise()) ) {
843  for(unsigned int a = 0; a < nuEffPur.size(); ++a) {
844  MCTListCellHit.insert(nuEffPur[a].neutrinoInt);
845  if(nuEffPur[a].neutrinoInt->NeutrinoSet())
846  NuListCellHit.insert(nuEffPur[a].neutrinoInt);
847  }
848  }
849 
850 
851 
852  // BEWARE OF SEGFAULTS!!! If this MCTruth did NOT come from a neutrino,
853  // then we must skip the next section to avoid segfaluts.
854  if(nuEffPur[0].neutrinoInt->NeutrinoSet()) {
855 
856  nu1X = Nu.Vx();
857  nu1Y = Nu.Vy();
858  nu1Z = Nu.Vz();
859  nu1E = Nu.E();
860  nu1pdg = Nu.PdgCode();
861 
862  nu1CCNC = nu.CCNC();
863  nu1IntType = nu.InteractionType();
864 
865  } // end if this MCTruth = neutrino
866 
867  nu1Eff = nuEffPur[0].efficiency;
868  nu1Pur = nuEffPur[0].purity;
869 
870  nu1VisEsl = nuEffPur[0].energySlice;
871  nu1VisEtot = nuEffPur[0].energyTotal;
872 
873  if(nuEffPur.size() > 1) {
874  if(nuEffPur[1].neutrinoInt->NeutrinoSet()) {
875 
876  nu = nuEffPur[1].neutrinoInt->GetNeutrino();
877  Nu = nu.Nu();
878 
879  nu2X = Nu.Vx();
880  nu2Y = Nu.Vy();
881  nu2Z = Nu.Vz();
882  nu2E = Nu.E();
883  nu2pdg = Nu.PdgCode();
884 
885  nu2CCNC = nu.CCNC();
886  nu2IntType = nu.InteractionType();
887 
888  } // end if this MCTruth = neutrino
889 
890  nu2Eff = nuEffPur[1].efficiency;
891  nu2Pur = nuEffPur[1].purity;
892  nu2VisEsl = nuEffPur[1].energySlice;
893  nu2VisEtot = nuEffPur[1].energyTotal;
894  }
895  if(nuEffPur.size() > 2) {
896  if(nuEffPur[2].neutrinoInt->NeutrinoSet()) {
897 
898  nu = nuEffPur[2].neutrinoInt->GetNeutrino();
899  Nu = nu.Nu();
900 
901  nu3X = Nu.Vx();
902  nu3Y = Nu.Vy();
903  nu3Z = Nu.Vz();
904  nu3E = Nu.E();
905  nu3pdg = Nu.PdgCode();
906 
907  nu3CCNC = nu.CCNC();
908  nu3IntType = nu.InteractionType();
909 
910  } // end if this MCTruth = neutrino
911 
912  nu3Eff = nuEffPur[2].efficiency;
913  nu3Pur = nuEffPur[2].purity;
914  nu3VisEsl = nuEffPur[2].energySlice;
915  nu3VisEtot = nuEffPur[2].energyTotal;
916  }
917  if(nuEffPur.size() > 3) {
918  if(nuEffPur[3].neutrinoInt->NeutrinoSet()) {
919 
920  nu = nuEffPur[3].neutrinoInt->GetNeutrino();
921  Nu = nu.Nu();
922 
923  nu4X = Nu.Vx();
924  nu4Y = Nu.Vy();
925  nu4Z = Nu.Vz();
926  nu4E = Nu.E();
927  nu4pdg = Nu.PdgCode();
928 
929  nu4CCNC = nu.CCNC();
930  nu4IntType = nu.InteractionType();
931 
932  } // end if this MCTruth = neutrino
933  nu4Eff = nuEffPur[3].efficiency;
934  nu4Pur = nuEffPur[3].purity;
935  nu4VisEsl = nuEffPur[3].energySlice;
936  nu4VisEtot = nuEffPur[3].energyTotal;
937  }
938  if(nuEffPur.size() > 4) {
939  if(nuEffPur[4].neutrinoInt->NeutrinoSet()) {
940 
941  nu = nuEffPur[4].neutrinoInt->GetNeutrino();
942  Nu = nu.Nu();
943 
944  nu5X = Nu.Vx();
945  nu5Y = Nu.Vy();
946  nu5Z = Nu.Vz();
947  nu5E = Nu.E();
948  nu5pdg = Nu.PdgCode();
949 
950  nu5CCNC = nu.CCNC();
951  nu5IntType = nu.InteractionType();
952 
953  } // end if this MCTruth = neutrino
954  nu5Eff = nuEffPur[4].efficiency;
955  nu5Pur = nuEffPur[4].purity;
956  nu5VisEsl = nuEffPur[4].energySlice;
957  nu5VisEtot = nuEffPur[4].energyTotal;
958  }
959  } // end if nuEffPur.size() > 0
960 
961  NumMatches = nuEffPur.size();
962 
963  // get info about primary particle contributing to the slice
964  const std::vector< const sim::Particle * > particles = BT->HitsToParticle((*slices)[i].AllCells());
965 
966  if(particles.size() > 0) {
967  pdgPartMost = particles[0]->PdgCode();
968 
969  // decide if michel e is present or primary in this slice
970  if(pdgPartMost == 11 && particles[0]->Process() == "Decay") IsMichelPrimary = 1;
971 
972  for(unsigned int a = 0; a < particles.size(); ++a) {
973  if(particles[a]->PdgCode() == 11 && particles[a]->Process() == "Decay") IsMichelPresent = 1;
974  }
975 
976  } // end if(particles.size() > 0)
977 
978  // decide if this is a "late" slice
979  if(nuEffPur.size() > 0) {
980  if(nuEffPur[0].neutrinoInt->NParticles() > 0) {
981  const simb::MCParticle p = nuEffPur[0].neutrinoInt->GetParticle(0);
982  double T = p.T();
983  if( (T-Tave) > 200.0) IsLateSlice = 1;
984  }
985  }
986 
987  if(NumHitsX >= 3 && NumHitsY >= 3 &&
988  nu1Eff >= 0.9 && nu1Pur >= 0.9) NgoodSL++;
989 
990  // count the percentage of this slice that is true noise
991  double Nhit = 0.0;
992  double NhitPE = 0.0;
993  double NNhit = 0.0;
994  double NNhitPE = 0.0;
995 
996  for(unsigned int hi = 0; hi < (*slices)[i].NCell(); ++hi) {
997 
998  art::Ptr<rb::CellHit> hit = (*slices)[i].Cell(hi);
999 
1000  Nhit += 1.0;
1001  NhitPE += hit->PE();
1002 
1003  if(BT->IsNoise(hit)) {
1004  NNhit += 1.0;
1005  NNhitPE += hit->PE();
1006  }
1007 
1008  } // end loop over hits in this slice
1009 
1010  PercentNoise = NNhit/Nhit;
1011  PercentNoiseByPE = NNhitPE/NhitPE;
1012 
1013  } // end if(IsMC)
1014 
1015 
1016 
1017 
1018 
1019 
1020  // fill slice vectors with values
1021 
1022  IsNoiseVec .push_back(IsNoise);
1023  NumHitsXVec.push_back(NumHitsX);
1024  NumHitsYVec.push_back(NumHitsY);
1025  slEVec .push_back(slE);
1026  slPEVec .push_back(slPE);
1027 
1028  PaveVec .push_back(Pave);
1029  TaveVec .push_back(Tave);
1030  CXaveVec.push_back(CXave);
1031  CYaveVec.push_back(CYave);
1032 
1033  PsdVec .push_back(Psd);
1034  TsdVec .push_back(Tsd);
1035  CXsdVec.push_back(CXsd);
1036  CYsdVec.push_back(CYsd);
1037 
1038  PhiVec .push_back(Phi);
1039  ThiVec .push_back(Thi);
1040  CXhiVec.push_back(CXhi);
1041  CYhiVec.push_back(CYhi);
1042  PloVec .push_back(Plo);
1043  TloVec .push_back(Tlo);
1044  CXloVec.push_back(CXlo);
1045  CYloVec.push_back(CYlo);
1046 
1047  nu1XVec.push_back(nu1X);
1048  nu1YVec.push_back(nu1Y);
1049  nu1ZVec.push_back(nu1Z);
1050  nu1EVec.push_back(nu1E);
1051  nu2XVec.push_back(nu2X);
1052  nu2YVec.push_back(nu2Y);
1053  nu2ZVec.push_back(nu2Z);
1054  nu2EVec.push_back(nu2E);
1055  nu3XVec.push_back(nu3X);
1056  nu3YVec.push_back(nu3Y);
1057  nu3ZVec.push_back(nu3Z);
1058  nu3EVec.push_back(nu3E);
1059  nu4XVec.push_back(nu4X);
1060  nu4YVec.push_back(nu4Y);
1061  nu4ZVec.push_back(nu4Z);
1062  nu4EVec.push_back(nu4E);
1063  nu5XVec.push_back(nu5X);
1064  nu5YVec.push_back(nu5Y);
1065  nu5ZVec.push_back(nu5Z);
1066  nu5EVec.push_back(nu5E);
1067 
1068  nu1pdgVec .push_back(nu1pdg);
1069  nu1CCNCVec .push_back(nu1CCNC);
1070  nu1IntTypeVec .push_back(nu1IntType);
1071  nu1VisEslVec .push_back(nu1VisEsl);
1072  nu1VisEtotVec .push_back(nu1VisEtot);
1073  nu2pdgVec .push_back(nu2pdg);
1074  nu2CCNCVec .push_back(nu2CCNC);
1075  nu2IntTypeVec .push_back(nu2IntType);
1076  nu2VisEslVec .push_back(nu2VisEsl);
1077  nu2VisEtotVec .push_back(nu2VisEtot);
1078  nu3pdgVec .push_back(nu3pdg);
1079  nu3CCNCVec .push_back(nu3CCNC);
1080  nu3IntTypeVec .push_back(nu3IntType);
1081  nu3VisEslVec .push_back(nu3VisEsl);
1082  nu3VisEtotVec .push_back(nu3VisEtot);
1083  nu4pdgVec .push_back(nu4pdg);
1084  nu4CCNCVec .push_back(nu4CCNC);
1085  nu4IntTypeVec .push_back(nu4IntType);
1086  nu4VisEslVec .push_back(nu4VisEsl);
1087  nu4VisEtotVec .push_back(nu4VisEtot);
1088  nu5pdgVec .push_back(nu5pdg);
1089  nu5CCNCVec .push_back(nu5CCNC);
1090  nu5IntTypeVec .push_back(nu5IntType);
1091  nu5VisEslVec .push_back(nu5VisEsl);
1092  nu5VisEtotVec .push_back(nu5VisEtot);
1093 
1094  pdgPartMostVec.push_back(pdgPartMost);
1095  NumMatchesVec .push_back(NumMatches);
1096 
1097  nu1EffVec.push_back(nu1Eff);
1098  nu1PurVec.push_back(nu1Pur);
1099  nu2EffVec.push_back(nu2Eff);
1100  nu2PurVec.push_back(nu2Pur);
1101  nu3EffVec.push_back(nu3Eff);
1102  nu3PurVec.push_back(nu3Pur);
1103  nu4EffVec.push_back(nu4Eff);
1104  nu4PurVec.push_back(nu4Pur);
1105  nu5EffVec.push_back(nu5Eff);
1106  nu5PurVec.push_back(nu5Pur);
1107 
1108  IsMichelPrimaryVec.push_back(IsMichelPrimary);
1109  IsMichelPresentVec.push_back(IsMichelPresent);
1110  IsLateSliceVec .push_back(IsLateSlice);
1111 
1112  PercentNoiseVec .push_back(PercentNoise);
1113  PercentNoiseByPEVec.push_back(PercentNoiseByPE);
1114 
1115  // fill histos
1116  if((*slices)[i].IsNoise()) {
1117 
1118  if(NumHitsX+NumHitsY == 0.0)
1119  std::cout << "\n\nEmpty Noise Slice! " << evt.id() << "\n\n\n";
1120  fNumHitsNoise->Fill(NumHitsX+NumHitsY);
1121  fPCXaveNoise ->Fill(Pave,CXave);
1122  fPCYaveNoise ->Fill(Pave,CYave);
1123  fTaveNoise ->Fill(Tave);
1124 
1125  NhitNoise += NumHitsX+NumHitsY;
1126  NhitTotal += NumHitsX+NumHitsY;
1127 
1128  }
1129  else {
1130 
1131  fNumHits->Fill(NumHitsX+NumHitsY);
1132  fComp ->Fill(nu1Eff);
1133  fPur ->Fill(nu1Pur);
1134  fCompPur->Fill(nu1Pur,nu1Eff);
1135  fPCXave ->Fill(Pave,CXave);
1136  fPCYave ->Fill(Pave,CYave);
1137  fPCXdiff->Fill(Phi-Plo,CXhi-CXlo);
1138  fPCYdiff->Fill(Phi-Plo,CYhi-CYlo);
1139  fTave ->Fill(Tave);
1140  fTsd ->Fill(Tsd);
1141  fSlDur ->Fill(Thi-Tlo);
1142 
1143  NhitTotal += NumHitsX+NumHitsY;
1144 
1145  /*
1146  if(Thi-Tlo >= 1000.0 && Thi-Tlo <= 1070.0)
1147  std::cout << "\nevent = " << evt.id() << " slice = " << i << " time = " << Thi-Tlo;
1148  */
1149 
1150  if(fIsMC) {
1151  fNSvsVE->Fill(nu1VisEtot);
1152  if(nu1Eff >= 0.9 && nu1Pur >= 0.9) fNGvsVE->Fill(nu1VisEtot);
1153  }
1154 
1155  }
1156 
1157  } // end loop over slices (i)
1158 
1159 
1160 
1161  // accumulate the set of MCTruths that contributed FLSHits
1162  if(fIsMC) {
1163 
1164  // Get the particle list vector
1166  evt.getByLabel("geantgen",parts);
1167 
1168  // Retrive the associated MCTruth(s)
1169  art::FindManyP<simb::MCTruth> truthAssn(parts,evt,"geantgen");
1170 
1171  // Loop through MCTruth object to accumulate list of keepers
1172  for(size_t p = 0; p < parts->size(); ++p) {
1173  for(size_t t = 0; t < truthAssn.at(p).size(); ++t) {
1174  art::Ptr<simb::MCTruth> mcTru = truthAssn.at(p)[t];
1175  MCTListFLSHitALL.insert(mcTru);
1176  if(mcTru->NeutrinoSet())
1177  NuListFLSHitALL.insert(mcTru);
1178  }
1179  }
1180 
1181  }
1182 
1183 
1184 
1185  /// TODO: Add vectors of info about each neutrino?
1186 
1187 
1188 
1189  fNumNu->Fill(NuListFLSHitALL.size());
1190  fNumMCT->Fill(MCTListFLSHitALL.size());
1191 
1192  NumSlice = NSlice;
1193  fNumSlices->Fill(NSlice);
1194 
1195  fNumMCTSlices->Fill(MCTListCellHit.size(),NSlice);
1196 
1197  NumNoise = NNoiseSlice;;
1198  fNumNoise->Fill(NumNoise);
1199 
1200  NumNuFLSHitALL = NuListFLSHitALL.size();
1201  NumNuCellHitALL = NuListCellHitALL.size();
1202  NumNuCellHit = NuListCellHit.size();
1203 
1204  NumMCTFLSHitALL = MCTListFLSHitALL.size();
1205  NumMCTCellHitALL = MCTListCellHitALL.size();
1206  NumMCTCellHit = MCTListCellHit.size();
1207 
1208  fEvtNoiseFraction->Fill((double)NhitNoise/(double)NhitTotal);
1209 
1210  fOutTree->Fill();
1211 
1212  /*
1213  std::cout << "\n\nFSL = " << NuListFLSHitALL.size() << "\t"
1214  << "cellALL = " << NuListCellHitALL.size() << "\t"
1215  << "cell = " << NuListCellHit.size() << "\n";
1216  std::cout << "FSL = " << MCTListFLSHitALL.size() << "\t"
1217  << "cellALL = " << MCTListCellHitALL.size() << "\t"
1218  << "cell = " << MCTListCellHit.size() << "\n\n";
1219  */
1220 
1221  return;
1222  }
1223 
1224  //----------------------------------------------------------------------
1225  std::vector<cheat::NeutrinoEffPur> SlicerAna::sortEffPur(const std::vector<cheat::NeutrinoEffPur>& nuEP)
1226  {
1227  std::vector<cheat::NeutrinoEffPur> results;
1228 
1229  // only push in stuff with non-zero purity
1230  for(unsigned int i = 0; i < nuEP.size(); ++i) {
1231  if(nuEP[i].purity > 0.0) results.push_back(nuEP[i]);
1232  }
1233 
1234  // sort by purity (bubble sort)
1235 
1236  if(results.size() > 0) {
1237  bool done = false;
1238  while(!done) {
1239  done = true;
1240  for(unsigned int i = 1; i < results.size(); ++i) {
1241  if(results[i-1].purity < results[i].purity) {
1242  cheat::NeutrinoEffPur temp = results[i];
1243  results[i] = results[i-1];
1244  results[i-1] = temp;
1245  done = false;
1246  } // end if
1247  } // end for i
1248  } // end while !done
1249  } // end if(results.size() > 0)
1250 
1251  return results;
1252 
1253  }
1254 
1255  //----------------------------------------------------------------------
1257  {
1258  IsNoiseVec.clear();
1259  NumHitsXVec.clear();
1260  NumHitsYVec.clear();
1261  slEVec.clear();
1262  slPEVec.clear();
1263 
1264  PaveVec.clear();
1265  TaveVec.clear();
1266  CXaveVec.clear();
1267  CYaveVec.clear();
1268 
1269  PsdVec.clear();
1270  TsdVec.clear();
1271  CXsdVec.clear();
1272  CYsdVec.clear();
1273 
1274  PhiVec.clear();
1275  PloVec.clear();
1276  ThiVec.clear();
1277  TloVec.clear();
1278  CXhiVec.clear();
1279  CXloVec.clear();
1280  CYhiVec.clear();
1281  CYloVec.clear();
1282 
1283  nu1XVec.clear();
1284  nu1YVec.clear();
1285  nu1ZVec.clear();
1286  nu1EVec.clear();
1287  nu2XVec.clear();
1288  nu2YVec.clear();
1289  nu2ZVec.clear();
1290  nu2EVec.clear();
1291  nu3XVec.clear();
1292  nu3YVec.clear();
1293  nu3ZVec.clear();
1294  nu3EVec.clear();
1295  nu4XVec.clear();
1296  nu4YVec.clear();
1297  nu4ZVec.clear();
1298  nu4EVec.clear();
1299  nu5XVec.clear();
1300  nu5YVec.clear();
1301  nu5ZVec.clear();
1302  nu5EVec.clear();
1303 
1304  nu1pdgVec.clear();
1305  nu1CCNCVec.clear();
1306  nu1IntTypeVec.clear();
1307  nu1VisEslVec.clear();
1308  nu1VisEtotVec.clear();
1309  nu2pdgVec.clear();
1310  nu2CCNCVec.clear();
1311  nu2IntTypeVec.clear();
1312  nu2VisEslVec.clear();
1313  nu2VisEtotVec.clear();
1314  nu3pdgVec.clear();
1315  nu3CCNCVec.clear();
1316  nu3IntTypeVec.clear();
1317  nu3VisEslVec.clear();
1318  nu3VisEtotVec.clear();
1319  nu4pdgVec.clear();
1320  nu4CCNCVec.clear();
1321  nu4IntTypeVec.clear();
1322  nu4VisEslVec.clear();
1323  nu4VisEtotVec.clear();
1324  nu5pdgVec.clear();
1325  nu5CCNCVec.clear();
1326  nu5IntTypeVec.clear();
1327  nu5VisEslVec.clear();
1328  nu5VisEtotVec.clear();
1329 
1330  pdgPartMostVec.clear();
1331  NumMatchesVec.clear();
1332 
1333  nu1EffVec.clear();
1334  nu1PurVec.clear();
1335  nu2EffVec.clear();
1336  nu2PurVec.clear();
1337  nu3EffVec.clear();
1338  nu3PurVec.clear();
1339  nu4EffVec.clear();
1340  nu4PurVec.clear();
1341  nu5EffVec.clear();
1342  nu5PurVec.clear();
1343 
1344  IsMichelPrimaryVec.clear();
1345  IsMichelPresentVec.clear();
1346  IsLateSliceVec.clear();
1347 
1348  PercentNoiseVec.clear();
1349  PercentNoiseByPEVec.clear();
1350 
1351  }
1352 
1353 } // end namespace slicer
1354 ////////////////////////////////////////////////////////////////////////
1355 
1356 
1357 
1359 
1360 namespace slicer
1361 {
1363 }
double E(const int i=0) const
Definition: MCParticle.h:232
std::vector< int > nu2pdgVec
PDG code for second nu associated with slice.
float TNS() const
Definition: CellHit.h:46
std::vector< int > nu5pdgVec
PDG code for fifth nu associated with slice.
std::vector< double > nu4VisEtotVec
total visible energy for the fourth best matched nu in the entire event
std::vector< double > slEVec
slice.TotalGeV()
SubRunNumber_t subRun() const
Definition: Event.h:72
std::string fCAFLabel
label for the process that made the standard records
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
std::vector< double > nu4VisEslVec
visible energy for the fourth best matched nu in the slice
double GetTimeRes(rb::CellHit const &cellhit)
std::vector< double > nu2VisEtotVec
total visible energy for the second best matched nu in the entire event
std::vector< unsigned int > IsMichelPresentVec
did a Michel e contribute to this slice?
int NumMCTCellHit
total number of MCTruths in non-noise slices counted by CellHits
std::vector< float > nu2YVec
y location for second nu associated with slice
std::vector< double > nu5VisEtotVec
total visible energy for the fifth best matched nu in the entire event
std::vector< int > pdgPartMostVec
the PDG code for the particle that contributed the most energy to this slice
std::vector< double > nu2VisEslVec
visible energy for the second best matched nu in the slice
void analyze(const art::Event &evt)
unsigned short Plane() const
Definition: CellHit.h:39
std::vector< double > nu3EffVec
efficiency for third nu associated with slice
std::vector< double > PaveVec
average plane number for hits in a non-noise slice
geo::View_t View() const
Definition: CellHit.h:41
std::vector< int > NumHitsYVec
number of Y hits per slice
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH1 * ratio(TH1 *h1, TH1 *h2)
std::vector< int > nu5IntTypeVec
Interaction type for fifth neutrino associated with slice.
int run
run number
std::vector< int > NumHitsXVec
number of X hits per slice
Vertical planes which measure X.
Definition: PlaneGeo.h:28
std::vector< float > nu4EVec
Energy for fourth nu associated with slice.
std::vector< double > CYhiVec
highest Y cell in the slice
std::vector< float > nu5YVec
y location for fifth nu associated with slice
std::vector< double > PercentNoiseVec
Percentage of slice that is true noise (by hit)
std::vector< unsigned int > IsMichelPrimaryVec
is this slice a Michel e? (the primary contribution by energy was from a michel e) ...
std::vector< double > PloVec
lowest plane in the slice
std::vector< int > nu4IntTypeVec
Interaction type for fourth neutrino associated with slice.
std::vector< double > TloVec
lowest plane in the slice
int NumMCTFLSHitALL
total number of MCTruths counted by FLSHits
recovalid::CAFCutter fCAFCutter
the CAFCutter helper class
DEFINE_ART_MODULE(TestTMapFile)
std::vector< float > nu3EVec
Energy for third nu associated with slice.
virtual void endSubRun(const art::SubRun &sr)
Particle class.
std::vector< float > nu1YVec
y location for primary nu associated with slice
std::vector< float > nu2ZVec
z location for second nu associated with slice
std::vector< std::vector< cheat::NeutrinoEffPur > > SlicesToMCTruthsTable(const std::vector< const rb::Cluster * > &sliceList) const
Given ALL the slices in an event, including the noise slice, returns a vector of vector of structures...
std::vector< double > nu3VisEtotVec
total visible energy for the third best matched nu in the entire event
std::vector< double > nu1EffVec
efficiency for primary nu associated with slice
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
std::vector< int > nu3pdgVec
PDG code for third nu associated with slice.
bool IsNoise(const art::Ptr< rb::CellHit > &hit) const
Is this hit not associated with any particles?
std::vector< double > nu5VisEslVec
visible energy for the fifth best matched nu in the slice
std::vector< double > PsdVec
st. dev. of plane number for hits in a non-noise slice
std::vector< float > nu2XVec
x location for second nu associated with slice
int NumMCTCellHitALL
total number of MCTruths counted by CellHits
std::vector< double > nu3VisEslVec
visible energy for the third best matched nu in the slice
std::vector< float > nu5XVec
x location for fifth nu associated with slice
std::vector< int > nu1CCNCVec
CCNC for primary nu associated with slice.
std::vector< double > nu1VisEslVec
visible energy for the best matched nu in the slice
std::vector< float > nu5ZVec
z location for fifth nu associated with slice
unsigned short Cell() const
Definition: CellHit.h:40
std::vector< double > CXsdVec
st. dev. cell number for X hits in a non-noise slice
int InteractionType() const
Definition: MCNeutrino.h:150
bool passCuts(int cut, const caf::StandardRecord *sr)
Definition: CAFCutter.cxx:37
TSpline3 hi("hi", xhi, yhi, 18,"0")
std::vector< double > nu1VisEtotVec
total visible energy for the best matched nu in the entire event
int NumNoise
number of noise slices
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
int NumNuFLSHitALL
number of nu counted by FLSHits
std::vector< double > nu5PurVec
purity for fifth nu associated with slice
std::vector< double > TsdVec
st. dev. time for hits in a non-noise slice
int fCutType
what cuts should we apply (see CAFCutter.h)
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
int NumNuCellHit
number of nu in non-noise slices counted by CellHits
int NumSlice
number of non-noise slices
std::vector< double > PercentNoiseByPEVec
Percentage of slice that is true noise (weighted by PE)
std::vector< float > nu1EVec
Energy for primary nu associated with slice.
std::vector< double > nu2PurVec
purity for second nu associated with slice
const double a
std::vector< int > NumMatchesVec
number of MCTruth matches per slice
std::vector< float > nu5EVec
Energy for fifth nu associated with slice.
std::vector< float > nu2EVec
Energy for second nu associated with slice.
std::vector< double > CYaveVec
average cell number for Y hits in a non-noise slice
T get(std::string const &key) const
Definition: ParameterSet.h:231
bool fApplyCAFCuts
should we apply CAF cuts?
int evt
double T(const int i=0) const
Definition: MCParticle.h:223
#define pot
std::vector< double > nu3PurVec
purity for third nu associated with slice
float PE() const
Definition: CellHit.h:42
std::vector< unsigned int > IsLateSliceVec
is the average time for this slice > 200 ns after the time of the neutrino interaction?
std::string fSlicerLabel
label for the process that made the slices
SlicerAna(fhicl::ParameterSet const &pset)
caf::StandardRecord * sr
const double j
Definition: BetheBloch.cxx:29
fvar< T > Phi(const fvar< T > &x)
Definition: Phi.hpp:12
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
std::vector< double > nu4PurVec
purity for fourth nu associated with slice
std::vector< double > nu2EffVec
efficiency for second nu associated with slice
std::vector< float > slPEVec
sum of PE for all hits in the slice
Definition: run.py:1
std::vector< int > nu5CCNCVec
CCNC for fifth nu associated with slice.
OStream cout
Definition: OStream.cxx:6
std::vector< double > CXaveVec
average cell number for X hits in a non-noise slice
bool fIsMC
is the file being processed MC?
std::vector< double > CXloVec
lowest X cell in the slice
EventNumber_t event() const
Definition: EventID.h:116
void reconfigure(const fhicl::ParameterSet &pset)
std::vector< int > nu3CCNCVec
CCNC for third nu associated with slice.
std::vector< int > nu1IntTypeVec
Interaction type for primary neutrino associated with slice.
std::vector< double > TaveVec
average time for hits in a non-noise slice
double Vx(const int i=0) const
Definition: MCParticle.h:220
std::vector< double > CYsdVec
st. dev. cell number for Y hits in a non-noise slice
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
T * make(ARGS...args) const
std::vector< int > nu2IntTypeVec
Interaction type for second neutrino associated with slice.
std::vector< int > nu3IntTypeVec
Interaction type for third neutrino associated with slice.
std::vector< unsigned int > IsNoiseVec
0 = non-noise slice, 1 = noise slice
int subrun
subrun number
std::vector< float > nu3YVec
y location for third nu associated with slice
int NumNuCellHitALL
number of nu counted by CellHits
std::vector< int > nu1pdgVec
PDG code for primary nu associated with slice.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::vector< float > nu3ZVec
z location for third nu associated with slice
Definition: event.h:1
std::vector< float > nu1ZVec
z location for primary nu associated with slice
std::vector< float > nu1XVec
x location for primary nu associated with slice
std::vector< cheat::NeutrinoEffPur > sortEffPur(const std::vector< cheat::NeutrinoEffPur > &nuEP)
std::vector< float > nu4YVec
y location for fourth nu associated with slice
double Vz(const int i=0) const
Definition: MCParticle.h:222
std::vector< float > nu4ZVec
z location for fourth nu associated with slice
std::vector< double > nu4EffVec
efficiency for fourth nu associated with slice
std::vector< double > CYloVec
lowest Y cell in the slice
std::vector< double > ThiVec
highest time in the slice
Class to help Slicer4D manage the collection of hits.
std::vector< float > nu4XVec
x location for fourth nu associated with slice
std::vector< double > PhiVec
highest plane in the slice
Helper class for Reco Validation modules.
Definition: CAFCutter.h:58
std::vector< double > nu1PurVec
purity for primary nu associated with slice
std::vector< double > CXhiVec
highest X cell in the slice
std::vector< double > nu5EffVec
efficiency for fifth nu associated with slice
double T
Definition: Xdiff_gwt.C:5
bool NeutrinoSet() const
Definition: MCTruth.h:77
double totgoodpot
normalized by 10^12 POT
Definition: POTSum.h:28
RunNumber_t run() const
Definition: Event.h:77
std::vector< float > nu3XVec
x location for third nu associated with slice
std::vector< int > nu4CCNCVec
CCNC for fourth nu associated with slice.
Event generator information.
Definition: MCNeutrino.h:18
void efficiency()
Definition: efficiency.C:58
EventID id() const
Definition: Event.h:56
std::vector< int > nu4pdgVec
PDG code for fourth nu associated with slice.
double Vy(const int i=0) const
Definition: MCParticle.h:221
int event
event number
float fTotalPOT
what is the POT in this file?
bool failedToGet() const
Definition: Handle.h:196
std::vector< int > nu2CCNCVec
CCNC for second nu associated with slice.
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.