CosmicBeamComparison_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // Class: CosmicBeamComparison
3 // Module Type: analyzer
4 // File: CosmicBeamComparison_module.cc
5 //
6 // Medbh Campbell: This is essentially a copy of SkimmerAna_module.cc but with some
7 // changes/additions, namely a couple of trees have been added to the output which
8 // contain useful information for analysing cosmic rejection cuts (by comparing
9 // cosmic data and beam mc)
10 ////////////////////////////////////////////////////////////////////////
11 
23 #include "fhiclcpp/ParameterSet.h"
27 
28 #include <memory>
29 #include <set>
30 #include <limits>
31 
32 #include "TH1.h"
33 #include "TH2.h"
34 #include "TTree.h"
35 
36 #include "NumuSandbox/NumuSandObj.h"//included to get X/Y/Z min/max 2
37 #include "RecoBase/Cluster.h"
38 #include "RecoBase/Prong.h"
39 #include "RecoBase/Shower.h"
40 #include "RecoBase/Vertex.h"
41 #include "CosRej/NueCosRej.h"
42 #include "CosRej/CosRejObj.h"
43 #include "ShowerLID/ShowerLID.h"
44 #include "ShowerLID/EventLID.h"
45 #include "LEM/PIDDetails.h"
46 #include "CVN/func/Result.h"
47 #include "Preselection/Veto.h"
49 #include "SummaryData/SpillData.h"
50 #include "SummaryData/POTSum.h"
52 #include "NovaDAQConventions/DAQConventions.h"
58 #include "Geometry/Geometry.h"
59 
60 namespace skim {
61  class CosmicBeamComparison;
62 }
63 
65 public:
66  explicit CosmicBeamComparison(fhicl::ParameterSet const & p);
67  virtual ~CosmicBeamComparison();
68 
69  // Plugins should not be copied or assigned.
74 
75  // Required functions.
76  void analyze(art::Event const& e) override;
77 
78  // Selected optional functions.
79  void reconfigure(fhicl::ParameterSet const & p);
80 
81  void beginJob() override;
82  void beginRun(art::Run const& r) override;
83  void beginSubRun(art::SubRun const& sr) override;
84 
85 private:
86 
87  void FillNueHistograms (art::Event const& e);
88  void FillNumuHistograms(art::Event const& e);
89 
90  std::string fSkimmedLabel; ///< label of module creating skimmed file
91  std::string fCosmicExposureLabel; ///< label of module creating cosmic exposure
92  std::string fPOTSumLabel; ///< label of module creating pot sum
93  std::string fNumuInstance; ///< label for numu products
94  std::string fNueInstance; ///< label for nue products
95  bool fCompareNumuCuts; ///< bool for comparing the numu cuts
96  bool fCompareNueCuts; ///< bool for comparing the nue cuts
97  novadaq::cnv::DetId fDetId; ///< the id of the detector we are looking at
100 
101  TTree *fNumuTree; ///< Tree to hold values of cut parameters for selected events
102  TTree *fNueTree; ///< Tree to hold values of cut parameters for selected events
103  TTree *fCosmicTree;///< Tree filled directly from cosrej etc. objects (bypasses parameters struct)
104  TTree *fHitTree;///< Tree to hold information about every hit in a kalman track
105  TTree *fCosmicHitTree;///< Tree to hold information about every hit in a kalman track
106 
107  NumuCutParameters fNumuCutParams; ///< from the ParametersNumu object
108  NueCutParameters fNueCutParams; ///< from the ParametersNumu object
109  NumuCutParameters fNumuSRParams; ///< from the ParametersNumu object
110  NueCutParameters fNueSRParams; ///< from the ParametersNumu object
111  EventInfo fEventInfo; ///< keep track of run and subrun, and event
112 
113  // POT and spill information
114  TH1F *fPOTByRun;
117 
118  // numu sanity check histograms, one for each cut
119  TH1F *fNeutrinoE;
120  TH1F *fReMId;
121  TH1F *fSliceHits;
136  TH1F *fFDKalSpeed;
137  TH1F *fFDCosSpeed;
138  TH1F *fFDDeltaT;
139  TH1F *fFDDeltaX;
140  TH1F *fFDDeltaY;
141  TH1F *fFDDeltaZ;
142  TH1F *fFDSpeed;
144  TH2F *fFDDTvsDY;
146  TH1F *fFDContPID;
147  TH1F *fFDSliceHits;
149 
157 
158  // nue sanity check histograms
172  TH1F* fFDPt;
173 
175  TH1F* fLEM;
176  TH1F* fCVN;
177  TH1F* fLID;
178  TH1F* fNDVtxX;
179  TH1F* fNDVtxY;
180  TH1F* fNDVtxZ;
181  TH1F* fNDNShowers;
182  TH1F* fNDShwVtxX;
183  TH1F* fNDShwVtxY;
184  TH1F* fNDShwVtxZ;
185  TH1F* fNDShwEndX;
186  TH1F* fNDShwEndY;
187  TH1F* fNDShwEndZ;
198 
199  // general sanity check histograms
200  //TH2F*
201 
202  //Values for CosmicTree (cr=cosmic rejection)
203  float crNeutrinoE = 0;
204  double crReMIdVal;
205  float crSliceHits;
214  TVector3 crSliceMinXYZ;
215  TVector3 crSliceMaxXYZ;
222  int crSliceNHitsNearEdge0;//number of slice hits within 0 cells of the edge
232  float crKalAngle;//cosine of angle between kalman track and beam
233  float crCosAngle;//cosine of angle between cosmic track and beam
242  float crConCosPID;
243  TVector3 crTrackStart;//kalman track start
244  TVector3 crTrackStop;//kalman track stop
245  int crTrackHits;//number of hits in the track
246  int crTrackNHitsNearEdge0;//number of track hits within 0 cells of the edge
258  int crTrackXViewHits;//number of hits in xview cells in the track
260  float crTrackTotalADC;//sum of adc of all hits
261  float crTrackTotalGeV;//simple sum of estimated gev of all hits
262  float crTrackCalorimetricEnergy;//simple estimate of neutrino energy
263  float crTrackMinTNS;//minimum time of all hits in track
264  float crTrackMaxTNS;//maximum time of all hits in track
265  TVector3 crTrackMinXYZ;
266  TVector3 crTrackMaxXYZ;
271  float crTrackTotalLength;//total length (cm) of all track segments
272  TVector3 crKalDir;
273  float crTrackStartTime;//time of kalman track start (NOT the hit with the minimum time)
275  TVector3 crCosTrackStart;//cosmic track start
276  TVector3 crCosTrackStop;//cosmic track stop
277  TVector3 crCosTrackDir;
278  int crCosTrackHits;//number of hits in the cosmic track
281  int crCosTrackNHitsNearEdge0;//number of CosTrack hits within 0 cells of the edge
291  int crCosTrackXViewHits;//number of hits in xview cells in the cosmic track
293  float crCosTrackTotalADC;//sum of adc of all hits
294  float crCosTrackTotalGeV;//simple sum of estimated gev of all hits
295  float crCosTrackCalorimetricEnergy;//simple estimate of neutrino energy
296  float crCosTrackMinTNS;//minimum time of all hits in cosmic track
297  float crCosTrackMaxTNS;//maximum time of all hits in cosmic track
304  float crCosTrackTotalLength;//total length (cm) of all cosmic track segments
305  float crCosTrackStartTime;//time of cosmic track start (NOT the hit with the minimum time)
307  float crKalSpeed;//kalman speed from hough fit of hit timings (neg if track travelling in -z dir)
308  float crCosSpeed;
310  float crKalDirScore;//hough fit: chisq difference between +/- c assumption in kalman track
311  float crCosScoreDiff;//hough fit: chisq difference between free fit and +/- c
312  float crKalScoreDiff;//hough fit: chisq difference between free fit and +/- c
313  float crKalSlope;//upmu (least squares) timing fit: fit speed
314  float crCosSlope;
315  float crKalChiSq;//upmu (least squares) timing fit: chisq
316  float crCosChiSq;
317  float crKalChiDiff;//upmu (least squares) timing fit: chisq diff between +/- c assumptions (pos if track is backwards)
319  float crNCloseLiveDiblocks;//to see if it would pass CheckSliceQuality
321  float crScatt;
322  float crQePIDVal;
324  float crXmin2;//second smallest X position of any hit
325  float crYmin2;
326  float crZmin2;
327  float crXmax2;//second smallest X position of any hit
328  float crYmax2;
329  float crZmax2;
330  float crscatllh;//Log-likelihood value from scattering angle (used to calculate ReMId)
331  float crdedxllh;//Log-likelihood value from dE/dx (used to calculate ReMId)
332  float crmeasfrac;//Fraction of planes used to calculate dedxllh
339  TVector3 crTruthNuVertex;
341  int crRun;//information to locate event(s) that pass cuts (useful for event display etc.)
342  int crSubrun;
343  int crEvent;
344 
345  //variables for HitTree with an entry for each hit in Kalman track
346  int htIndex;
348  TVector3 htUncalibXYZ;
349  double htUncalibW;
350  float htCalibTNS;
351  TVector3 htCalibXYZ;
352  float htCalibPE;
355  TVector3 htDir;
356  TVector3 htTrackXYZ;
357  int htPlane;
358  int htCell;
359  int htTrackIndex = 0;
360 
361  int chtIndex;
363  TVector3 chtUncalibXYZ;
364  double chtUncalibW;
367  TVector3 chtDir;
368  TVector3 chtTrackXYZ;
369  int chtPlane;
370  int chtCell;
371 
372 };
373 
374 //------------------------------------------------------------------------------
376 : EDAnalyzer(p)
377 {
378  this->reconfigure(p);
379 
380  return;
381 }
382 
383 //------------------------------------------------------------------------------
385 {
386 }
387 
388 //------------------------------------------------------------------------------
390 {
391  // make histograms for each selection criteria
393 
394  std::string infost("run/I:subrun/I:event/I");
395 
396  fPOTByRun = tfs->make<TH1F>("potByRun", ";Run Number;POT (#times 10^{15})", 10000, 10000., 20000.);
397  fGoodPOTByRun = tfs->make<TH1F>("goodPOTByRun", ";Run Number;POT (#times 10^{15})", 10000, 10000., 20000.);
398  fCosmicExposureByRun = tfs->make<TH1F>("cosmicExposureByRun", ";Run Number;Exposure (s)", 10000, 10000., 20000.);
399 
400  if( !fNumuInstance.empty() ){
401  fNeutrinoE = tfs->make<TH1F>("neutrinoE", ";E_{#nu} (GeV);Events", 100, 0., 10. );
402  fReMId = tfs->make<TH1F>("remid", ";ReMId;Events", 200, -1., 1. );
403  fSliceHits = tfs->make<TH1F>("sliceHits", ";Hits In Slice;Events", 100, 0., 400. );
404  fSliceContigPlanes = tfs->make<TH1F>("sliceContigPlanes", ";Contiguous Planes in Slice;Events", 50, 0., 50. );
405  fCellsFromEdge = tfs->make<TH1F>("cellsFromEdge", ";Cells From Edge;Events", 50, 0., 50. );
406  fQePId1Track = tfs->make<TH1F>("qepid1Track", ";QePId_{1 track};Events", 200, -1., 1. );
407  fQePId2Track = tfs->make<TH1F>("qepid2Track", ";QePId_{2 tracks};Events", 200, -1., 1. );
408  fFDKalmanFwdCell = tfs->make<TH1F>("fdKalmanFwdCell", ";Kalman Foward Cell;Events", 105, -5., 100. );
409  fFDKalmanBakCell = tfs->make<TH1F>("fdKalmanBakCell", ";Kalman Backward Cell;Events", 105, -5., 100. );
410  fFDCosFwdCell = tfs->make<TH1F>("fdCosFwdCell", ";cos Foward Cell;Events", 105, -5., 100. );
411  fFDCosBakCell = tfs->make<TH1F>("fdCosBakCell", ";cos Backward Cell;Events", 105, -5., 100. );
412  fFDPlanesToFront = tfs->make<TH1F>("fdPlanesToFront", ";Planes to Front;Events", 100, 0., 100. );
413  fFDPlanesToBack = tfs->make<TH1F>("fdPlanesToBack", ";Planes to Back;Events", 100, 0., 100. );
414  fFDKalmanAngle = tfs->make<TH1F>("fdKalmanAngle", ";Kalman Angle;Events", 200, -4., 4. );
415  fFDContPID = tfs->make<TH1F>("fdContPID", ";ContPID;Events", 200, -1., 1. );
416  fFDSliceHits = tfs->make<TH1F>("fdSliceHits", ";Slice Hits;Events", 500, 0., 500. );
417  fFDKalmanDirX = tfs->make<TH1F>("fdKalmanDirX", ";X component of Kalman Direction;Events", 200, -1.1,1.1 );
418  fFDKalmanDirY = tfs->make<TH1F>("fdKalmanDirY", ";Y component of Kalman Direction;Events", 200, -1.1,1.1 );
419  fFDKalmanDirZ = tfs->make<TH1F>("fdKalmanDirZ", ";Z component of Kalman Direction;Events", 200, -1.1,1.1 );
420  fFDKalSpeed = tfs->make<TH1F>("fdKalSpeed", ";Inverse speed from timing fit of Kalman track (ns/cm);Events", 200, -1.1, 1.1);
421  fFDCosSpeed = tfs->make<TH1F>("fdCosSpeed", ";Inverse speed from timing fit of Cosmic track (ns/cm);Events", 200, -1.1, 1.1);
422  fFDDeltaT = tfs->make<TH1F>("fdDeltaT", ";Diff between max and min hit time (ns);Events", 200, -1000,1000 );
423  fFDDeltaX = tfs->make<TH1F>("fdDeltaX", ";Diff between start and end X pos (cm);Events", 400, -2000,2000 );
424  fFDDeltaY = tfs->make<TH1F>("fdDeltaY", ";Diff between start and end Y pos (cm);Events", 400, -2000,2000 );
425  fFDDeltaZ = tfs->make<TH1F>("fdDeltaZ", ";Diff between start and end Z pos (cm);Events", 700, -1000,6000 );
426  fFDSpeed = tfs->make<TH1F>("fdSpeed", ";Speed (m/s);Events", 1000, 0, 10000000 );
427  fFDTrackLength = tfs->make<TH1F>("fdTracklength", ";Tracklength (cm);Events", 500, 0, 5000 );
428  fFDDTvsDY = tfs->make<TH2F>("fdDtVsDY", ";Diff between max and min hit time (ns);Diff between start and end Y pos (cm)", 200, 0, 1000, 400, -2000, 2000);
429  fFDCosSpeedVsScore = tfs->make<TH2F>("fdCosSpeedVsScore", ";Inverse speed of Cosmic track (ns/cm);Hough fit score in forward dir - reverse dir", 200, -0.1, 0.1, 200, -1, 1);
430  fFDDiffLastFirstLivePlanes = tfs->make<TH1F>("fddifflastfirstliveplane", ";Last good diblock in track's section of mask - first", 16, -1, 15);
431 
432  fNDMaxPlaneHit = tfs->make<TH1F>("ndMaxPlaneHit", ";Max Hit Plane;Events", 220, 0., 220. );
433  fNDMinPlaneHit = tfs->make<TH1F>("ndMinPlaneHit", ";Min Hit Plane;Events", 220, 0., 220. );
434  fNDTrackStartZ = tfs->make<TH1F>("ndTrackStartZ", ";Track Starting Z Pos (cm);Events", 200, 0., 2000.);
435  fNDTrackStopZ = tfs->make<TH1F>("ndTrackStopZ", ";Track Ending Z Pos (cm);Events", 200, 0., 2000.);
436  fNDKalmanFwdCell = tfs->make<TH1F>("ndKalmanFwdCell", ";Kalman Foward Cell;Events", 100, 0., 100. );
437  fNDKalmanBakCell = tfs->make<TH1F>("ndKalmanBakCell", ";Kalman Backward Cell;Events", 100, 0., 100. );
438  fNDHadCalEnergySum = tfs->make<TH1F>("ndHadCalEnergySum", ";Hadronic Energy (GeV);Events", 100, 0., 10. );
439 
440  std::string paramst("fBestPIDTrack/F:fNeutrinoE/F:fNDHadronicCal/F:");
441  paramst.append("fReMIdVal/D:fQePIDVal/F:fQePIDNTracks/F:");
442  paramst.append("fSliceHits/F:fSliceContigPlanes/F:fSliceMaxPlane/F:fSliceMinPlane/F:fSliceCellsFromEdge/F:");
443  paramst.append("fPlanesToFront/F:fPlanesToBack/F:fCosRejKalFwdCell/F:fCosRejKalBakCell/F:fCosRejAngleKal/F:");
444  paramst.append("fCosRejCosBakCell/F:fCosRejCosFwdCell/F:fCosRejConCosPID/F:fCosRejKalBakCellND/F:fCosRejKalFwdCellND/F:");
445 
446  std::string srst = paramst;
447  srst.append("fCosRejKalYPosAtTrans/F:fTrackStartZ/F:fTrackStopZ/F:fTrackCCE/F");
448 
449  paramst.append("fCosRejKalYPosAtTrans/F:fTrackStart:fTrackStop:fTrackCCE/F:");
450  paramst.append("fTrackStartTime/F:fTrackStopTime/F:");
451  paramst.append("fCosRejKalSpeed/F:fCosRejCosSpeed/F:fCosRejCosDirScore/F:");
452  paramst.append("fDiffLastFirstLiveDiblock/F");
453  fNumuTree = tfs->make<TTree>("Numu", "Skimmed Numu Interaction Cut Parameters");
454  fNumuTree->Branch("info", &fEventInfo, infost .c_str());
455  fNumuTree->Branch("params", &fNumuCutParams, paramst.c_str());
456  fNumuTree->Branch("sr", &fNumuSRParams, srst.c_str());//paramst.c_str());
457  }
458 
459  fCosmicTree = tfs->make<TTree>("FDCosmics", "Information for cosmic rej");
460  fCosmicTree->Branch("NeutrinoE", &crNeutrinoE);
461  fCosmicTree->Branch("ReMId", &crReMIdVal);
462  fCosmicTree->Branch("SliceHits", &crSliceHits);
463  fCosmicTree->Branch("SliceXViewHits", &crSliceXViewHits);
464  fCosmicTree->Branch("SliceYViewHits", &crSliceYViewHits);
465  fCosmicTree->Branch("SliceContigPlanes", &crSliceContigPlanes);
466  fCosmicTree->Branch("SliceMaxPlane", &crSliceMaxPlane);
467  fCosmicTree->Branch("SliceMinPlane", &crSliceMinPlane);
468  fCosmicTree->Branch("SliceTotalADC", &crSliceTotalADC);
469  fCosmicTree->Branch("SliceTotalGeV", &crSliceTotalGeV);
470  fCosmicTree->Branch("SliceCalorimetricEnergy", &crSliceCalorimetricEnergy);
471  fCosmicTree->Branch("SliceMinXYZ", &crSliceMinXYZ);
472  fCosmicTree->Branch("SliceMaxXYZ", &crSliceMaxXYZ);
473  fCosmicTree->Branch("SliceMinTNS", &crSliceMinTNS);
474  fCosmicTree->Branch("SliceMaxTNS", &crSliceMaxTNS);
475  fCosmicTree->Branch("CellsFromEdge", &crCellsFromEdge);
476  fCosmicTree->Branch("SliceClosestEdge", &crSliceClosestEdge);
477  fCosmicTree->Branch("PlanesToFront", &crPlanesToFront);
478  fCosmicTree->Branch("PlanesToBack", &crPlanesToBack);
479  fCosmicTree->Branch("SliceNHitsNearEdge0", &crSliceNHitsNearEdge0);
480  fCosmicTree->Branch("SliceNHitsNearEdge1", &crSliceNHitsNearEdge1);
481  fCosmicTree->Branch("SliceNHitsNearEdge2", &crSliceNHitsNearEdge2);
482  fCosmicTree->Branch("SliceNHitsNearEdge3", &crSliceNHitsNearEdge3);
483  fCosmicTree->Branch("SliceNHitsNearEdge4", &crSliceNHitsNearEdge4);
484  fCosmicTree->Branch("SliceNHitsNearEdge5", &crSliceNHitsNearEdge5);
485  fCosmicTree->Branch("SliceNHitsNearEdge6", &crSliceNHitsNearEdge6);
486  fCosmicTree->Branch("SliceNHitsNearEdge7", &crSliceNHitsNearEdge7);
487  fCosmicTree->Branch("SliceNHitsNearEdge8", &crSliceNHitsNearEdge8);
488  fCosmicTree->Branch("SliceNHitsNearEdge9", &crSliceNHitsNearEdge9);
489  fCosmicTree->Branch("KalAngle", &crKalAngle);
490  fCosmicTree->Branch("CosAngle", &crCosAngle);
491  fCosmicTree->Branch("KalFwdCell", &crKalFwdCell);
492  fCosmicTree->Branch("KalBakCell", &crKalBakCell);
493  fCosmicTree->Branch("CosFwdCell", &crCosFwdCell);
494  fCosmicTree->Branch("CosBakCell", &crCosBakCell);
495  fCosmicTree->Branch("KalFwdDist", &crKalFwdDist);
496  fCosmicTree->Branch("KalBakDist", &crKalBakDist);
497  fCosmicTree->Branch("CosFwdDist", &crCosFwdDist);
498  fCosmicTree->Branch("CosBakDist", &crCosBakDist);
499  fCosmicTree->Branch("ConCosPID", &crConCosPID);
500  fCosmicTree->Branch("KalDir", &crKalDir);
501  fCosmicTree->Branch("NTracks3D", &crNTracks3D);
502  fCosmicTree->Branch("TrackStart", &crTrackStart);
503  fCosmicTree->Branch("TrackStop", &crTrackStop);
504  fCosmicTree->Branch("TrackStartTime", &crTrackStartTime);
505  fCosmicTree->Branch("TrackStopTime", &crTrackStopTime);
506  fCosmicTree->Branch("TrackHits", &crTrackHits);
507  fCosmicTree->Branch("TrackXViewHits", &crTrackXViewHits);
508  fCosmicTree->Branch("TrackYViewHits", &crTrackYViewHits);
509  fCosmicTree->Branch("TrackCellsFromEdge", &crTrackCellsFromEdge);
510  fCosmicTree->Branch("TrackClosestEdge", &crTrackClosestEdge);
511  fCosmicTree->Branch("TrackNHitsNearEdge0", &crTrackNHitsNearEdge0);
512  fCosmicTree->Branch("TrackNHitsNearEdge1", &crTrackNHitsNearEdge1);
513  fCosmicTree->Branch("TrackNHitsNearEdge2", &crTrackNHitsNearEdge2);
514  fCosmicTree->Branch("TrackNHitsNearEdge3", &crTrackNHitsNearEdge3);
515  fCosmicTree->Branch("TrackNHitsNearEdge4", &crTrackNHitsNearEdge4);
516  fCosmicTree->Branch("TrackNHitsNearEdge5", &crTrackNHitsNearEdge5);
517  fCosmicTree->Branch("TrackNHitsNearEdge6", &crTrackNHitsNearEdge6);
518  fCosmicTree->Branch("TrackNHitsNearEdge7", &crTrackNHitsNearEdge7);
519  fCosmicTree->Branch("TrackNHitsNearEdge8", &crTrackNHitsNearEdge8);
520  fCosmicTree->Branch("TrackNHitsNearEdge9", &crTrackNHitsNearEdge9);
521  fCosmicTree->Branch("TrackTotalADC", &crTrackTotalADC);
522  fCosmicTree->Branch("TrackTotalGeV", &crTrackTotalGeV);
523  fCosmicTree->Branch("TrackCalorimetricEnergy", &crTrackCalorimetricEnergy);
524  fCosmicTree->Branch("TrackMinXYZ", &crTrackMinXYZ);
525  fCosmicTree->Branch("TrackMaxXYZ", &crTrackMaxXYZ);
526  fCosmicTree->Branch("TrackEarliestHitPos", &crTrackEarliestHitPos);
527  fCosmicTree->Branch("TrackLatestHitPos", &crTrackLatestHitPos);
528  fCosmicTree->Branch("TrackMinTNS", &crTrackMinTNS);
529  fCosmicTree->Branch("TrackMaxTNS", &crTrackMaxTNS);
530  fCosmicTree->Branch("TrackMaxPlane", &crTrackMaxPlane);
531  fCosmicTree->Branch("TrackMinPlane", &crTrackMinPlane);
532  fCosmicTree->Branch("TrackTotalLength", &crTrackTotalLength);
533  fCosmicTree->Branch("CosTrackStart", &crCosTrackStart);
534  fCosmicTree->Branch("CosTrackStop", &crCosTrackStop);
535  fCosmicTree->Branch("CosTrackStartTime", &crCosTrackStartTime);
536  fCosmicTree->Branch("CosTrackStopTime", &crCosTrackStopTime);
537  fCosmicTree->Branch("CosTrackHits", &crCosTrackHits);
538  fCosmicTree->Branch("CosTrackXViewHits", &crCosTrackXViewHits);
539  fCosmicTree->Branch("CosTrackYViewHits", &crCosTrackYViewHits);
540  fCosmicTree->Branch("CosTrackCellsFromEdge", &crCosTrackCellsFromEdge);
541  fCosmicTree->Branch("CosTrackClosestEdge", &crCosTrackClosestEdge);
542  fCosmicTree->Branch("CosTrackNHitsNearEdge0", &crCosTrackNHitsNearEdge0);
543  fCosmicTree->Branch("CosTrackNHitsNearEdge1", &crCosTrackNHitsNearEdge1);
544  fCosmicTree->Branch("CosTrackNHitsNearEdge2", &crCosTrackNHitsNearEdge2);
545  fCosmicTree->Branch("CosTrackNHitsNearEdge3", &crCosTrackNHitsNearEdge3);
546  fCosmicTree->Branch("CosTrackNHitsNearEdge4", &crCosTrackNHitsNearEdge4);
547  fCosmicTree->Branch("CosTrackNHitsNearEdge5", &crCosTrackNHitsNearEdge5);
548  fCosmicTree->Branch("CosTrackNHitsNearEdge6", &crCosTrackNHitsNearEdge6);
549  fCosmicTree->Branch("CosTrackNHitsNearEdge7", &crCosTrackNHitsNearEdge7);
550  fCosmicTree->Branch("CosTrackNHitsNearEdge8", &crCosTrackNHitsNearEdge8);
551  fCosmicTree->Branch("CosTrackNHitsNearEdge9", &crCosTrackNHitsNearEdge9);
552  fCosmicTree->Branch("CosTrackTotalADC", &crCosTrackTotalADC);
553  fCosmicTree->Branch("CosTrackTotalGeV", &crCosTrackTotalGeV);
554  fCosmicTree->Branch("CosTrackCalorimetricEnergy", &crCosTrackCalorimetricEnergy);
555  fCosmicTree->Branch("CosTrackMinXYZ", &crCosTrackMinXYZ);
556  fCosmicTree->Branch("CosTrackMaxXYZ", &crCosTrackMaxXYZ);
557  fCosmicTree->Branch("CosTrackEarliestHitPos", &crCosTrackEarliestHitPos);
558  fCosmicTree->Branch("CosTrackLatestHitPos", &crCosTrackLatestHitPos);
559  fCosmicTree->Branch("CosTrackMinTNS", &crCosTrackMinTNS);
560  fCosmicTree->Branch("CosTrackMaxTNS", &crCosTrackMaxTNS);
561  fCosmicTree->Branch("CosTrackMaxPlane", &crCosTrackMaxPlane);
562  fCosmicTree->Branch("CosTrackMinPlane", &crCosTrackMinPlane);
563  fCosmicTree->Branch("CosTrackTotalLength", &crCosTrackTotalLength);
564  fCosmicTree->Branch("Scatt", &crScatt);
565  fCosmicTree->Branch("QePIDVal", &crQePIDVal);
566  fCosmicTree->Branch("QePIDNTracks", &crQePIDNTracks);
567  fCosmicTree->Branch("KalSpeed", &crKalSpeed);
568  fCosmicTree->Branch("CosSpeed", &crCosSpeed);
569  fCosmicTree->Branch("KalDirScore", &crKalDirScore);
570  fCosmicTree->Branch("CosDirScore", &crCosDirScore);
571  fCosmicTree->Branch("KalScoreDiff", &crKalScoreDiff);
572  fCosmicTree->Branch("CosScoreDiff", &crCosScoreDiff);
573  fCosmicTree->Branch("KalSlope", &crKalSlope);
574  fCosmicTree->Branch("CosSlope", &crCosSlope);
575  fCosmicTree->Branch("KalChiSq", &crKalChiSq);
576  fCosmicTree->Branch("CosChiSq", &crCosChiSq);
577  fCosmicTree->Branch("KalChiDiff", &crKalChiDiff);
578  fCosmicTree->Branch("CosChiDiff", &crCosChiDiff);
579  fCosmicTree->Branch("NCloseLiveDiblocks", &crNCloseLiveDiblocks);
580  fCosmicTree->Branch("Xmin2", &crXmin2);
581  fCosmicTree->Branch("Ymin2", &crYmin2);
582  fCosmicTree->Branch("Zmin2", &crZmin2);
583  fCosmicTree->Branch("Xmax2", &crXmax2);
584  fCosmicTree->Branch("Ymax2", &crYmax2);
585  fCosmicTree->Branch("Zmax2", &crZmax2);
586  fCosmicTree->Branch("scatllh", &crscatllh);
587  fCosmicTree->Branch("dedxllh", &crdedxllh);
588  fCosmicTree->Branch("measfrac", &crmeasfrac);
589  fCosmicTree->Branch("TruthNNu", &crTruthNNu);
590  fCosmicTree->Branch("TruthNuNHitsSlc", &crTruthNuNHitsSlc);
591  fCosmicTree->Branch("TruthNuNHitsTot", &crTruthNuNHitsTot);
592  fCosmicTree->Branch("TruthNuIsCC", &crTruthNuIsCC);
593  fCosmicTree->Branch("TruthNuVertex", &crTruthNuVertex);
594  fCosmicTree->Branch("TruthNuPDG", &crTruthNuPDG);
595  fCosmicTree->Branch("TruthNuPDGOrig", &crTruthNuPDGOrig);
596  fCosmicTree->Branch("TruthNuVertex", &crTruthNuVertex);
597  fCosmicTree->Branch("TruthNuTime", &crTruthNuTime);
598  fCosmicTree->Branch("Run", &crRun);
599  fCosmicTree->Branch("Subrun", &crSubrun);
600  fCosmicTree->Branch("Event", &crEvent);
601 
602  fHitTree = tfs->make<TTree>("HitTree", "Entry per hit in Kalman track");
603  fHitTree->Branch("Index", &htIndex);
604  fHitTree->Branch("TrackIndex", &htTrackIndex);
605  fHitTree->Branch("UncalibTNS", &htUncalibTNS);
606  fHitTree->Branch("UncalibXYZ", &htUncalibXYZ);
607  fHitTree->Branch("CalibTNS", &htCalibTNS);
608  fHitTree->Branch("CalibXYZ", &htCalibXYZ);
609  fHitTree->Branch("CalibPE", &htCalibPE);
610  fHitTree->Branch("W", &htUncalibW);
611  fHitTree->Branch("DistFromStart", &htDistFromStart);
612  fHitTree->Branch("DistFromStop", &htDistFromStop);
613  fHitTree->Branch("Dir", &htDir);
614  fHitTree->Branch("TrackXYZ", &htTrackXYZ);
615  fHitTree->Branch("Plane", &htPlane);
616  fHitTree->Branch("Cell", &htCell);
617 
618  fCosmicHitTree = tfs->make<TTree>("CosmicHitTree", "Entry per hit in Cosmic track");
619  fCosmicHitTree->Branch("Index", &chtIndex);
620  fCosmicHitTree->Branch("TNS", &chtUncalibTNS);
621  fCosmicHitTree->Branch("XYZ", &chtUncalibXYZ);
622  fCosmicHitTree->Branch("W", &chtUncalibW);
623  fCosmicHitTree->Branch("DistFromStart", &chtDistFromStart);
624  fCosmicHitTree->Branch("DistFromStop", &chtDistFromStop);
625  fCosmicHitTree->Branch("Dir", &chtDir);
626  fCosmicHitTree->Branch("TrackXYZ", &chtTrackXYZ);
627  fCosmicHitTree->Branch("Plane", &chtPlane);
628  fCosmicHitTree->Branch("Cell", &chtCell);
629 
630  if( !fNueInstance.empty() ){
631  fLEM = tfs->make<TH1F>("lem", ";LEM;Slices", 100, 0., 1. );
632  fCVN = tfs->make<TH1F>("cvn", ";CVN;Slices", 100, 0., 1. );
633  fLID = tfs->make<TH1F>("lid", ";LID;Slices", 240, -0.1, 1.1);
634  fCalorimetricEnergy = tfs->make<TH1F>("calorimetricE", ";Calorimetric E (GeV)", 100, 0., 50.);
635 
636  fFDSliceHitsPerPlane = tfs->make<TH1F>("fdSliceHitsPerPlane", ";Slice Hits/Plane;Slices", 100, 0., 50.);
637  fFDSliceHitsXView = tfs->make<TH1F>("fdSliceHitsXView", ";Slice Hits (X view);Slices", 100, 0., 100.);
638  fFDSliceHitsYView = tfs->make<TH1F>("fdSliceHitsYView", ";Slice Hits (Y view);Slices", 100, 0., 100.);
639  fFDSliceHitsViewAsym = tfs->make<TH1F>("fdSliceHitsViewAsym", ";Slice Hits View Asymmetry;Slices", 200, -10., 10.);
640  fFDSliceCosAngleShw = tfs->make<TH1F>("fdSliceCosAngleShw", ";cos(#theta_{showers});Slices", 100, -1., 1.);
641  fFDSliceDistShwVtx = tfs->make<TH1F>("fdSliceDistShwVtx", ";#DeltaL_{shower,vertex};Slices", 100, 0., 400.);
642  fFDSliceFracShwHits = tfs->make<TH1F>("fdSliceFracShwHits", ";Hits_{shower}/Hits_{slice};Slices", 50, 0., 2.);
643  fFDVtxDistToEastWall = tfs->make<TH1F>("fdVtxDistToEastWall", ";Distance to East Wall;Showers", 100, 0., 100.);
644  fFDVtxDistToWestWall = tfs->make<TH1F>("fdVtxDistToWestWall", ";Distance to West Wall;Showers", 100, 0., 100.);
645  fFDVtxDistToTop = tfs->make<TH1F>("fdVtxDistToTop", ";Distance to Top;Showers", 100, 0., 400.);
646  fFDVtxDistToBottom = tfs->make<TH1F>("fdVtxDistToBottom", ";Distance to Bottom;Showers", 100, 0., 400.);
647  fFDVtxDistToFront = tfs->make<TH1F>("fdVtxDistToFront", ";Distance to Front;Showers", 100, 0., 100.);
648  fFDVtxDistToBack = tfs->make<TH1F>("fdVtxDistToBack", ";Distance to Back;Showers", 100, 0., 400.);
649  fFDPt = tfs->make<TH1F>("fdPt", ";p_{t};Showers", 50, 0., 2.);
650 
651  // now the ND histograms
652  fNDNShowers = tfs->make<TH1F>("ndNshowers", ";Nnumber of Showers;Slices", 10, 0., 10.);
653  fNDShowerCalE = tfs->make<TH1F>("ndshowercale", ";Shower E_{calorimetric};Slices", 120, 0., 6.);
654  fNDShowerCalEfraction = tfs->make<TH1F>("ndshowercalefrac", ";Fraction of Slice E in Shower;Slices", 110, 0., 1.1);
655  fNDSliceNoShowerCalE = tfs->make<TH1F>("ndnoshowercale", ";(Slice -Shower) E_{calorimetric};Slices",120, 0., 6.);
656  fNDShowerDistToVtx = tfs->make<TH1F>("ndNshowergap", ";Shower Distance to Vertex;Slices", 25, 0., 200.);
657  fNDVtxX = tfs->make<TH1F>("ndVtxX", ";x_{vertex} (cm);Showers", 200, -200., 200.);
658  fNDVtxY = tfs->make<TH1F>("ndVtxY", ";y_{vertex} (cm);Showers", 200, -200., 200.);
659  fNDVtxZ = tfs->make<TH1F>("ndVtxZ", ";z_{vertex} (cm);Showers", 700, 0., 1400.);
660  fNDShwVtxX = tfs->make<TH1F>("ndShwVtxX", ";x_{vertex} (cm);Showers", 200, -200., 200.);
661  fNDShwVtxY = tfs->make<TH1F>("ndShwVtxY", ";y_{vertex} (cm);Showers", 200, -200., 200.);
662  fNDShwVtxZ = tfs->make<TH1F>("ndShwVtxZ", ";z_{vertex} (cm);Showers", 700, 0., 1400.);
663  fNDShwEndX = tfs->make<TH1F>("ndShwEndX", ";x_{vertex} (cm);Showers", 200, -200., 200.);
664  fNDShwEndY = tfs->make<TH1F>("ndShwEndY", ";y_{vertex} (cm);Showers", 200, -200., 200.);
665  fNDShwEndZ = tfs->make<TH1F>("ndShwEndZ", ";z_{vertex} (cm);Showers", 700, 0., 1400.);
666  fNDSliceHits = tfs->make<TH1F>("ndSliceHits", ";Hits in Slice;Slices", 500, 0., 500.);
667  fNDCalorimetricE = tfs->make<TH1F>("ndCalorimetricE", ";E_{calorimetric} (GeV);Showers", 120, 0., 6.);
668  fNDCalEperHit = tfs->make<TH1F>("ndCalEperHit", ";E_{calorimetric}/ hit ;Slices", 100, 0., 0.1);
669  fNDShowerLength = tfs->make<TH1F>("ndShowerLength", ";L_{shower} (cm);Showers", 700, 0., 1400.);
670  fNDSliceNHitsPerPlane = tfs->make<TH1F>("ndSliceHitsPerPlane", ";Slice Hits/Plane;Slices", 100, 0., 50.);
671  fNDPlanesToFront = tfs->make<TH1F>("ndPlanesToFront", ";Planes to Front;Events", 600, 0., 600.);
672 
673  std::string paramst("fNumPlane/F:fPlanesToFront/F:fCellsPerPlane/F:fNumXCell/F:fNumYCell/F:");
674  paramst.append("fSliceHits/F:fHitAsymmetry/F:fNumShowers/F:fShowerNumXCell/F:fShowerNumYCell/F:");
675  paramst.append("fCosShowers/F:fShowerDirX/F:fShowerDirY/F:fShowerDirZ/F:");
676  paramst.append("fShowerCalE/F:fShowerLength/F:fFracShowerHits/F:");
677  paramst.append("fShowerVertexX/F:fShowerVertexY/F:fShowerVertexZ/F:");
678  paramst.append("fShowerEndX/F:fShowerEndY/F:fShowerEndZ/F:");
679  paramst.append("fShowerMinX/F:fShowerMinY/F:fShowerMaxX/F:fShowerMaxY/F:fShowerMinVtxZ/F:");
680  paramst.append("fShowerMinEndZ/F:fShowerMaxVtxZ/F:fShowerMaxEndZ/F:fShowerVtxDist/F:fShowerPt/F:");
681  paramst.append("fSliceMinPlane/F:fVertexX/F:fVertexY/F:fVertexZ/F:fVertexMaxW/F:fCalorimetricE/F:fNueEnergy/F:");
682  paramst.append("fMinWestDist/F:fMinEastDist/F:fMinBotDist/F:fMinTopDist/F:fMinFrontDist/F:fMinBackDist:");
683  paramst.append("fLIDVal/F:fLEMVal/F:fCVNVal/F");
684  fNueTree = tfs->make<TTree>("Nue", "Skimmed Nue Interaction Cut Parameters");
685  fNueTree->Branch("info", &fEventInfo, infost .c_str());
686  fNueTree->Branch("params", &fNueCutParams, paramst.c_str());
687  fNueTree->Branch("sr", &fNueSRParams, paramst.c_str());
688 
689  }
690 
691  // scan all the parameters in an interactive ROOT session
692  // 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");
693 
694  // compare the values from skimming vs CAF
695  // 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");
696 
697 // 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");
698 
699 
700  return;
701 }
702 
703 //------------------------------------------------------------------------------
705 {
706  fDetId = fDetService->DetId();
707 
708  return;
709 }
710 
711 //------------------------------------------------------------------------------
713 {
714  // get the CosmicExposure and POTSum out of the subrun
715  //art::Handle<sumdata::CosmicExposure> ce;
716  //art::Handle<sumdata::POTSum> pot;
717  //sr.getByLabel(fCosmicExposureLabel, ce );
718  //sr.getByLabel(fPOTSumLabel, pot);
719 
720 // if( ce.isValid() ){
721 // LOG_DEBUG("CosmicBeamComparison") << sr.run() << " " << ce->totlivetime;
722 // fCosmicExposureByRun->Fill(sr.run(), ce->totlivetime);
723 // }
724 // else LOG_DEBUG("CosmicBeamComparison") << "No cosmic exposure found";
725 // if( pot.isValid() ){
726 // LOG_DEBUG("CosmicBeamComparison") << sr.run() << " " << pot->totgoodpot;
727 // fPOTByRun->Fill(sr.run(), 1.e-16 * pot->totgoodpot);
728 // }
729 // else LOG_DEBUG("CosmicBeamComparison") << "No pot summary found";
730 
731  std::vector<art::Handle<sumdata::CosmicExposure>> ce;
732  std::vector<art::Handle<sumdata::POTSum>> pot;
733  sr.getManyByType(ce );
734  sr.getManyByType(pot);
735  if( ce.size() > 0 ){
736  LOG_DEBUG("CosmicBeamComparison") << sr.run() << " " << ce[0]->totlivetime;
737  fCosmicExposureByRun->Fill(sr.run(), ce[0]->totlivetime);
738  }
739  else LOG_DEBUG("CosmicBeamComparison") << "No cosmic exposure found";
740  if( pot.size() > 0 ){
741  LOG_DEBUG("CosmicBeamComparison") << sr.run() << " " << pot[0]->totgoodpot;
742  fPOTByRun->Fill(sr.run(), 1.e-3 * pot[0]->totpot);
743  fGoodPOTByRun->Fill(sr.run(), 1.e-3 * pot[0]->totgoodpot);
744  }
745  else LOG_DEBUG("CosmicBeamComparison") << "No pot summary found";
746 
747 
748  return;
749 }
750 
751 //------------------------------------------------------------------------------
753 {
754  fEventInfo.run = e.run();
755  fEventInfo.subrun = e.subRun();
756  fEventInfo.event = e.id().event();
757 
758  LOG_DEBUG("CosmicBeamComparison")
759  << e.run() << " " << e.subRun() << " " << e.id().event();
760 
761  if( !fNumuInstance.empty() ) this->FillNumuHistograms(e);
762  if( !fNueInstance.empty() ) this->FillNueHistograms(e);
763 
764  return;
765 }
766 
767 //------------------------------------------------------------------------------
769 {
770 
771  fSkimmedLabel = p.get<std::string >("SkimmedLabel");
772  fCosmicExposureLabel = p.get<std::string >("CosmicExposureLabel", "ifdbspillinfo");
773  fPOTSumLabel = p.get<std::string >("POTSumLabel", "exposure" );
774  fCompareNumuCuts = p.get<bool >("CompareNumuCuts", false);
775  fCompareNueCuts = p.get<bool >("CompareNueCuts", false);
776  auto instanceLabels = p.get<std::vector<std::string> >("Instances");
777 
778  for(auto itr : instanceLabels){
779  if(itr.find("numu") != std::string::npos) fNumuInstance = itr;
780  if(itr.find("nue") != std::string::npos) fNueInstance = itr;
781  }
782 
783  return;
784 }
785 
786 //------------------------------------------------------------------------------
788 {
789  // get the relevant data products out of the event record
790  // need the slice, the vertex for the slice, and the showers in the slice
792 
794  e.getByLabel(tag, sliceHandle);
795  if( sliceHandle.failedToGet() ){
796  LOG_WARNING("CosmicBeamComparison") << "No slices passing nue selection in this event";
797  return;
798  }
799 
800  // Get the standard record to compare the values used for the cuts
801  //art::FindOne<caf::StandardRecord> fosr(sliceHandle, e, tag);
802 
803  art::FindMany<rb::Shower> fms (sliceHandle, e, tag);
804  art::FindMany<rb::Prong> fmp (sliceHandle, e, tag);
805  art::FindMany<rb::Vertex> fmv (sliceHandle, e, tag);
806  art::FindMany<cosrej::NueCosRej> fmc (sliceHandle, e, tag);
807  art::FindMany<slid::EventLID> fmlid (sliceHandle, e, tag);
808  art::FindMany<lem::PIDDetails> fmlem (sliceHandle, e, tag);
809  art::FindMany<cvn::Result> fmcvn (sliceHandle, e, tag);
810  art::FindMany<cvn::Result> fmocvn (sliceHandle, e, tag);
811  art::FindMany<presel::Veto> fmnueveto(sliceHandle, e, tag);
812 
813 
815 
816  for(size_t s = 0; s < sliceHandle->size(); ++s){
817 
818  std::vector<rb::Shower const*> showers;
819  std::vector<rb::Prong const*> prongs;
820  std::vector<rb::Vertex const*> vertices;
821  std::vector<slid::ShowerLID const*> slids;
822  std::vector<slid::EventLID const*> elids;
823  std::vector<lem::PIDDetails const*> lems;
824  std::vector<cvn::Result const*> cvns;
825  std::vector<presel::Veto const*> nuevetos;
826  std::vector<cosrej::NueCosRej const*> cosrejs;
827 
828  su.FillSliceVector(s, fmv, vertices);
829  su.FillSliceVector(s, fmc, cosrejs);
830  su.FillSliceVector(s, fms, showers);
831  su.FillSliceVector(s, fmp, prongs);
832  su.FillSliceVector(s, fmlid, elids);
833  su.FillSliceVector(s, fmlem, lems);
834  su.FillSliceVector(s, fmcvn, cvns);
835  su.FillSliceVector(s, fmnueveto, nuevetos);
836  su.FillSliceVector(e, tag, showers, slids);
837 
838  if(vertices.size() != 1 ||
839  cosrejs .size() != 1 ||
840  lems .size() != 1 ||
841  cvns .size() != 1 ||
842  elids .size() != 1 ||
843  showers .size() < 1 ||
844  prongs .size() < 1 ||
845  slids .size() != showers.size()){
846  LOG_DEBUG("CosmicBeamComparison")
847  << "Wrong number of vertices: " << vertices.size()
848  << " : 1 or cosrejs: " << cosrejs .size()
849  << " : 1 or lids: " << elids .size()
850  << " : 1 or lems: " << lems .size()
851  << " : 1 or cvns: " << cvns .size()
852  << " : 1 or showers: " << showers .size()
853  << " : or prongs: " << prongs .size()
854  << " or slids: " << slids .size()
855  << " for the skimmed slice, that should never happen";
856  continue;
857  }
858 
859  // make the parameters object
860  skim::ParametersNue params((*sliceHandle)[s],
861  *(vertices[0]),
862  *(cosrejs[0]),
863  showers,
864  prongs,
865  slids,
866  *(elids[0]),
867  *(lems[0]),
868  *(cvns[0]),
869  *(nuevetos[0]));
870 
871 
872  //skim::ParametersNue sr(fosr.at(s).ref());
873 
874  fNueCutParams = params.ParametersStruct();
875  //fNueSRParams = sr .ParametersStruct();
876 
877  //if (fCompareNueCuts) params.Compare(fNueSRParams);
878 
879  fLEM ->Fill(params.LEMVal());
880  fCVN ->Fill(params.CVNVal());
881  fLID ->Fill(params.LIDVal());
882  fCalorimetricEnergy->Fill(params.CalorimetricE());
883 
885  fFDSliceHitsPerPlane->Fill(params.CellsPerPlane());
886  fFDSliceHitsXView ->Fill(params.ShowerNumXCell());
887  fFDSliceHitsYView ->Fill(params.ShowerNumYCell());
888  fFDSliceHitsViewAsym->Fill(params.HitAsymmetry());
889  fFDSliceCosAngleShw ->Fill(params.CosShowers());
890  fFDSliceDistShwVtx ->Fill(params.ShowerVtxDist());
891  fFDSliceFracShwHits ->Fill(params.FracShowerHits());
892  fFDVtxDistToEastWall->Fill(params.MinEastDist());
893  fFDVtxDistToWestWall->Fill(params.MinWestDist());
894  fFDVtxDistToTop ->Fill(params.MinTopDist());
895  fFDVtxDistToBottom ->Fill(params.MinBotDist());
896  fFDVtxDistToFront ->Fill(params.MinFrontDist());
897  fFDVtxDistToBack ->Fill(params.MinBackDist());
898  fFDPt ->Fill(params.ShowerPt());
899  }
900 
902  fNDVtxX ->Fill(params.Vertex().X() );
903  fNDVtxY ->Fill(params.Vertex().Y() );
904  fNDVtxZ ->Fill(params.Vertex().Z() );
905  fNDShwVtxX ->Fill(params.ShowerVertex().X());
906  fNDShwVtxY ->Fill(params.ShowerVertex().Y());
907  fNDShwVtxZ ->Fill(params.ShowerVertex().Z());
908  fNDShwEndX ->Fill(params.ShowerEnd().X() );
909  fNDShwEndY ->Fill(params.ShowerEnd().Y() );
910  fNDShwEndZ ->Fill(params.ShowerEnd().Z() );
911  fNDSliceHits ->Fill(params.SliceHits() );
912  fNDCalorimetricE ->Fill(params.CalorimetricE() );
913  fNDShowerLength ->Fill(params.ProngLength() );
914  fNDSliceNHitsPerPlane->Fill(params.CellsPerPlane() );
915  fNDPlanesToFront ->Fill(params.PlanesToFront() );
916  fNDNShowers ->Fill(params.NumShowers() );
917  fNDShowerDistToVtx ->Fill(params.ShowerVtxDist() );
918  fNDShowerCalE ->Fill(params.ShowerCalE() );
919  fNDShowerCalEfraction->Fill(params.ShowerCalE() / params.CalorimetricE());
920  fNDSliceNoShowerCalE ->Fill(params.CalorimetricE() - params.ShowerCalE() );
921  fNDCalEperHit ->Fill(params.CalorimetricE() / params.SliceHits() );
922  }
923 
924  fNueTree->Fill();
925 
926  } // end loop over skimmed slices
927 
928  return;
929 }
930 
931 //Stopgap
932 std::pair<int, int> calcFirstLastLivePlane(int plane,
933  std::bitset<14> binary,
935  //if not the FD, return full ND size
936  if (det != novadaq::cnv::kFARDET) return std::pair<int, int>(0,213);
937  int testDB = (plane / 64);
938  int minblock = testDB;
939  int maxblock = testDB;
940  //find block boundaries;
941  for (int i=testDB-1; i>=0; --i){
942  if (binary[i]) minblock=i;
943  else break;
944  }
945  for (int i=testDB+1; i<14; ++i){
946  if (binary[i]) maxblock=i;
947  else break;
948  }
949  return std::pair<int, int>(64*(minblock), 64*(maxblock+1)-1);
950 }
951 
952 //------------------------------------------------------------------------------
954 {
956  art::InputTag cosinputTag(fSkimmedLabel, fNumuInstance+"Cosmics");
957 
958  // get the relevant data products out of the event record
959  // need the slice, rb::Energy, remid::ReMId, cosrej::CosRejObj,
960  // qeef::QePId, numue::NumuE, rb::Track
961  auto sliceHandle = e.getValidHandle< std::vector<rb::Cluster> >(inputTag);
962  if( sliceHandle.failedToGet() ){
963  LOG_WARNING("CosmicBeamComparison") << "No slices passing numu selection in this event";
964  return;
965  }
966 
967  // Get the standard record to compare the values used for the cuts
968  //art::FindOne<caf::StandardRecord> fosr (sliceHandle, e, inputTag);
969 
970  art::FindManyP<rb::Track> fmt (sliceHandle, e, inputTag);
971  art::FindManyP<rb::Track> fmct (sliceHandle, e, cosinputTag);
972  art::FindMany<rb::Energy> fme (sliceHandle, e, inputTag);
973  art::FindMany<cosrej::CosRejObj> fmcro (sliceHandle, e, inputTag);
974  art::FindMany<qeef::QePId> fmqe (sliceHandle, e, inputTag);
975  art::FindMany<numue::NumuE> fmnue (sliceHandle, e, inputTag);
976  art::FindMany<cvn::Result> fmcvn (sliceHandle, e, inputTag);
977  art::FindMany<cvn::Result> fmocvn(sliceHandle, e, inputTag);
978 
979  art::FindMany<cosrej::NueCosRej> fmnuecro (sliceHandle, e, inputTag); // Nue CosRej added for Analysis 2017
980  art::FindMany<rb::Track> fmtks (sliceHandle, e, inputTag);
981  art::FindMany<rb::Shower> fms (sliceHandle, e, inputTag);
982  art::FindMany<rb::Prong> fmp (sliceHandle, e, inputTag);
983  art::FindMany<rb::Vertex> fmv (sliceHandle, e, inputTag);
984  art::FindMany<slid::EventLID> fmlid (sliceHandle, e, inputTag);
985 
986  /*
987  std::vector< art::Ptr<rb::Track> > trackPtrs;
988  std::vector< art::Ptr<rb::Track> > costrackPtrs;
989  std::vector<cosrej::CosRejObj const*> cosrej;
990  std::vector<qeef::QePId const*> qepid;
991  std::vector<numue::NumuE const*> numue;
992  */
993  /*
994  // std::vector<cosrej::NueCosRej const*> nuecosrej; // Nue CosRej added for Analysis 2017
995  std::vector<cvn::Result const*> cvns; /// 2017
996  std::vector<rb::Shower const*> showers;
997  std::vector<rb::Prong const*> prongs;
998  std::vector<rb::Vertex const*> vertices;
999  std::vector<slid::ShowerLID const*> slids;
1000  */
1001  //==> size_t bestPIDTrack = std::numeric_limits<size_t>::max();
1002 
1004 
1005  for(size_t s = 0; s < sliceHandle->size(); ++s){
1006 
1007 
1008  std::vector< art::Ptr<rb::Track> > trackPtrs;
1009  std::vector< art::Ptr<rb::Track> > costrackPtrs;
1010  std::vector<cosrej::CosRejObj const*> cosrej;
1011  std::vector<qeef::QePId const*> qepid;
1012  std::vector<numue::NumuE const*> numue;
1013 
1014  std::vector<cosrej::NueCosRej const*> nuecosrej; // Nue CosRej added for Analysis 2017
1015  std::vector<cvn::Result const*> cvns; /// 2017
1016  std::vector<cvn::Result const*> cvn2017s; /// 2018
1017  std::vector<rb::Track const*> tracks;
1018  std::vector<rb::Shower const*> showers;
1019  std::vector<rb::Prong const*> prongs;
1020  std::vector<rb::Vertex const*> vertices;
1021  std::vector<slid::ShowerLID const*> slids;
1022 
1023  size_t bestPIDTrack = std::numeric_limits<size_t>::max();
1024 
1025 // if(fosr.at(s).ref().hdr.ismc)//only fill truth information if it's there
1026 // {
1027 // crTruthNNu = fosr.at(s).ref().mc.allnus.size();
1028 // if(fosr.at(s).ref().mc.nu.size() > 0)
1029 // {
1030 // crTruthNuNHitsSlc = fosr.at(s).ref().mc.nu[0].nhitslc;
1031 // crTruthNuNHitsTot = fosr.at(s).ref().mc.nu[0].nhittot;
1032 // crTruthNuIsCC = fosr.at(s).ref().mc.nu[0].iscc;
1033 // crTruthNuPDG = fosr.at(s).ref().mc.nu[0].pdg;
1034 // crTruthNuPDGOrig = fosr.at(s).ref().mc.nu[0].pdgorig;
1035 // crTruthNuVertex = fosr.at(s).ref().mc.nu[0].vtx;
1036 // crTruthNuTime = fosr.at(s).ref().mc.nu[0].time;
1037 // }
1038 // }
1039 
1040 
1041 // fmt.get(s, trackPtrs);
1042 // fmct.get(s, costrackPtrs);
1043 
1044  // better have at least one track
1045  if(trackPtrs.size() < 1) continue;
1046 
1047  // get the best pid for the slice - this index will not be the
1048  // same as in the preskimmed slice
1049  bestPIDTrack = remid::HighestPIDTrack(trackPtrs, inputTag, e);
1050 
1051  fmcro.get(s, cosrej);
1052  fmnue.get(s, numue);
1053  fmnuecro.get(s, nuecosrej); // Nue CosRej added for Analysis 2017
1054  fmcvn.get( s, cvns ); // 2017
1055  fmocvn.get(s, cvn2017s); // 2018
1056 
1057  su.FillSliceVector(s, fmtks, tracks);
1058  su.FillSliceVector(s, fmv, vertices);
1059  su.FillSliceVector(s, fms, showers);
1060  su.FillSliceVector(s, fmp, prongs);
1061  su.FillSliceVector(e, inputTag, showers, slids);
1062 
1063 
1064  //fmqe .get(s, qepid);
1065  //qeef::QePId qep;
1066  //if( qepid.size() > 0 ) qep = *(qepid.front());
1067 
1068 
1069  if(vertices.size() != 1 ||
1070  cvns .size() != 1 ||
1071  cvn2017s.size() != 1 ||
1072  showers .size() < 1 ||
1073  prongs .size() < 1 ||
1074  tracks .size() < 1 ||
1075  slids .size() != showers.size()){
1076  LOG_DEBUG("CosmicBeamComparison")
1077  << "Wrong number of vertices: " << vertices.size()
1078  << " : 1 or cvns: " << cvns .size()
1079  << " : 1 or cvn2017s: " << cvn2017s.size()
1080  << " : 1 or showers: " << showers .size()
1081  << " : or prongs: " << prongs .size()
1082  << " : or tracks: " << tracks .size()
1083  << " or slids: " << slids .size()
1084  << " for the skimmed slice, that should never happen";
1085  continue;
1086  }
1087 
1088  // this is not terribly efficient, as we are making the FindOneP
1089  // for every slice in the event, but I don't see a good way around
1090  // that while still getting the bookkeeping right
1091  art::FindOne<rb::Energy> foe (trackPtrs, e, inputTag);
1092  art::FindOne<remid::ReMId> foid(trackPtrs, e, inputTag);
1093 
1094  //LOG_VERBATIM("CosmicBeamComparison") << "Filling params";
1095  skim::ParametersNumu params(bestPIDTrack,
1096  foe.at(bestPIDTrack).ref(),
1097  foid.at(bestPIDTrack).ref(),
1098  *(cosrej.front()),
1099  *(nuecosrej.front()),
1100  //qep,
1101  *(vertices[0]),
1102  tracks,
1103  showers,
1104  prongs,
1105  slids,
1106  *(costrackPtrs[0]), //Not sure what else to do here...
1107  *(numue.front()),
1108  sliceHandle->at(s),
1109  *(trackPtrs[bestPIDTrack]),
1110  *(cvns.front()),
1111  *(cvn2017s.front()),
1112  e.isRealData());
1113 
1114 // skim::ParametersNumu sr(fosr.at(s).ref(),
1115 // e.isRealData());
1116 
1117  fNumuCutParams = params.ParametersStruct();
1118  //LOG_VERBATIM("CosmicBeamComparison") << "fNumuCutParams.fReMIdVal = " << fNumuCutParams.fReMIdVal;
1119 // fNumuSRParams = sr .ParametersStruct();
1120 //
1121 // if (fCompareNumuCuts) params.Compare(fNumuSRParams);
1122 
1123  // Fill sanity check histograms
1124  fNeutrinoE ->Fill(params.NeutrinoE());
1125  fReMId ->Fill(params.ReMIdValue());
1126  fSliceHits ->Fill(params.SliceHits());
1127  fSliceContigPlanes->Fill(params.SliceContigPlanes());
1128  fCellsFromEdge ->Fill(params.SliceCellsFromEdge());
1129  if(params.QePIDNTracks() == 1) fQePId1Track->Fill(params.QePIDValue());
1130  if(params.QePIDNTracks() == 2) fQePId2Track->Fill(params.QePIDValue());
1131 
1133  fFDKalmanFwdCell->Fill(params.CosRejKalFwdCell());
1134  fFDKalmanBakCell->Fill(params.CosRejKalBakCell());
1135  fFDKalmanAngle ->Fill(params.CosRejAngleKal());
1136  fFDCosBakCell ->Fill(params.CosRejCosBakCell());
1137  fFDCosFwdCell ->Fill(params.CosRejCosFwdCell());
1138  fFDContPID ->Fill(params.CosRejConCosPID());
1139  fFDPlanesToBack ->Fill(params.PlanesToBack());
1140  fFDPlanesToFront->Fill(params.PlanesToFront());
1141  fFDSliceHits ->Fill(params.SliceHits());
1142  fFDKalmanDirX ->Fill(params.TrackDir().X());
1143  fFDKalmanDirY ->Fill(params.TrackDir().Y());
1144  fFDKalmanDirZ ->Fill(params.TrackDir().Z());
1145  fFDKalSpeed ->Fill(params.CosRejKalSpeed());
1146  fFDCosSpeed ->Fill(params.CosRejCosSpeed());
1147  fFDDeltaT ->Fill(params.TrackStopTime()-params.TrackStartTime());
1148  fFDDeltaX ->Fill(params.TrackStop().X()-params.TrackStart().X());
1149  fFDDeltaY ->Fill(params.TrackStop().Y()-params.TrackStart().Y());
1150  fFDDeltaZ ->Fill(params.TrackStop().Z()-params.TrackStart().Z());
1151  fFDDTvsDY ->Fill(params.TrackStopTime()-params.TrackStartTime(), params.TrackStop().Y()-params.TrackStart().Y());
1152  fFDCosSpeedVsScore->Fill(params.CosRejCosSpeed(),params.CosRejCosDirScore());
1153  fFDSpeed ->Fill(((params.TrackStop()-params.TrackStart()).Mag() * 1e-2) / ((params.TrackStopTime()-params.TrackStartTime()) * 1e-9));//Convert cm/ns to m/s
1154  fFDTrackLength ->Fill((params.TrackStop()-params.TrackStart()).Mag());
1155 
1156  //Fill CosmicTree - mostly the same as params, but ReMId works
1157  crNeutrinoE = params.NeutrinoE();
1158  crReMIdVal = params.ReMIdValue();
1159  crSliceHits = sliceHandle->at(s).NCell();
1160  crSliceContigPlanes = sliceHandle->at(s).MostContiguousPlanes(geo::kXorY);
1161  crSliceMaxPlane = sliceHandle->at(s).MaxPlane();
1162  crSliceMinPlane = sliceHandle->at(s).MinPlane();
1163  crSliceXViewHits = sliceHandle->at(s).NXCell();
1164  crSliceYViewHits = sliceHandle->at(s).NYCell();
1165  crSliceTotalADC = sliceHandle->at(s).TotalADC();
1166  crSliceTotalGeV = sliceHandle->at(s).TotalGeV();
1167  crSliceCalorimetricEnergy = sliceHandle->at(s).CalorimetricEnergy();
1168  crSliceMinXYZ = sliceHandle->at(s).MinXYZ();
1169  crSliceMaxXYZ = sliceHandle->at(s).MaxXYZ();
1170  crSliceMinTNS = sliceHandle->at(s).MinTNS();
1171  crSliceMaxTNS = sliceHandle->at(s).MaxTNS();
1172  crCellsFromEdge = params.SliceCellsFromEdge();
1173  crSliceClosestEdge = params.SliceClosestEdge();
1174  crSliceNHitsNearEdge0 = params.SliceNHitsNearEdge(0);
1175  crSliceNHitsNearEdge1 = params.SliceNHitsNearEdge(1);
1176  crSliceNHitsNearEdge2 = params.SliceNHitsNearEdge(2);
1177  crSliceNHitsNearEdge3 = params.SliceNHitsNearEdge(3);
1178  crSliceNHitsNearEdge4 = params.SliceNHitsNearEdge(4);
1179  crSliceNHitsNearEdge5 = params.SliceNHitsNearEdge(5);
1180  crSliceNHitsNearEdge6 = params.SliceNHitsNearEdge(6);
1181  crSliceNHitsNearEdge7 = params.SliceNHitsNearEdge(7);
1182  crSliceNHitsNearEdge8 = params.SliceNHitsNearEdge(8);
1183  crSliceNHitsNearEdge9 = params.SliceNHitsNearEdge(9);
1184  //crPlanesToFront = sr.PlanesToFront();//params.PlanesToFront()
1185  //crPlanesToBack = sr.PlanesToBack();//params.PlanesToBack() returns a lot of negative numbers, I'm not entirely convinced by it.
1186  crKalAngle = params.CosRejAngleKal();
1187  crCosAngle = cosrej.front()->AngleCos();
1188  crKalFwdCell = params.CosRejKalFwdCell();
1189  crKalBakCell = params.CosRejKalBakCell();
1190  crCosFwdCell = params.CosRejCosFwdCell();
1191  crCosBakCell = params.CosRejCosBakCell();
1192  crKalFwdDist = cosrej.front()->KalFwdDist();
1193  crKalBakDist = cosrej.front()->KalBakDist();
1194  crCosFwdDist = cosrej.front()->CosFwdDist();
1195  crCosBakDist = cosrej.front()->CosBakDist();
1196  crConCosPID = params.CosRejConCosPID();
1197  crKalDir = params.TrackDir();
1198  crTrackStart = params.TrackStart();
1199  crTrackStop = params.TrackStop();
1200  crTrackStartTime = params.TrackStartTime();
1201  crTrackStopTime = params.TrackStopTime();
1202  crTrackHits = (trackPtrs[bestPIDTrack])->NCell();
1203  crTrackCellsFromEdge = params.TrackNCellsFromEdge();
1204  crTrackClosestEdge = params.TrackClosestEdge();
1205  crTrackNHitsNearEdge0 = params.TrackNHitsNearEdge(0);
1206  crTrackNHitsNearEdge1 = params.TrackNHitsNearEdge(1);
1207  crTrackNHitsNearEdge2 = params.TrackNHitsNearEdge(2);
1208  crTrackNHitsNearEdge3 = params.TrackNHitsNearEdge(3);
1209  crTrackNHitsNearEdge4 = params.TrackNHitsNearEdge(4);
1210  crTrackNHitsNearEdge5 = params.TrackNHitsNearEdge(5);
1211  crTrackNHitsNearEdge6 = params.TrackNHitsNearEdge(6);
1212  crTrackNHitsNearEdge7 = params.TrackNHitsNearEdge(7);
1213  crTrackNHitsNearEdge8 = params.TrackNHitsNearEdge(8);
1214  crTrackNHitsNearEdge9 = params.TrackNHitsNearEdge(9);
1215  crTrackXViewHits = (trackPtrs[bestPIDTrack])->NXCell();
1216  crTrackYViewHits = (trackPtrs[bestPIDTrack])->NYCell();
1217  crTrackEarliestHitPos = params.TrackEarliestHitPos();
1218  crTrackLatestHitPos = params.TrackLatestHitPos();
1219  crTrackTotalADC = (trackPtrs[bestPIDTrack])->TotalADC();//sum of adc of all hits
1220  crTrackTotalGeV = (trackPtrs[bestPIDTrack])->TotalGeV();//simple sum of estimated gev of all hits
1221  crTrackCalorimetricEnergy = (trackPtrs[bestPIDTrack])->CalorimetricEnergy();//simple estimate of neutrino energy
1222  crTrackMinTNS = (trackPtrs[bestPIDTrack])->MinTNS();
1223  crTrackMaxTNS = (trackPtrs[bestPIDTrack])->MaxTNS();
1224  crTrackMinXYZ = (trackPtrs[bestPIDTrack])->MinXYZ();
1225  crTrackMaxXYZ = (trackPtrs[bestPIDTrack])->MaxXYZ();
1226  crTrackMinPlane = (trackPtrs[bestPIDTrack])->MinPlane();
1227  crTrackMaxPlane = (trackPtrs[bestPIDTrack])->MaxPlane();
1228  crTrackTotalLength = (trackPtrs[bestPIDTrack])->TotalLength();
1229  crCosTrackStart = params.CosmicTrackStart();
1230  crCosTrackStop = params.CosmicTrackStop();
1231  crCosTrackStartTime = params.CosTrackStartTime();
1232  crCosTrackStopTime = params.CosTrackStopTime();
1233  crCosTrackDir = params.CosmicTrackDir();
1234  crCosTrackHits = params.CosTrackHits();
1235  crCosTrackCellsFromEdge = params.CosTrackNCellsFromEdge();
1236  crCosTrackClosestEdge = params.CosTrackClosestEdge();
1237  crCosTrackNHitsNearEdge0 = params.CosTrackNHitsNearEdge(0);
1238  crCosTrackNHitsNearEdge1 = params.CosTrackNHitsNearEdge(1);
1239  crCosTrackNHitsNearEdge2 = params.CosTrackNHitsNearEdge(2);
1240  crCosTrackNHitsNearEdge3 = params.CosTrackNHitsNearEdge(3);
1241  crCosTrackNHitsNearEdge4 = params.CosTrackNHitsNearEdge(4);
1242  crCosTrackNHitsNearEdge5 = params.CosTrackNHitsNearEdge(5);
1243  crCosTrackNHitsNearEdge6 = params.CosTrackNHitsNearEdge(6);
1244  crCosTrackNHitsNearEdge7 = params.CosTrackNHitsNearEdge(7);
1245  crCosTrackNHitsNearEdge8 = params.CosTrackNHitsNearEdge(8);
1246  crCosTrackNHitsNearEdge9 = params.CosTrackNHitsNearEdge(9);
1247  crCosTrackXViewHits = params.CosTrackXViewHits();
1248  crCosTrackYViewHits = params.CosTrackXViewHits();
1249  crCosTrackTotalADC = params.CosTrackTotalADC();
1250  crCosTrackTotalGeV = params.CosTrackTotalGeV();
1251  crCosTrackCalorimetricEnergy = params.CosTrackTotalCalorimetricEnergy();
1252  crCosTrackMinTNS = params.CosTrackMinTNS();
1253  crCosTrackMaxTNS = params.CosTrackMaxTNS();
1254  crCosTrackMinXYZ = params.CosTrackMinXYZ();
1255  crCosTrackMaxXYZ = params.CosTrackMaxXYZ();
1256  crCosTrackEarliestHitPos = params.CosmicTrackEarliestHitPos();
1257  crCosTrackLatestHitPos = params.CosmicTrackLatestHitPos();
1258  crCosTrackMinPlane = params.CosTrackMinPlane();
1259  crCosTrackMaxPlane = params.CosTrackMaxPlane();
1260  crCosTrackTotalLength = (costrackPtrs[0])->TotalLength();
1261  crKalSpeed = params.CosRejKalSpeed();
1262  crCosSpeed = params.CosRejCosSpeed();
1263  crCosDirScore = params.CosRejCosDirScore();
1264  crKalDirScore = cosrej.front()->KDirScore();
1265  crCosScoreDiff = cosrej.front()->CScoreDiff();
1266  crKalScoreDiff = cosrej.front()->KScoreDiff();
1267  crKalSlope = cosrej.front()->KalSlope();
1268  crCosSlope = cosrej.front()->CosSlope();
1269  crKalChiSq = cosrej.front()->KalChisq();
1270  crCosChiSq = cosrej.front()->CosChisq();
1271  crKalChiDiff = cosrej.front()->KalChiDiff();
1272  crCosChiDiff = cosrej.front()->CosChiDiff();
1273  crNTracks3D = cosrej.front()->NTracks3D();
1274  crScatt = cosrej.front()->Scatt();
1275  //crQePIDVal = qep.Value();
1276  //crQePIDNTracks = qep.Ntrk();
1277 // crXmin2 = fosr.at(s).ref().sel.contain.xmin2;
1278 // crYmin2 = fosr.at(s).ref().sel.contain.ymin2;
1279 // crZmin2 = fosr.at(s).ref().sel.contain.zmin2;
1280 // crXmax2 = fosr.at(s).ref().sel.contain.xmax2;
1281 // crYmax2 = fosr.at(s).ref().sel.contain.ymax2;
1282 // crZmax2 = fosr.at(s).ref().sel.contain.zmax2;
1283 // crscatllh = fosr.at(s).ref().sel.remid.scatllh;
1284 // crdedxllh = fosr.at(s).ref().sel.remid.dedxllh;
1285  crmeasfrac = foid.at(bestPIDTrack).ref().MeasurementFraction();
1286 
1288  unsigned int dibmask = (rh->GoodDiBlockMask()&rh->GetConfiguration());
1289  std::bitset<14> binary(dibmask);
1290  std::pair<int,int> planesA = calcFirstLastLivePlane(params.SliceMinPlane(), binary, fDetId);
1291  std::pair<int,int> planesB = calcFirstLastLivePlane(params.SliceMaxPlane(), binary, fDetId);
1292  if ((planesA.first != planesB.first) || (planesA.second != planesB.second)) crNCloseLiveDiblocks = -1000;
1293  else crNCloseLiveDiblocks = (planesA.second - planesA.first + 1) * 0.015625;
1294 
1295  crRun = e.run();
1296  crSubrun = e.subRun();
1297  crEvent = e.id().event();
1298 
1299  fCosmicTree->Fill();
1300  }
1301 
1302  LOG_VERBATIM("CosmicBeamComparison") << "About to fill hittree";
1303  for(unsigned int hitnumber = 0; hitnumber < (trackPtrs[bestPIDTrack])->NCell(); hitnumber++){
1304  rb::RecoHit calibhit = (trackPtrs[bestPIDTrack])->RecoHit((trackPtrs[bestPIDTrack])->Cell(hitnumber));//Calibrated hit
1305  htIndex = hitnumber;
1306  if(htIndex==0) htTrackIndex++;
1307  htUncalibTNS = (trackPtrs[bestPIDTrack])->Cell(hitnumber)->TNS();
1308  htPlane = (trackPtrs[bestPIDTrack])->Cell(hitnumber)->Plane();
1309  htCell = (trackPtrs[bestPIDTrack])->Cell(hitnumber)->Cell();
1310  geo::View_t* view = 0;
1311  double hitpos[3];
1312  double dpos[3];//half dimensions of cell in xyz - need it for the cellinfo function
1313  //const geo::OfflineChan offline_channel(htPlane, htCell);
1314  fGeom->CellInfo(htPlane, htCell, view, hitpos, dpos);
1315  //offline_channel.GetCenter
1316  htUncalibXYZ.SetXYZ(hitpos[0], hitpos[1], hitpos[2]);
1317 
1318  htCalibTNS = calibhit.T();
1319  htCalibXYZ.SetX(calibhit.X());
1320  htCalibXYZ.SetY(calibhit.Y());
1321  htCalibXYZ.SetZ(calibhit.Z());
1322  htCalibPE = calibhit.PECorr();
1323 
1324  //*view = (trackPtrs[bestPIDTrack])->Cell(hitnumber)->View();
1325  //if(!view) LOG_VERBATIM("CosmicBeamComparison") << "View not assigned";
1326  //if(htPlane %2 == 0) htXYZ.SetY((trackPtrs[bestPIDTrack])->W((trackPtrs[bestPIDTrack]->Cell(hitnumber)).get()));//was if(*view == geo::kX)
1327  //else htXYZ.SetX((trackPtrs[bestPIDTrack])->W((trackPtrs[bestPIDTrack]->Cell(hitnumber)).get()));
1328  htUncalibW = (trackPtrs[bestPIDTrack])->W((trackPtrs[bestPIDTrack]->Cell(hitnumber)).get());
1329  htDistFromStart = (trackPtrs[bestPIDTrack])->DistanceFromStart(htCalibXYZ.Z());
1330  htDistFromStop = (trackPtrs[bestPIDTrack])->DistanceFromEnd(htCalibXYZ.Z());
1331  htDir = (trackPtrs[bestPIDTrack])->InterpolateDir(htCalibXYZ.Z());
1332  double closestx;
1333  double closesty;
1334  (trackPtrs[bestPIDTrack])->InterpolateXY(htCalibXYZ.Z(), closestx, closesty);
1335  htTrackXYZ.SetXYZ(closestx, closesty, htCalibXYZ.Z());
1336 
1337  fHitTree->Fill();
1338  }
1339 
1340  LOG_VERBATIM("CosmicBeamComparison") << "About to fill cosmic hittree";
1341  for(unsigned int hitnumber = 0; hitnumber < (costrackPtrs[0])->NCell(); hitnumber++){
1342  chtIndex = hitnumber;
1343  chtUncalibTNS = (costrackPtrs[0])->Cell(hitnumber)->TNS();
1344  chtPlane = (costrackPtrs[0])->Cell(hitnumber)->Plane();
1345  chtCell = (costrackPtrs[0])->Cell(hitnumber)->Cell();
1346  geo::View_t* view = 0;
1347  double hitpos[3];
1348  double dpos[3];//half dimensions of cell in xyz - need it for the cellinfo function
1349  //const geo::OfflineChan offline_channel(chtPlane, chtCell);
1350  fGeom->CellInfo(chtPlane, chtCell, view, hitpos, dpos);
1351  //offline_channel.GetCenter
1352  chtUncalibXYZ.SetXYZ(hitpos[0], hitpos[1], hitpos[2]);
1353  //*view = (costrackPtrs[0])->Cell(hitnumber)->View();
1354  //if(!view) LOG_VERBATIM("CosmicBeamComparison") << "View not assigned";
1355  //if(chtPlane %2 == 0) chtUncalibXYZ.SetY((costrackPtrs[0])->W((costrackPtrs[0]->Cell(hitnumber)).get()));//was if(*view == geo::kX)
1356  //else chtUncalibXYZ.SetX((costrackPtrs[0])->W((costrackPtrs[0]->Cell(hitnumber)).get()));
1357  chtUncalibW = (costrackPtrs[0])->W((costrackPtrs[0]->Cell(hitnumber)).get());
1358  chtDistFromStart = (costrackPtrs[0])->DistanceFromStart(chtUncalibXYZ.Z());
1359  chtDistFromStop = (costrackPtrs[0])->DistanceFromEnd(chtUncalibXYZ.Z());
1360  chtDir = (costrackPtrs[0])->InterpolateDir(chtUncalibXYZ.Z());
1361  double closestx;
1362  double closesty;
1363  (costrackPtrs[0])->InterpolateXY(chtUncalibXYZ.Z(), closestx, closesty);
1364  chtTrackXYZ.SetXYZ(closestx, closesty, chtUncalibXYZ.Z());
1365 
1366  fCosmicHitTree->Fill();
1367  }
1368 
1370  fNDMaxPlaneHit ->Fill(params.SliceMaxPlane());
1371  fNDMinPlaneHit ->Fill(params.SliceMinPlane());
1372  fNDTrackStartZ ->Fill(params.TrackStart().Z());
1373  fNDTrackStopZ ->Fill(params.TrackStop().Z());
1374  fNDKalmanBakCell ->Fill(params.CosRejKalBakCellND());
1375  fNDKalmanFwdCell ->Fill(params.CosRejKalFwdCellND());
1376  fNDHadCalEnergySum->Fill(params.NDHadronicCal());
1377  }
1378 
1379  fNumuTree->Fill();
1380 
1381  }// end loop over slices
1382 
1383  return;
1384 }
1385 
std::string fNumuInstance
label for numu products
float T() const
Definition: RecoHit.h:39
void beginRun(art::Run const &r) override
std::string fNueInstance
label for nue products
void FillNueHistograms(art::Event const &e)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
float const & PlanesToFront() const
Definition: ParametersNue.h:70
bool fCompareNueCuts
bool for comparing the nue cuts
novadaq::cnv::DetId DetId() const
What detector are we in?
float const & CalorimetricE() const
float const & ShowerNumYCell() const
Definition: ParametersNue.h:78
SubRunNumber_t subRun() const
Definition: Event.h:72
Energy estimators for CC events.
Definition: FillEnergies.h:7
TH2 * rh
Definition: drawXsec.C:5
float const & ShowerPt() const
Definition: ParametersNue.h:88
float const & CVNVal() const
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
float const & ProngLength() const
Definition: ParametersNue.h:80
TVector3 const & ShowerEnd() const
Definition: ParametersNue.h:85
std::string fCosmicExposureLabel
label of module creating cosmic exposure
const char * p
Definition: xmltok.h:285
art::ServiceHandle< ds::DetectorService > fDetService
detector service
size_type get(size_type i, reference item, data_reference data) const
Definition: FindMany.h:469
float Z() const
Definition: RecoHit.h:38
std::pair< int, int > calcFirstLastLivePlane(int plane, std::bitset< 14 > binary, novadaq::cnv::DetId det)
TFile * fmc
Definition: bdt_com.C:14
CosmicBeamComparison(fhicl::ParameterSet const &p)
float const & MinEastDist() const
bool isRealData() const
Definition: Event.h:83
DEFINE_ART_MODULE(TestTMapFile)
NumuCutParameters fNumuSRParams
from the ParametersNumu object
Module to create a summary of total POT seen in a job.
Definition: Evaluator.h:27
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
float const & CellsPerPlane() const
Definition: ParametersNue.h:71
Definition: Run.h:31
void FillSliceVector(size_t const &s, art::FindMany< T > &fmp, std::vector< T const * > &ptrMap)
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
TVector3 const & Vertex() const
Definition: ParametersNue.h:98
std::string fPOTSumLabel
label of module creating pot sum
std::string fSkimmedLabel
label of module creating skimmed file
float const & ShowerCalE() const
Definition: ParametersNue.h:81
const XML_Char * s
Definition: expat.h:262
TTree * fCosmicTree
Tree filled directly from cosrej etc. objects (bypasses parameters struct)
void FillNumuHistograms(art::Event const &e)
Far Detector at Ash River, MN.
void beginSubRun(art::SubRun const &sr) override
Result for CVN.
void analyze(art::Event const &e) override
bool fCompareNumuCuts
bool for comparing the numu cuts
float const & MinTopDist() const
novadaq::cnv::DetId fDetId
the id of the detector we are looking at
float const & HitAsymmetry() const
Definition: ParametersNue.h:75
art::ServiceHandle< geo::Geometry > fGeom
T get(std::string const &key) const
Definition: ParameterSet.h:231
TVector3 const & ShowerVertex() const
Definition: ParametersNue.h:86
TTree * fNumuTree
Tree to hold values of cut parameters for selected events.
NueCutParameters fNueCutParams
from the ParametersNumu object
#define pot
NueCutParameters const & ParametersStruct() const
Near Detector in the NuMI cavern.
void getManyByType(std::vector< Handle< PROD >> &results) const
Definition: DataViewImpl.h:446
caf::StandardRecord * sr
float const & MinBackDist() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
TTree * fHitTree
Tree to hold information about every hit in a kalman track.
Vertex location in position and time.
float const & LIDVal() const
#define LOG_WARNING(category)
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
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
CosmicBeamComparison & operator=(CosmicBeamComparison const &)=delete
NueCutParameters fNueSRParams
from the ParametersNumu object
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
float X() const
Definition: RecoHit.h:36
float Mag() const
float const & ShowerNumXCell() const
Definition: ParametersNue.h:77
EventInfo fEventInfo
keep track of run and subrun, and event
TRandom3 r(0)
float const & MinBotDist() const
float const & SliceHits() const
Definition: ParametersNue.h:74
float Y() const
Definition: RecoHit.h:37
#define LOG_VERBATIM(category)
float PECorr() const
Definition: RecoHit.cxx:47
float const & FracShowerHits() const
Definition: ParametersNue.h:83
ValidHandle< PROD > getValidHandle(InputTag const &tag) const
int GoodDiBlockMask(int subrun=-1, bool reload=false)
float const & LEMVal() const
Float_t e
Definition: plot.C:35
RunNumber_t run() const
Definition: Event.h:77
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
float const & MinFrontDist() const
NumuCutParameters fNumuCutParams
from the ParametersNumu object
#define W(x)
TTree * fCosmicHitTree
Tree to hold information about every hit in a kalman track.
void reconfigure(fhicl::ParameterSet const &p)
float const & NumShowers() const
Definition: ParametersNue.h:76
TTree * fNueTree
Tree to hold values of cut parameters for selected events.
EventID id() const
Definition: Event.h:56
Encapsulate the geometry of one entire detector (near, far, ndos)
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
bool failedToGet() const
Definition: Handle.h:196
RunNumber_t run() const
Definition: SubRun.h:49
float const & MinWestDist() const