SkimmerAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: SkimmerAna
3 // Module Type: analyzer
4 // File: SkimmerAna_module.cc
5 //
6 // Generated at Mon Aug 3 10:52:56 2015 by Brian Rebel using artmod
7 // from cetpkgsupport v1_08_06.
8 ////////////////////////////////////////////////////////////////////////
9 
19 #include "art_root_io/TFileService.h"
20 #include "art_root_io/TFileDirectory.h"
21 #include "fhiclcpp/ParameterSet.h"
25 
26 #include <memory>
27 #include <set>
28 #include <limits>
29 
30 #include "TH1.h"
31 #include "TH2.h"
32 #include "TTree.h"
33 
34 #include "RecoBase/Cluster.h"
35 #include "RecoBase/Prong.h"
36 #include "RecoBase/Shower.h"
37 #include "RecoBase/Vertex.h"
38 #include "CosRej/NueCosRej.h"
39 #include "CosRej/CosRejObj.h"
40 #include "ShowerLID/ShowerLID.h"
41 #include "ShowerLID/EventLID.h"
42 #include "LEM/PIDDetails.h"
43 #include "CVN/func/Result.h"
44 #include "Preselection/Veto.h"
46 #include "SummaryData/SpillData.h"
47 #include "SummaryData/POTSum.h"
49 #include "NovaDAQConventions/DAQConventions.h"
55 
56 namespace skim {
57  class SkimmerAna;
58 }
59 
61 public:
62  explicit SkimmerAna(fhicl::ParameterSet const & p);
63  virtual ~SkimmerAna();
64 
65  // Plugins should not be copied or assigned.
66  SkimmerAna(SkimmerAna const &) = delete;
67  SkimmerAna(SkimmerAna &&) = delete;
68  SkimmerAna & operator = (SkimmerAna const &) = delete;
69  SkimmerAna & operator = (SkimmerAna &&) = delete;
70 
71  // Required functions.
72  void analyze(art::Event const& e) override;
73 
74  // Selected optional functions.
75  void reconfigure(fhicl::ParameterSet const & p);
76 
77  void beginJob() override;
78  void beginRun(art::Run const& r) override;
79  void beginSubRun(art::SubRun const& sr) override;
80 
81 private:
82 
83  void FillNueHistograms (art::Event const& e);
84  void FillNumuHistograms(art::Event const& e);
85 
86  std::string fSkimmedLabel; ///< label of module creating skimmed file
87  std::string fCosmicExposureLabel; ///< label of module creating cosmic exposure
88  std::string fPOTSumLabel; ///< label of module creating pot sum
89  std::string fNumuInstance; ///< label for numu products
90  std::string fNueInstance; ///< label for nue products
91  bool fCompareNumuCuts; ///< bool for comparing the numu cuts
92  bool fCompareNueCuts; ///< bool for comparing the nue cuts
93  novadaq::cnv::DetId fDetId; ///< the id of the detector we are looking at
95 
96  TTree *fNumuTree; ///< Tree to hold values of cut parameters for selected events
97  TTree *fNueTree; ///< Tree to hold values of cut parameters for selected events
98 
99  NumuCutParameters fNumuCutParams; ///< from the ParametersNumu object
100  NueCutParameters fNueCutParams; ///< from the ParametersNumu object
101  NumuCutParameters fNumuSRParams; ///< from the ParametersNumu object
102  NueCutParameters fNueSRParams; ///< from the ParametersNumu object
103  EventInfo fEventInfo; ///< keep track of run and subrun, and event
104 
105  // POT and spill information
106  TH1F *fPOTByRun;
109 
110  // numu sanity check histograms, one for each cut
111  TH1F *fNeutrinoE;
112  TH1F *fReMId;
113  TH1F *fSliceHits;
125  TH1F *fFDContPID;
127 
135 
136  // nue sanity check histograms
150  TH1F* fFDPt;
151 
153  TH1F* fLEM;
154  TH1F* fCVN;
155  TH1F* fLID;
156  TH1F* fNDVtxX;
157  TH1F* fNDVtxY;
158  TH1F* fNDVtxZ;
159  TH1F* fNDNShowers;
160  TH1F* fNDShwVtxX;
161  TH1F* fNDShwVtxY;
162  TH1F* fNDShwVtxZ;
163  TH1F* fNDShwEndX;
164  TH1F* fNDShwEndY;
165  TH1F* fNDShwEndZ;
176 
177  // general sanity check histograms
178  //TH2F*
179 
180 };
181 
182 //------------------------------------------------------------------------------
184 : EDAnalyzer(p)
185 {
186  this->reconfigure(p);
187 
188  return;
189 }
190 
191 //------------------------------------------------------------------------------
193 {
194 }
195 
196 //------------------------------------------------------------------------------
198 {
199 
200  // make histograms for each selection criteria
202 
203  std::string infost("run/I:subrun/I:event/I");
204 
205  fPOTByRun = tfs->make<TH1F>("potByRun", ";Run Number;POT (#times 10^{15})", 10000, 10000., 20000.);
206  fGoodPOTByRun = tfs->make<TH1F>("goodPOTByRun", ";Run Number;POT (#times 10^{15})", 10000, 10000., 20000.);
207  fCosmicExposureByRun = tfs->make<TH1F>("cosmicExposureByRun", ";Run Number;Exposure (s)", 10000, 10000., 20000.);
208 
209  if( !fNumuInstance.empty() ){
210  fNeutrinoE = tfs->make<TH1F>("neutrinoE", ";E_{#nu} (GeV);Events", 100, 0., 10. );
211  fReMId = tfs->make<TH1F>("remid", ";ReMId;Events", 200, -1., 1. );
212  fSliceHits = tfs->make<TH1F>("sliceHits", ";Hits In Slice;Events", 100, 0., 400. );
213  fSliceContigPlanes = tfs->make<TH1F>("sliceContigPlanes", ";Contiguous Planes in Slice;Events", 50, 0., 50. );
214  fCellsFromEdge = tfs->make<TH1F>("cellsFromEdge", ";Cells From Edge;Events", 50, 0., 50. );
215  fQePId1Track = tfs->make<TH1F>("qepid1Track", ";QePId_{1 track};Events", 200, -1., 1. );
216  fQePId2Track = tfs->make<TH1F>("qepid2Track", ";QePId_{2 tracks};Events", 200, -1., 1. );
217  fFDKalmanFwdCell = tfs->make<TH1F>("fdKalmanFwdCell", ";Kalman Foward Cell;Events", 100, 0., 100. );
218  fFDKalmanBakCell = tfs->make<TH1F>("fdKalmanBakCell", ";Kalman Backward Cell;Events", 100, 0., 100. );
219  fFDCosFwdCell = tfs->make<TH1F>("fdCosFwdCell", ";cos Foward Cell;Events", 100, 0., 100. );
220  fFDCosBakCell = tfs->make<TH1F>("fdCosBakCell", ";cos Backward Cell;Events", 100, 0., 100. );
221  fFDPlanesToFront = tfs->make<TH1F>("fdPlanesToFront", ";Planes to Front;Events", 100, 0., 100. );
222  fFDPlanesToBack = tfs->make<TH1F>("fdPlanesToBack", ";Planes to Back;Events", 100, 0., 100. );
223  fFDKalmanAngle = tfs->make<TH1F>("fdKalmanAngle", ";Kalman Angle;Events", 200, -4., 4. );
224  fFDContPID = tfs->make<TH1F>("fdContPID", ";ContPID;Events", 200, -1., 1. );
225  fFDSliceHits = tfs->make<TH1F>("fdSliceHits", ";Slice Hits;Events", 500, 0., 500. );
226 
227  fNDMaxPlaneHit = tfs->make<TH1F>("ndMaxPlaneHit", ";Max Hit Plane;Events", 220, 0., 220. );
228  fNDMinPlaneHit = tfs->make<TH1F>("ndMinPlaneHit", ";Min Hit Plane;Events", 220, 0., 220. );
229  fNDTrackStartZ = tfs->make<TH1F>("ndTrackStartZ", ";Track Starting Z Pos (cm);Events", 200, 0., 2000.);
230  fNDTrackStopZ = tfs->make<TH1F>("ndTrackStopZ", ";Track Ending Z Pos (cm);Events", 200, 0., 2000.);
231  fNDKalmanFwdCell = tfs->make<TH1F>("ndKalmanFwdCell", ";Kalman Foward Cell;Events", 100, 0., 100. );
232  fNDKalmanBakCell = tfs->make<TH1F>("ndKalmanBakCell", ";Kalman Backward Cell;Events", 100, 0., 100. );
233  fNDHadCalEnergySum = tfs->make<TH1F>("ndHadCalEnergySum", ";Hadronic Energy (GeV);Events", 100, 0., 10. );
234 
235  std::string paramst("fBestPIDTrack/F:fNeutrinoE/F:fNDHadronicCal/F:");
236  paramst.append("fReMIdVal/D:fQePIDVal/F:fQePIDNTracks/F:");
237  paramst.append("fSliceHits/F:fSliceContigPlanes/F:fSliceMaxPlane/F:fSliceMinPlane/F:fSliceCellsFromEdge/F:");
238  paramst.append("fPlanesToFront/F:fPlanesToBack/F:fCosRejKalFwdCell/F:fCosRejKalBakCell/F:fCosRejAngleKal/F:");
239  paramst.append("fCosRejCosBakCell/F:fCosRejCosFwdCell/F:fCosRejConCosPID/F:fCosRejKalBakCellND/F:fCosRejKalFwdCellND/F:");
240  paramst.append("fCosRejKalYPosAtTrans/F:fTrackStartZ/F:fTrackStopZ/F:fTrackCCE/F");
241  paramst.append("fMinWestDist/F:fMinEastDist/F:fMinBotDist/F:fMinTopDist/F:fMinFrontDist/F:fMinBackDist:"); /// 2017
242  fNumuTree = tfs->make<TTree>("Numu", "Skimmed Numu Interaction Cut Parameters");
243  fNumuTree->Branch("info", &fEventInfo, infost .c_str());
244  fNumuTree->Branch("params", &fNumuCutParams, paramst.c_str());
245  fNumuTree->Branch("sr", &fNumuSRParams, paramst.c_str());
246  }
247 
248  if( !fNueInstance.empty() ){
249  fLEM = tfs->make<TH1F>("lem", ";LEM;Slices", 100, 0., 1. );
250  fCVN = tfs->make<TH1F>("cvn", ";CVN;Slices", 100, 0., 1. );
251  fLID = tfs->make<TH1F>("lid", ";LID;Slices", 240, -0.1, 1.1);
252  fCalorimetricEnergy = tfs->make<TH1F>("calorimetricE", ";Calorimetric E (GeV)", 100, 0., 50.);
253 
254  fFDSliceHitsPerPlane = tfs->make<TH1F>("fdSliceHitsPerPlane", ";Slice Hits/Plane;Slices", 100, 0., 50.);
255  fFDSliceHitsXView = tfs->make<TH1F>("fdSliceHitsXView", ";Slice Hits (X view);Slices", 100, 0., 100.);
256  fFDSliceHitsYView = tfs->make<TH1F>("fdSliceHitsYView", ";Slice Hits (Y view);Slices", 100, 0., 100.);
257  fFDSliceHitsViewAsym = tfs->make<TH1F>("fdSliceHitsViewAsym", ";Slice Hits View Asymmetry;Slices", 200, -10., 10.);
258  fFDSliceCosAngleShw = tfs->make<TH1F>("fdSliceCosAngleShw", ";cos(#theta_{showers});Slices", 100, -1., 1.);
259  fFDSliceDistShwVtx = tfs->make<TH1F>("fdSliceDistShwVtx", ";#DeltaL_{shower,vertex};Slices", 100, 0., 400.);
260  fFDSliceFracShwHits = tfs->make<TH1F>("fdSliceFracShwHits", ";Hits_{shower}/Hits_{slice};Slices", 50, 0., 2.);
261  fFDVtxDistToEastWall = tfs->make<TH1F>("fdVtxDistToEastWall", ";Distance to East Wall;Showers", 100, 0., 100.);
262  fFDVtxDistToWestWall = tfs->make<TH1F>("fdVtxDistToWestWall", ";Distance to West Wall;Showers", 100, 0., 100.);
263  fFDVtxDistToTop = tfs->make<TH1F>("fdVtxDistToTop", ";Distance to Top;Showers", 100, 0., 400.);
264  fFDVtxDistToBottom = tfs->make<TH1F>("fdVtxDistToBottom", ";Distance to Bottom;Showers", 100, 0., 400.);
265  fFDVtxDistToFront = tfs->make<TH1F>("fdVtxDistToFront", ";Distance to Front;Showers", 100, 0., 100.);
266  fFDVtxDistToBack = tfs->make<TH1F>("fdVtxDistToBack", ";Distance to Back;Showers", 100, 0., 400.);
267  fFDPt = tfs->make<TH1F>("fdPt", ";p_{t};Showers", 50, 0., 2.);
268 
269  // now the ND histograms
270  fNDNShowers = tfs->make<TH1F>("ndNshowers", ";Nnumber of Showers;Slices", 10, 0., 10.);
271  fNDShowerCalE = tfs->make<TH1F>("ndshowercale", ";Shower E_{calorimetric};Slices", 120, 0., 6.);
272  fNDShowerCalEfraction = tfs->make<TH1F>("ndshowercalefrac", ";Fraction of Slice E in Shower;Slices", 110, 0., 1.1);
273  fNDSliceNoShowerCalE = tfs->make<TH1F>("ndnoshowercale", ";(Slice -Shower) E_{calorimetric};Slices",120, 0., 6.);
274  fNDShowerDistToVtx = tfs->make<TH1F>("ndNshowergap", ";Shower Distance to Vertex;Slices", 25, 0., 200.);
275  fNDVtxX = tfs->make<TH1F>("ndVtxX", ";x_{vertex} (cm);Showers", 200, -200., 200.);
276  fNDVtxY = tfs->make<TH1F>("ndVtxY", ";y_{vertex} (cm);Showers", 200, -200., 200.);
277  fNDVtxZ = tfs->make<TH1F>("ndVtxZ", ";z_{vertex} (cm);Showers", 700, 0., 1400.);
278  fNDShwVtxX = tfs->make<TH1F>("ndShwVtxX", ";x_{vertex} (cm);Showers", 200, -200., 200.);
279  fNDShwVtxY = tfs->make<TH1F>("ndShwVtxY", ";y_{vertex} (cm);Showers", 200, -200., 200.);
280  fNDShwVtxZ = tfs->make<TH1F>("ndShwVtxZ", ";z_{vertex} (cm);Showers", 700, 0., 1400.);
281  fNDShwEndX = tfs->make<TH1F>("ndShwEndX", ";x_{vertex} (cm);Showers", 200, -200., 200.);
282  fNDShwEndY = tfs->make<TH1F>("ndShwEndY", ";y_{vertex} (cm);Showers", 200, -200., 200.);
283  fNDShwEndZ = tfs->make<TH1F>("ndShwEndZ", ";z_{vertex} (cm);Showers", 700, 0., 1400.);
284  fNDSliceHits = tfs->make<TH1F>("ndSliceHits", ";Hits in Slice;Slices", 500, 0., 500.);
285  fNDCalorimetricE = tfs->make<TH1F>("ndCalorimetricE", ";E_{calorimetric} (GeV);Showers", 120, 0., 6.);
286  fNDCalEperHit = tfs->make<TH1F>("ndCalEperHit", ";E_{calorimetric}/ hit ;Slices", 100, 0., 0.1);
287  fNDShowerLength = tfs->make<TH1F>("ndShowerLength", ";L_{shower} (cm);Showers", 700, 0., 1400.);
288  fNDSliceNHitsPerPlane = tfs->make<TH1F>("ndSliceHitsPerPlane", ";Slice Hits/Plane;Slices", 100, 0., 50.);
289  fNDPlanesToFront = tfs->make<TH1F>("ndPlanesToFront", ";Planes to Front;Events", 600, 0., 600.);
290 
291  std::string paramst("fNumPlane/F:fPlanesToFront/F:fCellsPerPlane/F:fNumXCell/F:fNumYCell/F:");
292  paramst.append("fSliceHits/F:fHitAsymmetry/F:fNumShowers/F:fShowerNumXCell/F:fShowerNumYCell/F:");
293  paramst.append("fCosShowers/F:fShowerDirX/F:fShowerDirY/F:fShowerDirZ/F:");
294  paramst.append("fShowerCalE/F:fShowerLength/F:fFracShowerHits/F:");
295  paramst.append("fShowerVertexX/F:fShowerVertexY/F:fShowerVertexZ/F:");
296  paramst.append("fShowerEndX/F:fShowerEndY/F:fShowerEndZ/F:");
297  paramst.append("fShowerMinX/F:fShowerMinY/F:fShowerMaxX/F:fShowerMaxY/F:fShowerMinVtxZ/F:");
298  paramst.append("fShowerMinEndZ/F:fShowerMaxVtxZ/F:fShowerMaxEndZ/F:fShowerVtxDist/F:fShowerPt/F:");
299  paramst.append("fSliceMinPlane/F:fVertexX/F:fVertexY/F:fVertexZ/F:fVertexMaxW/F:fCalorimetricE/F:fNueEnergy/F:");
300  paramst.append("fMinWestDist/F:fMinEastDist/F:fMinBotDist/F:fMinTopDist/F:fMinFrontDist/F:fMinBackDist:");
301  paramst.append("fLIDVal/F:fLEMVal/F:fCVNVal/F");
302  fNueTree = tfs->make<TTree>("Nue", "Skimmed Nue Interaction Cut Parameters");
303  fNueTree->Branch("info", &fEventInfo, infost .c_str());
304  fNueTree->Branch("params", &fNueCutParams, paramst.c_str());
305  fNueTree->Branch("sr", &fNueSRParams, paramst.c_str());
306 
307  }
308 
309  // scan all the parameters in an interactive ROOT session
310  // Numu->Scan("info.run:info.event:params.fBestPIDTrack:params.fNeutrinoE:params.fNDHadronicCal:params.fReMIdVal:params.fQePIDVal:params.fQePIDNTracks:params.fSliceHits:params.fSliceContigPlanes:params.fSliceMaxPlane:params.fSliceMinPlane:params.fSliceCellsFromEdge:params.fPlanesToFront:params.fPlanesToBack:params.fCosRejKalFwdCell:params.fCosRejKalBakCell:params.fCosRejAngleKal:params.fCosRejCosBakCell:params.fCosRejCosFwdCell:params.fCosRejConCosPID:params.fCosRejKalBakCellND:params.fCosRejKalFwdCellND:params.fCosRejKalYPosAtTrans:params.fTrackStartZ:params.fTrackStopZ:params.fTrackCCE");
311 
312  // compare the values from skimming vs CAF
313  // Numu->Scan("info.run:info.event:params.fBestPIDTrack/sr.fBestPIDTrack:params.fNeutrinoE/sr.fNeutrinoE:params.fNDHadronicCal/sr.fNDHadronicCal:params.fReMIdVal/sr.fReMIdVal:params.fQePIDVal/sr.fQePIDVal:params.fQePIDNTracks/sr.fQePIDNTracks:params.fSliceHits/sr.fSliceHits:params.fSliceContigPlanes/sr.fSliceContigPlanes:params.fSliceMaxPlane/sr.fSliceMaxPlane:params.fSliceMinPlane/sr.fSliceMinPlane:params.fSliceCellsFromEdge/sr.fSliceCellsFromEdge:params.fPlanesToFront/sr.fPlanesToFront:params.fPlanesToBack/sr.fPlanesToBack:params.fCosRejKalFwdCell/sr.fCosRejKalFwdCell:params.fCosRejKalBakCell/sr.fCosRejKalBakCell:params.fCosRejAngleKal/sr.fCosRejAngleKal:params.fCosRejCosBakCell/sr.fCosRejCosBakCell:params.fCosRejCosFwdCell/sr.fCosRejCosFwdCell:params.fCosRejConCosPID/sr.fCosRejConCosPID:params.fCosRejKalBakCellND/sr.fCosRejKalBakCellND:params.fCosRejKalFwdCellND/sr.fCosRejKalFwdCellND:params.fCosRejKalYPosAtTrans/sr.fCosRejKalYPosAtTrans:params.fTrackStartZ/sr.fTrackStartZ:params.fTrackStopZ/sr.fTrackStopZ:params.fTrackCCE/sr.fTrackCCE");
314 
315 // Nue->Scan("info.run:info.event:params.fNumPlane:params.fCellsPerPlane:fNumXCell:params.fNumYCell:params.fSliceHits:params.fHitAsymmetry:params.fCosShowers:params.fShowerDirX:params.fShowerDirY:params.fShowerDirZ:params.fShowerLength:params.fFracShowerHits:params.fShowerVertexX:params.fShowerVertexY:params.fShowerVertexZ:params.fShowerEndX:params.fShowerEndY:params.fShowerEndZ:params.fShowerMaxW:params.fShowerMaxW:params.fShowerMinVtxZ:params.fShowerMinEndZ:params.fShowerMaxVtxZ:params.fShowerMaxEndZ:params.fShowerVtxDist:params.fShowerPt:params.fVertexX:params.fVertexY:params.fVertexZ:params.fVertexMaxW:params.fCalorimetricE:params.fMinWestDist:params.fMinEastDist:params.fMinBotDist:params.fMinTopDist:params.fMinFrontDist:params.fMinBackDist");
316 
317 
318  return;
319 }
320 
321 //------------------------------------------------------------------------------
323 {
324  fDetId = fDetService->DetId();
325 
326  return;
327 }
328 
329 //------------------------------------------------------------------------------
331 {
332  // get the CosmicExposure and POTSum out of the subrun
333  //art::Handle<sumdata::CosmicExposure> ce;
334  //art::Handle<sumdata::POTSum> pot;
335  //sr.getByLabel(fCosmicExposureLabel, ce );
336  //sr.getByLabel(fPOTSumLabel, pot);
337 
338 // if( ce.isValid() ){
339 // MF_LOG_DEBUG("SkimmerAna") << sr.run() << " " << ce->totlivetime;
340 // fCosmicExposureByRun->Fill(sr.run(), ce->totlivetime);
341 // }
342 // else MF_LOG_DEBUG("SkimmerAna") << "No cosmic exposure found";
343 // if( pot.isValid() ){
344 // MF_LOG_DEBUG("SkimmerAna") << sr.run() << " " << pot->totgoodpot;
345 // fPOTByRun->Fill(sr.run(), 1.e-16 * pot->totgoodpot);
346 // }
347 // else MF_LOG_DEBUG("SkimmerAna") << "No pot summary found";
348 
349  std::vector<art::Handle<sumdata::CosmicExposure>> ce;
350  std::vector<art::Handle<sumdata::POTSum>> pot;
351  sr.getManyByType(ce );
352  sr.getManyByType(pot);
353  if( ce.size() > 0 ){
354  MF_LOG_DEBUG("SkimmerAna") << sr.run() << " " << ce[0]->totlivetime;
355  fCosmicExposureByRun->Fill(sr.run(), ce[0]->totlivetime);
356  }
357  else MF_LOG_DEBUG("SkimmerAna") << "No cosmic exposure found";
358  if( pot.size() > 0 ){
359  MF_LOG_DEBUG("SkimmerAna") << sr.run() << " " << pot[0]->totgoodpot;
360  fPOTByRun->Fill(sr.run(), 1.e-3 * pot[0]->totpot);
361  fGoodPOTByRun->Fill(sr.run(), 1.e-3 * pot[0]->totgoodpot);
362  }
363  else MF_LOG_DEBUG("SkimmerAna") << "No pot summary found";
364 
365 
366  return;
367 }
368 
369 //------------------------------------------------------------------------------
371 {
372  fEventInfo.run = e.run();
373  fEventInfo.subrun = e.subRun();
374  fEventInfo.event = e.id().event();
375 
376  // MF_LOG_DEBUG("SkimmerAna")
377  MF_LOG_VERBATIM("SkimmerAna") <<" ===========>>>>> SkimmerAna "
378  << e.run() << " " << e.subRun() << " " << e.id().event();
379 
380  if( !fNumuInstance.empty() ) this->FillNumuHistograms(e);
381  if( !fNueInstance.empty() ) this->FillNueHistograms(e);
382 
383  return;
384 }
385 
386 //------------------------------------------------------------------------------
388 {
389 
390  fSkimmedLabel = p.get<std::string >("SkimmedLabel");
391  fCosmicExposureLabel = p.get<std::string >("CosmicExposureLabel", "ifdbspillinfo");
392  fPOTSumLabel = p.get<std::string >("POTSumLabel", "exposure" );
393  fCompareNumuCuts = p.get<bool >("CompareNumuCuts", false);
394  fCompareNueCuts = p.get<bool >("CompareNueCuts", false);
395  auto instanceLabels = p.get<std::vector<std::string> >("Instances");
396 
397  for(auto itr : instanceLabels){
398  if(itr.find("numu") != std::string::npos) fNumuInstance = itr;
399  if(itr.find("nue") != std::string::npos) fNueInstance = itr;
400  }
401 
402  return;
403 }
404 
405 //------------------------------------------------------------------------------
407 {
408  // get the relevant data products out of the event record
409  // need the slice, the vertex for the slice, and the showers in the slice
411  auto sliceHandle = e.getValidHandle< std::vector<rb::Cluster> >(tag);
412  if( sliceHandle.failedToGet() ){
413  MF_LOG_WARNING("SkimmerAna") << "No slices passing nue selection in this event";
414  return;
415  }
416 
417  art::FindMany<rb::Shower> fms (sliceHandle, e, tag);
418  art::FindMany<rb::Prong> fmp (sliceHandle, e, tag);
419  art::FindMany<rb::Vertex> fmv (sliceHandle, e, tag);
420  art::FindMany<cosrej::NueCosRej> fmc (sliceHandle, e, tag);
421  art::FindMany<slid::EventLID> fmlid (sliceHandle, e, tag);
422  art::FindMany<lem::PIDDetails> fmlem (sliceHandle, e, tag);
423  art::FindMany<cvn::Result> fmcvn (sliceHandle, e, tag);
424  art::FindMany<presel::Veto> fmnueveto(sliceHandle, e, tag);
425 
426 
428 
429  for(size_t s = 0; s < sliceHandle->size(); ++s){
430 
431  std::vector<rb::Shower const*> showers;
432  std::vector<rb::Prong const*> prongs;
433  std::vector<rb::Vertex const*> vertices;
434  std::vector<slid::ShowerLID const*> slids;
435  std::vector<slid::EventLID const*> elids;
436  std::vector<lem::PIDDetails const*> lems;
437  std::vector<cvn::Result const*> cvns;
438  std::vector<presel::Veto const*> nuevetos;
439  std::vector<cosrej::NueCosRej const*> cosrejs;
440 
441  su.FillSliceVector(s, fmv, vertices);
442  su.FillSliceVector(s, fmc, cosrejs);
443  su.FillSliceVector(s, fms, showers);
444  su.FillSliceVector(s, fmp, prongs);
445  su.FillSliceVector(s, fmlid, elids);
446  su.FillSliceVector(s, fmlem, lems);
447  su.FillSliceVector(s, fmcvn, cvns);
448  su.FillSliceVector(s, fmnueveto, nuevetos);
449  su.FillSliceVector(e, tag, showers, slids);
450 
451  if(vertices.size() != 1 ||
452  cosrejs .size() != 1 ||
453  lems .size() != 1 ||
454  cvns .size() != 1 ||
455  elids .size() != 1 ||
456  showers .size() < 1 ||
457  prongs .size() < 1 ||
458  slids .size() != showers.size()){
459  MF_LOG_DEBUG("SkimmerAna")
460  << "Wrong number of vertices: " << vertices.size()
461  << " : 1 or cosrejs: " << cosrejs .size()
462  << " : 1 or lids: " << elids .size()
463  << " : 1 or lems: " << lems .size()
464  << " : 1 or cvns: " << cvns .size()
465  << " : 1 or showers: " << showers .size()
466  << " : or prongs: " << prongs .size()
467  << " or slids: " << slids .size()
468  << " for the skimmed slice, that should never happen";
469  continue;
470  }
471 
472  // make the parameters object
473  skim::ParametersNue params((*sliceHandle)[s],
474  *(vertices[0]),
475  *(cosrejs[0]),
476  showers,
477  prongs,
478  slids,
479  *(elids[0]),
480  *(lems[0]),
481  *(cvns[0]),
482  *(nuevetos[0]));
483 
484 
485  //skim::ParametersNue sr(fosr.at(s).ref());
486 
487  fNueCutParams = params.ParametersStruct();
488  //fNueSRParams = sr .ParametersStruct();
489 
490  //if (fCompareNueCuts) params.Compare(fNueSRParams);
491 
492  fLEM ->Fill(params.LEMVal());
493  fCVN ->Fill(params.CVNVal());
494  fLID ->Fill(params.LIDVal());
495  fCalorimetricEnergy->Fill(params.CalorimetricE());
496 
498  fFDSliceHitsPerPlane->Fill(params.CellsPerPlane());
499  fFDSliceHitsXView ->Fill(params.ShowerNumXCell());
500  fFDSliceHitsYView ->Fill(params.ShowerNumYCell());
501  fFDSliceHitsViewAsym->Fill(params.HitAsymmetry());
502  fFDSliceCosAngleShw ->Fill(params.CosShowers());
503  fFDSliceDistShwVtx ->Fill(params.ShowerVtxDist());
504  fFDSliceFracShwHits ->Fill(params.FracShowerHits());
505  fFDVtxDistToEastWall->Fill(params.MinEastDist());
506  fFDVtxDistToWestWall->Fill(params.MinWestDist());
507  fFDVtxDistToTop ->Fill(params.MinTopDist());
508  fFDVtxDistToBottom ->Fill(params.MinBotDist());
509  fFDVtxDistToFront ->Fill(params.MinFrontDist());
510  fFDVtxDistToBack ->Fill(params.MinBackDist());
511  fFDPt ->Fill(params.ShowerPt());
512  }
513 
515  fNDVtxX ->Fill(params.Vertex().X() );
516  fNDVtxY ->Fill(params.Vertex().Y() );
517  fNDVtxZ ->Fill(params.Vertex().Z() );
518  fNDShwVtxX ->Fill(params.ShowerVertex().X());
519  fNDShwVtxY ->Fill(params.ShowerVertex().Y());
520  fNDShwVtxZ ->Fill(params.ShowerVertex().Z());
521  fNDShwEndX ->Fill(params.ShowerEnd().X() );
522  fNDShwEndY ->Fill(params.ShowerEnd().Y() );
523  fNDShwEndZ ->Fill(params.ShowerEnd().Z() );
524  fNDSliceHits ->Fill(params.SliceHits() );
525  fNDCalorimetricE ->Fill(params.CalorimetricE() );
526  fNDShowerLength ->Fill(params.ProngLength() );
527  fNDSliceNHitsPerPlane->Fill(params.CellsPerPlane() );
528  fNDPlanesToFront ->Fill(params.PlanesToFront() );
529  fNDNShowers ->Fill(params.NumShowers() );
530  fNDShowerDistToVtx ->Fill(params.ShowerVtxDist() );
531  fNDShowerCalE ->Fill(params.ShowerCalE() );
532  fNDShowerCalEfraction->Fill(params.ShowerCalE() / params.CalorimetricE());
533  fNDSliceNoShowerCalE ->Fill(params.CalorimetricE() - params.ShowerCalE() );
534  fNDCalEperHit ->Fill(params.CalorimetricE() / params.SliceHits() );
535  }
536 
537  fNueTree->Fill();
538 
539  } // end loop over skimmed slices
540 
541  return;
542 }
543 
544 //------------------------------------------------------------------------------
546 {
548 
549  // get the relevant data products out of the event record
550  // need the slice, rb::Energy, remid::ReMId, cosrej::CosRejObj,
551  // qeef::QePId, numue::NumuE, rb::Track
552  auto sliceHandle = e.getValidHandle< std::vector<rb::Cluster> >(inputTag);
553  if( sliceHandle.failedToGet() ){
554  MF_LOG_WARNING("SkimmerAna") << "No slices passing numu selection in this event";
555  return;
556  }
557 
558  MF_LOG_VERBATIM("SkimmerAna") << "Start Filling Numu Histos ... " ;
559 
560 
561  art::FindManyP<rb::Track> fmt (sliceHandle, e, inputTag);
562  art::FindMany<rb::Energy> fme (sliceHandle, e, inputTag);
563  art::FindMany<cosrej::CosRejObj> fmcro (sliceHandle, e, inputTag);
564  art::FindMany<rb::Shower> fms (sliceHandle, e, inputTag);
565  art::FindMany<rb::Prong> fmp (sliceHandle, e, inputTag);
566  art::FindMany<rb::Vertex> fmv (sliceHandle, e, inputTag);
567  art::FindMany<rb::Track> fmtks (sliceHandle, e, inputTag);
568  art::FindMany<cosrej::NueCosRej> fmnuecro(sliceHandle, e, inputTag); // Nue CosRej added for Analysis 2017
569  art::FindMany<qeef::QePId> fmqe (sliceHandle, e, inputTag);
570  art::FindMany<cvn::Result> fmcvn (sliceHandle, e, inputTag);
571  art::FindMany<cvn::Result> fmocvn (sliceHandle, e, inputTag);
572  art::FindMany<numue::NumuE> fmnue (sliceHandle, e, inputTag);
573 
574 
576 
577  size_t bestPIDTrack = std::numeric_limits<size_t>::max();
578 
579  for(size_t s = 0; s < sliceHandle->size(); ++s){
580 
581  std::vector< art::Ptr<rb::Track> > trackPtrs;
582  std::vector<cosrej::CosRejObj const*> cosrej;
583  std::vector<qeef::QePId const*> qepid;
584  std::vector<cvn::Result const*> cvns;
585  std::vector<cvn::Result const*> cvn2017s;
586  std::vector<numue::NumuE const*> numue;
587 
588  std::vector<cosrej::NueCosRej const*> nuecosrej; // Nue CosRej added for Analysis 2017
589  std::vector<rb::Track const*> tracks; /// 2017
590  std::vector<rb::Shower const*> showers; /// 2017
591  std::vector<rb::Prong const*> prongs; /// 2017
592  std::vector<rb::Vertex const*> vertices; /// 2017
593  std::vector<slid::ShowerLID const*> slids; /// 2017
594 
595  fmt.get(s, trackPtrs);
596 
597  // better have at least one track
598  if(trackPtrs.size() < 1 ) continue;
599 
600  // get the best pid for the slice - this index will not be the
601  // same as in the preskimmed slice
602  bestPIDTrack = remid::HighestPIDTrack(trackPtrs, inputTag, e);
603 
604  fmcro.get(s, cosrej);
605  fmnue.get(s, numue);
606  fmqe .get(s, qepid);
607 
608  fmnuecro.get(s, nuecosrej); // Nue CosRej added for Analysis 2017
609  fmcvn.get(s, cvns);
610  su.FillSliceVector(s, fmtks, tracks);
611  su.FillSliceVector(s, fmv, vertices);
612  su.FillSliceVector(s, fms, showers);
613  su.FillSliceVector(s, fmp, prongs);
614  su.FillSliceVector(e, inputTag, showers, slids);
615 
616  if(vertices.size() != 1 ||
617  cvns .size() != 1 ||
618  cvn2017s.size() != 1 ||
619  showers .size() < 1 ||
620  prongs .size() < 1 ||
621  tracks .size() < 1 ||
622  slids .size() != showers.size()){
623  MF_LOG_DEBUG("SkimmerAna")
624  << "Wrong number of vertices: " << vertices.size()
625  << " : 1 or cvns: " << cvns .size()
626  << " : 1 or cvn2017s: " << cvn2017s.size()
627  << " : 1 or showers: " << showers .size()
628  << " : or prongs: " << prongs .size()
629  << " : or tracks: " << tracks .size()
630  << " or slids: " << slids .size()
631  << " for the skimmed slice, that should never happen";
632  continue;
633  }
634 
635  qeef::QePId qep;
636  if( qepid.size() > 0 ) qep = *(qepid.front());
637 
638  // this is not terribly efficient, as we are making the FindOneP
639  // for every slice in the event, but I don't see a good way around
640  // that while still getting the bookkeeping right
641  art::FindOne<rb::Energy> foe (trackPtrs, e, inputTag);
642  art::FindOne<remid::ReMId> foid(trackPtrs, e, inputTag);
643 
644  MF_LOG_VERBATIM("SkimmerAna") << "skim::ParametersNumu params ...";
645  skim::ParametersNumu params(bestPIDTrack,
646  foe.at(bestPIDTrack).ref(),
647  foid.at(bestPIDTrack).ref(),
648  *(cosrej.front()),
649  *(nuecosrej.front()),
650  *(vertices[0]),
651  tracks,
652  showers,
653  prongs,
654  slids,
655  qep,
656  *(numue.front()),
657  sliceHandle->at(s),
658  *(trackPtrs[bestPIDTrack]),
659  *(cvns[0]),
660  *(cvn2017s[0]),
661  e.isRealData());
662 
663 // skim::ParametersNumu sr(fosr.at(s).ref(),
664 // e.isRealData());
665 
666  fNumuCutParams = params.ParametersStruct();
667 // fNumuSRParams = sr .ParametersStruct();
668 //
669 // if (fCompareNumuCuts) params.Compare(fNumuSRParams);
670 
671  // Fill sanity check histograms
672  fNeutrinoE ->Fill(params.NeutrinoE());
673  fReMId ->Fill(params.ReMIdValue());
674  fSliceHits ->Fill(params.SliceHits());
675  fSliceContigPlanes->Fill(params.SliceContigPlanes());
676  fCellsFromEdge ->Fill(params.SliceCellsFromEdge());
677  if(params.QePIDNTracks() == 1) fQePId1Track->Fill(params.QePIDValue());
678  if(params.QePIDNTracks() == 2) fQePId2Track->Fill(params.QePIDValue());
679 
681  fFDKalmanFwdCell->Fill(params.CosRejKalFwdCell());
682  fFDKalmanBakCell->Fill(params.CosRejKalBakCell());
683  fFDKalmanAngle ->Fill(params.CosRejAngleKal());
684  fFDCosBakCell ->Fill(params.CosRejCosBakCell());
685  fFDCosFwdCell ->Fill(params.CosRejCosFwdCell());
686  fFDContPID ->Fill(params.CosRejConCosPID());
687  fFDPlanesToBack ->Fill(params.PlanesToBack());
688  fFDPlanesToFront->Fill(params.PlanesToFront());
689  fFDSliceHits ->Fill(params.SliceHits());
690  }
691 
693  fNDMaxPlaneHit ->Fill(params.SliceMaxPlane());
694  fNDMinPlaneHit ->Fill(params.SliceMinPlane());
695  fNDTrackStartZ ->Fill(params.TrackStartZ());
696  fNDTrackStopZ ->Fill(params.TrackStopZ());
697  fNDKalmanBakCell ->Fill(params.CosRejKalBakCellND());
698  fNDKalmanFwdCell ->Fill(params.CosRejKalFwdCellND());
699  fNDHadCalEnergySum->Fill(params.NDHadronicCal());
700  }
701 
702  MF_LOG_VERBATIM("SkimmerAna") << " ==============>>>>>> SkimmerAna Fill NumuTree ...";
703  fNumuTree->Fill();
704 
705  }// end loop over slices
706 
707  return;
708 }
709 
void beginJob() override
float const & PlanesToFront() const
Definition: ParametersNue.h:70
novadaq::cnv::DetId DetId() const
What detector are we in?
float const & CalorimetricE() const
float const & ShowerNumYCell() const
Definition: ParametersNue.h:78
void beginRun(art::Run const &r) override
Energy estimators for CC events.
Definition: FillEnergies.h:7
float const & ShowerPt() const
Definition: ParametersNue.h:88
float const & CVNVal() const
SkimmerAna(fhicl::ParameterSet const &p)
TTree * fNumuTree
Tree to hold values of cut parameters for selected events.
NueCutParameters fNueCutParams
from the ParametersNumu object
float const & ProngLength() const
Definition: ParametersNue.h:80
TVector3 const & ShowerEnd() const
Definition: ParametersNue.h:85
void reconfigure(fhicl::ParameterSet const &p)
const char * p
Definition: xmltok.h:285
SkimmerAna & operator=(SkimmerAna const &)=delete
size_type get(size_type i, reference item, data_reference data) const
Definition: FindMany.h:455
TFile * fmc
Definition: bdt_com.C:14
void FillNueHistograms(art::Event const &e)
float const & MinEastDist() const
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:27
DEFINE_ART_MODULE(TestTMapFile)
Module to create a summary of total POT seen in a job.
Definition: Evaluator.h:27
std::string fNueInstance
label for nue products
novadaq::cnv::DetId fDetId
the id of the detector we are looking at
float const & CellsPerPlane() const
Definition: ParametersNue.h:71
Definition: Run.h:21
void FillSliceVector(size_t const &s, art::FindMany< T > &fmp, std::vector< T const * > &ptrMap)
EventInfo fEventInfo
keep track of run and subrun, and event
TVector3 const & Vertex() const
Definition: ParametersNue.h:98
NumuCutParameters fNumuCutParams
from the ParametersNumu object
float const & ShowerCalE() const
Definition: ParametersNue.h:81
const XML_Char * s
Definition: expat.h:262
bool isRealData() const
Far Detector at Ash River, MN.
Result for CVN.
TTree * fNueTree
Tree to hold values of cut parameters for selected events.
std::string fNumuInstance
label for numu products
float const & MinTopDist() const
float const & HitAsymmetry() const
Definition: ParametersNue.h:75
T get(std::string const &key) const
Definition: ParameterSet.h:231
TVector3 const & ShowerVertex() const
Definition: ParametersNue.h:86
#define pot
NueCutParameters const & ParametersStruct() const
Near Detector in the NuMI cavern.
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
Definition: DataViewImpl.h:491
void getManyByType(std::vector< Handle< PROD >> &results) const
Definition: DataViewImpl.h:558
std::string fPOTSumLabel
label of module creating pot sum
caf::StandardRecord * sr
SubRunNumber_t subRun() const
float const & MinBackDist() const
RunNumber_t run() const
Vertex location in position and time.
float const & LIDVal() const
art::ServiceHandle< ds::DetectorService > fDetService
detector service
void FillNumuHistograms(art::Event const &e)
void beginSubRun(art::SubRun const &sr) override
EventNumber_t event() const
Definition: EventID.h:116
float const & CosShowers() const
Definition: ParametersNue.h:79
Cosmic Rejection PIDs for Numu analysis.
Definition: FillParentInfo.h:9
std::string fSkimmedLabel
label of module creating skimmed file
#define MF_LOG_VERBATIM(category)
#define MF_LOG_DEBUG(id)
float const & ShowerNumXCell() const
Definition: ParametersNue.h:77
NueCutParameters fNueSRParams
from the ParametersNumu object
TRandom3 r(0)
bool fCompareNumuCuts
bool for comparing the numu cuts
float const & MinBotDist() const
size_type get(size_type i, reference item, data_reference data) const
Definition: FindManyP.h:455
float const & SliceHits() const
Definition: ParametersNue.h:74
float const & FracShowerHits() const
Definition: ParametersNue.h:83
#define MF_LOG_WARNING(category)
float const & LEMVal() const
Float_t e
Definition: plot.C:35
std::string fCosmicExposureLabel
label of module creating cosmic exposure
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:66
float const & MinFrontDist() const
float const & NumShowers() const
Definition: ParametersNue.h:76
EventID id() const
unsigned int HighestPIDTrack(const std::vector< art::Ptr< rb::Track > > &sliceTracks, const std::string &remidModuleLabel, const art::Event &e)
Definition: ReMId.cxx:249
float const & ShowerVtxDist() const
Definition: ParametersNue.h:87
NumuCutParameters fNumuSRParams
from the ParametersNumu object
float const & MinWestDist() const
bool fCompareNueCuts
bool for comparing the nue cuts
void analyze(art::Event const &e) override
enum BeamMode string