KalmanTrackAna_module.cc
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////
2 // \file KalmanTrackAna.cxx
3 // \brief Analysis module to make validation plots for reconstructed tracks
4 // \version $Id: KalmanTrackAna_module.cc,v 1.2 2012/12/10 21:28:53 raddatz Exp $
5 // \author Nick Raddatz
6 ////////////////////////////////////////////////////////////////////////////
7 #include <string>
8 #include <vector>
9 
10 // ROOT includes
11 #include "TH1D.h"
12 #include "TH2D.h"
13 #include "TVector3.h"
14 #include "TTree.h"
15 
16 // Framework includes
23 
24 // NOvA includes
26 #include "Calibrator/Calibrator.h"
27 #include "Geometry/Geometry.h"
29 #include "RecoBase/CellHit.h"
30 #include "RecoBase/Cluster.h"
31 #include "RecoBase/Track.h"
32 #include "RecoBase/RecoHit.h"
33 #include "Simulation/Simulation.h"
34 #include "NovaDAQConventions/DAQConventions.h"
35 #include "MCCheater/BackTracker.h"
36 #include "RecoBase/FilterList.h"
37 
38 // Required for CAFAna Cuts
40 #include "CAFAna/Cuts/Cuts.h"
41 #include "CAFAna/Cuts/TruthCuts.h"
42 #include "CAFAna/Cuts/TimingCuts.h"
43 #include "CAFAna/Vars/Vars.h"
44 
47 
48 namespace trk {
49  // A module to analyze tracks produced by KalmanTrack module
51  public:
52 
53  explicit KalmanTrackAna(fhicl::ParameterSet const& pset);
54  virtual ~KalmanTrackAna();
55  void reconfigure(const fhicl::ParameterSet& pset);
56 
57  void beginJob();
58  void analyze(art::Event const& evt);
59 
60  private:
61  void CheckRecoTracks(std::vector< art::Ptr<rb::Track> > const& trkcol,
62  art::PtrVector<rb::CellHit> const& allhits);
63  std::set<int> FindVisibleProngs(const art::PtrVector<rb::CellHit>& allhits,
64  std::set<int> ids) const;
66  geo::View_t view) const;
67  double Containment(TVector3 start) const;
68 
69  void CheckAssociations(art::Ptr<rb::Cluster> slice,std::vector<art::Ptr<rb::Track> > tracks, bool filtered);
70 
71 
72  std::string fSliceLabel; ///< Where to find Slices to reconstruct
73  std::string fTrackLabel; ///< Where to find Tracks to reconstruct
74 
75  std::string fCAFLabel; ///< label for the process that made the standard records
76  bool fApplyCAFCuts; ///< should we apply the CAF cuts?
77  int fCutType; ///< what cuts to apply (see CAFCutter.h)
78 
79  recovalid::CAFCutter fCAFCutter; ///< the CAFCutter helper class
80 
81  // Individual Track Level Histograms
82  TH1D* fTrackLength; ///< Length of reconstructed 3d tracks
83  TH1D* fHitsPerTrack; ///< Number of hits in each reconstructed 3d track
84  TH1D* fTrackLength2d; ///< Length of reconstructed 3d tracks
85  TH1D* fHitsPerTrack2d; ///< Number of hits in each reconstructed 3d track
86  TH1D* fCosX; ///< Cosine theta of 3d tracks with respect to the X axis
87  TH1D* fCosY; ///< Cosine theta of 3D tracks with respect ot the Y axis
88  TH1D* fCosZ; ///< Cosine theta of 3d tracks with respect to the Z axis
89  TH1D* fCosX2d; ///< Cosine theta of 2d tracks with respect to the X axis
90  TH1D* fCosY2d; ///< Cosine theta of 2D tracks with respect ot the Y axis
91  TH1D* fCosZ2d; ///< Cosine theta of 2d tracks with respect to the Z axis
92  TH1D* fXPosStart; ///< x start position of tracks
93  TH1D* fYPosStart; ///< y start position of tracks
94  TH1D* fZPosStart; ///< z start position of tracks
95  TH1D* fXPosEnd; ///< x end position of tracks
96  TH1D* fYPosEnd; ///< y end position of tracks
97  TH1D* fZPosEnd; ///< z end position of tracks
98  TH1D* fXPos; ///< position of x track hits
99  TH1D* fYPos; ///< position of y track hits
100  TH1D* fXRecoPos; ///< reconstructed x position from yz track hit
101  TH1D* fYRecoPos; ///< reconstructed y position from xz track hit
102  TH2D* fXvYStart; ///< Position map for track start position
103  TH2D* fXvYEnd; ///< Position map for track end position
104  TH2D* fXvYReco; ///< X vs Y positions of track hits
105  TH2D* fZvXReco; ///< Z vs X of reconstructed track hit positions
106  TH2D* fZvYReco; ///< Z vs Y of reconstructed track hit positions
107  TH2D* fZvX; ///< Z vs X of track hit positions
108  TH2D* fZvY; ///< Z vs Y of track hit positions
109 
110  // The following histograms are for CAF validation purposes
111  TH1D* fCosXCAF; ///< Cosine theta of tracks with respect to the X axis for CAF comparison
112  TH1D* fCosYCAF; ///< Cosine theta of tracks with respect ot the Y axis for CAF comparison
113  TH1D* fCosZCAF; ///< Cosine theta of tracks with respect to the Z axis for CAF comparison
114  TH1D* fViewCAF; ///< View stored in CAF
115  TH1D* fLengthCAF; ///< Length of reconstructed tracks stored in CAF
116  TH1D* fHitsCAF; ///< Number of hits in each track stored in CAF
117  TH1D* fXPosStartCAF; ///< x start position of tracks in CAF
118  TH1D* fYPosStartCAF; ///< y start position of tracks in CAF
119  TH1D* fZPosStartCAF; ///< z start position of tracks in CAF
120  TH1D* fXPosEndCAF; ///< x end position of tracks in CAF
121  TH1D* fYPosEndCAF; ///< y end position of tracks in CAF
122  TH1D* fZPosEndCAF; ///< z end position of tracks in CAF
123 
124  // Event Level Histograms
125  TH1D* fTracksPerEvent; ///< number of 3d tracks found per event
126  TH1D* f2DTracksPerEvent; ///< number of 2d tracks per event
127  TH1D* fTracksPerSlice; ///< number of 3d tracks found per non noise slice
128  TH1D* f2DTracksPerSlice; ///< number of 2d tracks per non noise slice
129 
130 
131  // Track Hit Level Histograms
132  TH1D* fXHitsPerPlane; ///< number of hits found in each x plane
133  TH1D* fYHitsPerPlane; ///< number of hits found in each y plane
134  TH1D* fHitsPerPlane; ///< number of hits found in each plane
135  TH1D* fXHitsPerCell; ///< number of hits found in each x cell
136  TH1D* fYHitsPerCell; ///< number of hits found in each y cell
137  TH2D* fXPlaneVCell; ///< Plane vs Cell detector map for x view
138  TH2D* fYPlaneVCell; ///< Plane vs Cell detector map for y view
139  TH2D* fPlaneVCell; ///< Plane vs Cell map for whole detector
140 
141  // Tree holding truth information for determining things like completeness, purity, etc... of muons
142  TTree* fTruthTree;
143 
144  int fpdg;
145  float fE, fp, fvtxdst, feff, fpur, flen, fx, fy, fz, fdx, fdy, fdz, fenu;
146  int fmode;
147  int fccnc;
148  float fweight;
149  int fevent;
151  int fnhits;
152  int fview;
153  float frecolen;
154 
155  struct Counts
156  {
158  {
159  nHits[0] = nHits[1] = 0;
160  nOnlyHits[0] = nOnlyHits[1] = 0;
161  }
162 
163  int nHits[2]; // Indexed by view
164  int nOnlyHits[2]; // Hit by only one particle
165  };
166 
167  };
168 
169  //......................................................................
170 
172  EDAnalyzer(pset)
173  {
174  reconfigure(pset);
175  }
176 
177  //......................................................................
178 
180  {
181  }
182 
183  //.....................................................................
185  fSliceLabel = (pset.get< std::string >("SliceLabel") );
186  fTrackLabel = (pset.get< std::string >("TrackLabel") );
187  fCAFLabel = (pset.get< std::string >("CAFLabel") );
188  fApplyCAFCuts = (pset.get< bool > ("ApplyCAFCuts") );
189  fCutType = (pset.get< int > ("CutType") );
190  }
191 
192  //......................................................................
193 
195  {
196 
197  double maxLength = 3000.0;
198  int maxHits = 500;
199  double xMax = 1200.0;
200  double yMax = 1200.0;
201  double zMax = 6500.0;
202  double nCell = 384;
203  double nPlane = 960;
204  int xbins = 480;
205  int ybins = 480;
206  int zbins = 1300;
207 
209  // Individual Track Level Histograms
210  fTrackLength = tfs->make<TH1D>("trackLength",";Track length (cm);Number of tracks",600,0.0,maxLength);
211  fHitsPerTrack = tfs->make<TH1D>("hitsPerTrack",";Number of hits;Number of tracks", maxHits+1,-0.5,(double)maxHits+0.5);
212  fTrackLength2d = tfs->make<TH1D>("trackLength2d",";Track length (cm);Number of 2D tracks",600,0.0,maxLength);
213  fHitsPerTrack2d = tfs->make<TH1D>("hitsPerTrack2d",";Number of hits;Number of 2D tracks", maxHits+1,-0.5,(double)maxHits+0.5);
214  fCosX = tfs->make<TH1D>("dirCosX",";Track direction cos(#theta_{x});Number of tracks", 101,-1.0,1.0);
215  fCosY = tfs->make<TH1D>("dirCosY",";Track direction cos(#theta_{y});Number of tracks", 101,-1.0,1.0);
216  fCosZ = tfs->make<TH1D>("dirCosZ",";Track direction cos(#theta_{z});Number of tracks",100,-1.0,1.0);
217  fCosX2d = tfs->make<TH1D>("dirCosX2d",";Track direction cos(#theta_{x});Number of 2D tracks", 101,-1.0,1.0);
218  fCosY2d = tfs->make<TH1D>("dirCosY2d",";Track direction cos(#theta_{y});Number of tracks", 101,-1.0,1.0);
219  fCosZ2d = tfs->make<TH1D>("dirCosZ2d",";Track direction cos(#theta_{z});Number of tracks",100,-1.0,1.0);
220 
221  fXPosStart = tfs->make<TH1D>("xPosStart",";X start position (cm);Number of tracks",xbins,-xMax,xMax);
222  fYPosStart = tfs->make<TH1D>("yPosStart",";Y start position (cm);Number of tracks",ybins,-yMax,yMax);
223  fZPosStart = tfs->make<TH1D>("zPosStart",";Z start position (cm);Number of tracks",zbins,0.0,zMax);
224  fXPosEnd = tfs->make<TH1D>("xPosEnd",";X end position (cm);Number of tracks",xbins,-xMax,xMax);
225  fYPosEnd = tfs->make<TH1D>("yPosEnd",";y end position (cm);Number of tracks",ybins,-yMax,yMax);
226  fZPosEnd = tfs->make<TH1D>("zPosEnd",";Z end position (cm);Number of tracks",zbins,0.0,zMax);
227  fXvYStart = tfs->make<TH2D>("xyPosStart",";X start position (cm);Y start position (cm)",xbins,-xMax,xMax,ybins,-yMax,yMax );
228  fXvYEnd = tfs->make<TH2D>("xyPosEnd",";X end position (cm);Y end position (cm)",xbins,-xMax,xMax,ybins,-yMax ,yMax);
229  fXvYReco = tfs->make<TH2D>("xyRecoHits",";X position (cm);Y position (cm);Number of hits",xbins,-xMax,xMax,ybins,-yMax,yMax);
230  fZvXReco = tfs->make<TH2D>("zxRecoHits",";Z position (cm);X position (cm);Number of hits",zbins,-200.0,zMax,xbins,-xMax,xMax);
231  fZvYReco = tfs->make<TH2D>("zyRecoHits",";Z position (cm);Y position (cm);Number of hits",zbins,-200.0,zMax,ybins,-yMax,yMax);
232  fZvX = tfs->make<TH2D>("zxCellHits",";Z position (cm);X position (cm);Number of hits",zbins,-200.0,zMax,xbins,-xMax,xMax);
233  fZvY = tfs->make<TH2D>("zyCellHits",";Z position (cm);Y position (cm);Number of hits",zbins,-200.0,zMax,ybins,-yMax,yMax);
234 
235 
236  int ndirbins = 100;
237  double dirmin = -1.1;
238  double dirmax = 1.1;
239  int xbinscaf = 100;
240  int ybinscaf = 100;
241  int zbinscaf = 100;
242  double xmaxcaf = 1000.0;
243  double ymaxcaf = 1000.0;
244  double zmaxcaf = 6000.0;
245  int nlenbins = 100;
246  double nlenmin = -0.5;
247  double nlenmax = 5500;
248  int nhitsbins = 100;
249  double nhitmin = -0.5;
250  double nhitmax = 800.0;
251 
252  fCosXCAF = tfs->make<TH1D>("dirXCAF",";Cos #theta_{X}", ndirbins,dirmin,dirmax);
253  fCosYCAF = tfs->make<TH1D>("dirYCAF",";Cos #theta_{Y}", ndirbins,dirmin,dirmax);
254  fCosZCAF = tfs->make<TH1D>("dirZCAF",";Cos #theta_{Z}",ndirbins,dirmin,dirmax);
255  fViewCAF = tfs->make<TH1D>("viewCAF",";view;",3,-0.5,2.5);
256  fXPosStartCAF = tfs->make<TH1D>("startXCAF",";X Start Position (cm);Number of tracks",xbinscaf,-xmaxcaf,xmaxcaf);
257  fYPosStartCAF = tfs->make<TH1D>("startYCAF",";Y Start Position (cm);Number of tracks",ybinscaf,-ymaxcaf,ymaxcaf);
258  fZPosStartCAF = tfs->make<TH1D>("startZCAF",";Z Start Position (cm);Number of tracks",zbinscaf,0.0,zmaxcaf);
259  fXPosEndCAF = tfs->make<TH1D>("stopXCAF",";X End Position (cm);Number of tracks",xbinscaf,-xmaxcaf,xmaxcaf);
260  fYPosEndCAF = tfs->make<TH1D>("stopYCAF",";Y End Position (cm);Number of tracks",ybinscaf,-ymaxcaf,ymaxcaf);
261  fZPosEndCAF = tfs->make<TH1D>("stopZCAF",";Z End Position (cm);Number of tracks",zbinscaf,0.0,zmaxcaf);
262  fLengthCAF = tfs->make<TH1D>("lenCAF",";Track Length (cm);Number of tracks",nlenbins,nlenmin,nlenmax);
263  fHitsCAF = tfs->make<TH1D>("nhitCAF",";NHits;Number of tracks",nhitsbins,nhitmin,nhitmax);
264 
265  // Event Level Histograms
266  fTracksPerEvent = tfs->make<TH1D>("3dtracksPerEvent",";3D Tracks Per Event;Number of events",101,-0.5,100.5);
267  f2DTracksPerEvent = tfs->make<TH1D>("2dTracksPerEvent",";2D Tracks Per Event;Number of events",101,-0.5,100.5);
268  fTracksPerSlice = tfs->make<TH1D>("3dtracksPerSlice",";3D Tracks Per Event;Number of events",51,-0.5,50.5);
269  f2DTracksPerSlice = tfs->make<TH1D>("2dTracksPerSlice",";2D Tracks Per Event;Number of events",51,-0.5,50.5);
270 
271  //Hit Level Histograms
272  fXHitsPerPlane = tfs->make<TH1D>("xHitsPerPlane",";Plane Number;Number of hits",nPlane,-0.5,nPlane-0.5);
273  fYHitsPerPlane = tfs->make<TH1D>("yHitsPerPlane",";Plane Number;Number of hits",nPlane,-0.5,nPlane-0.5);
274  fHitsPerPlane = tfs->make<TH1D>("hitsPerPlane",";Plane Number;Number of hits",nPlane,-0.5,nPlane-0.5);
275  fXHitsPerCell = tfs->make<TH1D>("xHitsPerCell",";Cell Number;Number of hits",nCell,-0.5,nCell-0.5);
276  fYHitsPerCell = tfs->make<TH1D>("yHitsPerCell",";Cell Number;Number of hits",nCell,-0.5,nCell-0.5);
277  fXPos = tfs->make<TH1D>("xPos",";X Position (cm);Number of track hits", xbins, -xMax, xMax);
278  fYPos = tfs->make<TH1D>("yPos",";Y Position (cm);Number of track hits", ybins, -yMax, yMax);
279  fXRecoPos = tfs->make<TH1D>("xRecoPos",";X Position (cm);Number of track hits", xbins, -xMax, xMax);
280  fYRecoPos = tfs->make<TH1D>("yRecoPos",";Y Position (cm);Number of track hits", ybins, -yMax, yMax);
281  fXPlaneVCell = tfs->make<TH2D>("xPlaneVCell",";Plane number;Cell number",nPlane,-0.5,nPlane-0.5,nCell,-0.5,nCell-0.5);
282  fYPlaneVCell = tfs->make<TH2D>("yPlaneVCell",";Plane number;Cell number",nPlane,-0.5,nPlane-0.5,nCell,-0.5,nCell-0.5);
283  fPlaneVCell = tfs->make<TH2D>("planeVCell",";Plane number;Cell number", nPlane,-0.5,nPlane-0.5,nCell,-0.5,nCell-0.5);
284 
285 
286  // Truth Info
287  fTruthTree = tfs->make<TTree>("truthtree","");
288  fTruthTree->Branch("pdg",&fpdg,"pdg/I");
289  fTruthTree->Branch("E",&fE,"E/F");
290  fTruthTree->Branch("p",&fp,"P/F");
291  fTruthTree->Branch("vtxdst",&fvtxdst,"vtxdst/F");
292  fTruthTree->Branch("eff",&feff,"eff/F");
293  fTruthTree->Branch("pur",&fpur,"pur/F");
294  fTruthTree->Branch("len",&flen,"len/F");
295  fTruthTree->Branch("x",&fx,"x/F");
296  fTruthTree->Branch("y",&fy,"y/F");
297  fTruthTree->Branch("z",&fz,"z/F");
298  fTruthTree->Branch("dx",&fdx,"dx/F");
299  fTruthTree->Branch("dy",&fdy,"dy/F");
300  fTruthTree->Branch("dz",&fdz,"dz/F");
301  fTruthTree->Branch("Enu",&fenu,"Enu/F");
302  fTruthTree->Branch("Mode",&fmode,"Mode/I");
303  fTruthTree->Branch("weight",&fweight,"weight/F");
304  fTruthTree->Branch("event",&fevent,"event/I");
305  fTruthTree->Branch("passMuon",&fpassMuon,"passMuon/I");
306  fTruthTree->Branch("nhits",&fnhits,"nhits/I");
307  fTruthTree->Branch("view",&fview,"view/I");
308  fTruthTree->Branch("recolen",&frecolen,"recolen/F");
309  fTruthTree->Branch("ccnc",&fccnc,"ccnc/I");
310 
311  }
312 
313  //......................................................................
314 
315  ///
316  /// The analysis method is intended to perform checks of the
317  /// performance of the reconstruction for the purposes of
318  /// optimization, testing, and validation.
319  ///
320  /// \param evt - Read-only access to the event data
321  ///
322  // checks to make include -
323  // 1. are the direction cosines of the track correct
324  // 2. fraction of hits used in tracks
325  // 3. number of particles leaving tracks vs number tracks found
326  // 4. fraction of hits in a track that belong to a given particle
327 
329  {
330 
331  //get some information about the cell locations from the Geometry object.
333 
334  // get a backtracker service handle
336 
337  // get all the slices
339  evt.getByLabel(fSliceLabel,slicecol);
340 
341  art::PtrVector<rb::Cluster> slicelist;
342  for(unsigned int i = 0; i<slicecol->size();++i){
343  art::Ptr<rb::Cluster>slice(slicecol,i);
344  slicelist.push_back(slice);
345  }
346 
347  // setup the associations between tracks and slices
348  art::FindManyP<rb::Track> fmtrack(slicecol,evt,fTrackLabel);
349 
350  // Get the FMP for the Standard Records
351  art::FindManyP<caf::StandardRecord> recordFMP(slicecol, evt, fCAFLabel);
352 
353  double n3dTracksPerEvt = 0;
354  double n2dTracksPerEvt = 0;
355 
356  // loop over slices
357  for(unsigned int islice = 0; islice < slicelist.size(); ++islice){
358 
359  // Apply the CAF-level cuts
360  bool pass = true;
361  if( fApplyCAFCuts ) {
362  // get record associated with the slice
363  std::vector< art::Ptr<caf::StandardRecord> > records = recordFMP.at(islice);
364  if(records.size() > 0 && recordFMP.isValid()) {
365  pass = fCAFCutter.passCuts(fCutType, records[0].get());
366  }
367  else { pass = false; }
368  }
369 
370  if(!pass) continue;
371 
372  std::vector<art::Ptr<rb::Track> > tracks;
373  if(fmtrack.isValid()){
374  tracks = fmtrack.at(islice);
375  }
376 
377  bool filtered = rb::IsFiltered(evt,slicelist[islice]);
378 
379  CheckAssociations(slicelist[islice],tracks,filtered);
380 
381  if(slicelist[islice]->IsNoise()){ continue; }
382 
383  double n3dTracksPerSlice = 0;
384  double n2dTracksPerSlice = 0;
385 
386  // get truth info for plots
387  if(bt->HaveTruthInfo()){
388  CheckRecoTracks(tracks,slicelist[islice]->AllCells());
389  }
390 
391  //loop over tracks and make reco track histograms
392  for(unsigned int i=0;i<tracks.size();i++){
393  art::Ptr<rb::Track> track = tracks[i];
394  // check if the track is 2d or 3d
395  fViewCAF->Fill(track->View());
396  fCosXCAF->Fill(track->Dir().X());
397  fCosYCAF->Fill(track->Dir().Y());
398  fCosZCAF->Fill(track->Dir().Z());
399  fXPosStartCAF->Fill(track->Start().X());
400  fYPosStartCAF->Fill(track->Start().Y());
401  fZPosStartCAF ->Fill(track->Start().Z());
402  fXPosEndCAF->Fill(track->Stop().X());
403  fYPosEndCAF->Fill(track->Stop().Y());
404  fZPosEndCAF->Fill(track->Stop().Z());
405  fLengthCAF->Fill(track->TotalLength());
406  fHitsCAF->Fill(track->NCell());
407  if(track->Is3D()){
408  ++n3dTracksPerEvt;
409  ++n3dTracksPerSlice;
410  TVector3 start(track->TrajectoryPoint(0));
411  fXPosStart->Fill(start.X());
412  fYPosStart->Fill(start.Y());
413  fZPosStart->Fill(start.Z());
414  fXvYStart->Fill(start.X(),start.Y());
415  TVector3 end(track->Stop());
416  fXPosEnd->Fill(end.X());
417  fYPosEnd->Fill(end.Y());
418  fZPosEnd->Fill(end.Z());
419  fXvYEnd->Fill(end.X(),end.Y());
420  TVector3 length = end-start;
421  fHitsPerTrack->Fill(track->NCell());
422  fTrackLength->Fill(track->TotalLength());
423  double cosX = (end.X()-start.X())/length.Mag();
424  fCosX->Fill(cosX);
425  double cosY = (end.Y()-start.Y())/length.Mag();
426  fCosY->Fill(cosY);
427  double cosZ = (end.Z()-start.Z())/length.Mag();
428  fCosZ->Fill(cosZ);
429  for(unsigned int nx = 0; nx<track->NXCell();++nx){
430  art::Ptr<rb::CellHit> xcellhit = track->XCell(nx);
431  int cell = xcellhit->Cell();
432  int plane = xcellhit->Plane();
433  double xyz[3];
434  geom->Plane(plane)->Cell(cell)->GetCenter(xyz);
435  fZvX->Fill(xyz[2],xyz[0]);
436  fXPos->Fill(xyz[0]);
437  fXPlaneVCell->Fill(plane,cell);
438  fPlaneVCell->Fill(plane,cell);
439  fXHitsPerPlane->Fill(plane);
440  fHitsPerPlane->Fill(plane);
441  fXHitsPerCell->Fill(cell);
442 
443  //Get the recohit from this track cell hit
444  rb::RecoHit xReco = track->RecoHit(xcellhit);
445  fYRecoPos->Fill(xReco.Y());
446  fXvYReco->Fill(xReco.X(),xReco.Y());
447  fZvYReco->Fill(xReco.Z(),xReco.Y());
448  }
449 
450  for(unsigned int ny = 0; ny<track->NYCell();++ny){
451  art::Ptr<rb::CellHit> ycellhit = track->YCell(ny);
452  int cell = ycellhit->Cell();
453  int plane = ycellhit->Plane();
454  double xyz[3];
455  geom->Plane(plane)->Cell(cell)->GetCenter(xyz);
456  fZvY->Fill(xyz[2],xyz[1]);
457  fYPos->Fill(xyz[1]);
458  fYPlaneVCell->Fill(plane,cell);
459  fPlaneVCell->Fill(plane,cell);
460  fYHitsPerPlane->Fill(plane);
461  fHitsPerPlane->Fill(plane);
462  fYHitsPerCell->Fill(cell);
463 
464  //Get the recohit from this track cell hit
465  rb::RecoHit yReco = track->RecoHit(ycellhit);
466  fXRecoPos->Fill(yReco.X());
467  fXvYReco->Fill(yReco.X(),yReco.Y());
468  fZvXReco->Fill(yReco.Z(),yReco.X());
469  }
470 
471  }
472  else{
473  ++n2dTracksPerEvt;
474  ++n2dTracksPerSlice;
475  if(track->View() == geo::kX){
476  TVector3 start(track->TrajectoryPoint(0));
477  fXPosStart->Fill(start.X());
478  fZPosStart->Fill(start.Z());
479  TVector3 end(track->Stop());
480  fXPosEnd->Fill(end.X());
481  fZPosEnd->Fill(end.Z());
482  TVector3 length = end-start;
483  fHitsPerTrack2d->Fill(track->NCell());
484  fTrackLength2d->Fill(track->TotalLength());
485  double cosX2d = (end.X()-start.X())/length.Mag();
486  fCosX2d->Fill(cosX2d);
487  double cosZ2d = (end.Z()-start.Z())/length.Mag();
488  fCosZ2d->Fill(cosZ2d);
489  }
490  else if(track->View() == geo::kY){
491  TVector3 start(track->TrajectoryPoint(0));
492  fYPosStart->Fill(start.Y());
493  fZPosStart->Fill(start.Z());
494  TVector3 end(track->Stop());
495  fYPosEnd->Fill(end.Y());
496  fZPosEnd->Fill(end.Z());
497  TVector3 length = end-start;
498  fHitsPerTrack2d->Fill(track->NCell());
499  fTrackLength2d->Fill(track->TotalLength());
500  double cosY2d = (end.Y()-start.Y())/length.Mag();
501  fCosY2d->Fill(cosY2d);
502  double cosZ2d = (end.Z()-start.Z())/length.Mag();
503  fCosZ2d->Fill(cosZ2d);
504  }
505  }
506  }
507  fTracksPerSlice->Fill(n3dTracksPerSlice);
508  f2DTracksPerSlice->Fill(n2dTracksPerSlice);
509  }
510 
511  fTracksPerEvent->Fill(n3dTracksPerEvt);
512  f2DTracksPerEvent->Fill(n2dTracksPerEvt);
513 
514  }
515 
516  //-------------------------------------------------------------------
517  void KalmanTrackAna::CheckRecoTracks(std::vector< art::Ptr<rb::Track> > const& trkcol,
518  art::PtrVector<rb::CellHit> const& allhits)
519  {
520 
521  // grab the set of eve IDs for this plist
523  std::set<int> trackIDs = bt->GetTrackIDList();
524 
525  // Only keep the ones associated with not completely unreasonable tracks
526  // remove for now, take care of this with PassMuonCuts function
527  trackIDs = FindVisibleProngs(allhits, trackIDs);
528 
529  // Empty set means accept all PDGs. Otherwise, filter out the unwanted ones
530  std::set<int> newTrackIDs;
531  for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end(); ++it){
532  const int pdg = bt->TrackIDToParticle(*it)->PdgCode();
533  if(abs(pdg) != 13){ continue; }
534  // make sure this is a primary particle
535  if(bt->TrackIDToMotherParticle(*it)->TrackId() != *it){continue;}
536  // check that its a particle we care about
537  newTrackIDs.insert(*it);
538  } // end for it
539  trackIDs = newTrackIDs;
540 
541  std::map<int, int> parents;
542  for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end(); ++it){
543  const sim::Particle* mother = bt->TrackIDToParticle(*it);
544  if(!mother) continue; // Can this happen?
545  for(int i = 0; i < mother->NumberDaughters(); ++i){
546  const int id = mother->Daughter(i);
547  assert(bt->TrackIDToParticle(id));
548  if(abs(bt->TrackIDToParticle(id)->PdgCode()) == 11){ parents[id] = *it; }
549  } // end for i
550  } // end for it
551 
552  // Indexed by what view they're in. X, Y, or 3D
553  assert(geo::kXorY == 2); // Ensure enum still works how we think
554  std::map<int, double> bestPur[3], bestEff[3], bestLength[3], bestHits[3];
555 
556  // Go through each reconstructed object, accumulating the best scores seen
557  // for each trackID we're interested in.
558  for(size_t c = 0; c < trkcol.size(); ++c){
559  const rb::Track* trk = trkcol[c].get();
560  const geo::View_t view = trk->View();
561 
562  // get the hits associated with this object
564 
565  // use the cheat::BackTracker to find purity and efficiency for these
566  // hits for each particle ID
567  std::map<int, double> purMap, effMap;
568  bt->HitCollectionPurity(trackIDs, hits, &purMap, &parents);
569  bt->HitCollectionEfficiency(trackIDs, hits, allhits, view, &effMap);
570 
571  for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end(); ++it){
572  const double p = purMap[*it];
573  const double e = effMap[*it];
574 
575  if(e > bestEff[view][*it]) {
576  bestEff[view][*it] = e;
577  bestPur[view][*it] = p;
578 
579  // If we were told to fill a track length histogram then downcast
580  // should be safe.
581  bestLength[view][*it] = trk->TotalLength();
582  bestHits[view][*it] = trk->NCell();
583  }
584  } // end for it
585  } // end for c
586 
587  // For every track. If it was matched in 3D, fill the histograms.
588  // If there is no 3D match look for 2 2D matches and fill with weight 0.5
589  // If required we can add a flag to only try 3D or only try 2D.
590  for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end(); ++it){
591  // Try 3D first
592  for(int view = geo::kXorY; view >= geo::kX; --view){
593  // If no 3D match, continue to two 2D options
594  if(view == geo::kXorY && bestEff[view].find(*it) == bestEff[view].end()) continue;
595 
596  const double weight = (view == geo::kXorY) ? 1 : 0.5;
597 
598  const double trueLength = TotalLengthInDetector(bt->TrackIDToParticle(*it), geo::View_t(view));
599 
600  const sim::Particle* part = bt->TrackIDToParticle(*it);
601  fE = part->E();
602  fp = part->P();
603  TVector3 truestart;
604  truestart.SetXYZ(part->Vx(),part->Vy(),part->Vz());
605  fx = part->Vx();
606  fy = part->Vy();
607  fz = part->Vz();
608  fvtxdst = Containment(truestart);
609  fpdg = part->PdgCode();
610  fdx = part->Px()/part->P();
611  fdy = part->Py()/part->P();
612  fdz = part->Pz()/part->P();
613  fnhits = bestHits[view][*it];
614  feff = bestEff[view][*it];
615  fpur = bestPur[view][*it];
616  frecolen = bestLength[view][*it];
617  flen = trueLength;
618  fweight = weight;
619  fpassMuon = bt->PassMuonCuts(*it,allhits);
620  fview = view;
621  fmode = -1;
622  fenu = -1;
623 
624  // check if this is coming from a neutrino
625  const art::Ptr<simb::MCTruth> mctruth = bt->ParticleToMCTruth(part);
626  if(mctruth->NeutrinoSet()){
627  // get the neutrino if it exists
628  const simb::MCNeutrino neutrino = mctruth->GetNeutrino();
629  fmode = neutrino.Mode();
630  fenu = neutrino.Nu().E();
631  fccnc = neutrino.CCNC();
632  }
633 
634  fTruthTree->Fill();
635 
636  // If we did find a 3D match, don't try 2D
637  if(bestEff[geo::kXorY].find(*it) != bestEff[geo::kXorY].end()) break;
638  } // end for view
639  } // end for it
640  }
641 
642  //-------------------------------------------------------------------
644  std::set<int> ids) const
645  {
646  std::map<int, Counts> counts;
647 
649 
650  // loop over all hits and check if this track id contributed to them
651  for(unsigned int i = 0; i < allhits.size(); ++i){
652  //get vector of track id's that contribute to this hit
653  const art::Ptr<rb::CellHit> hit = allhits[i];
654  const std::vector<cheat::TrackIDE> trackIDs = bt->HitToTrackIDE(hit);
655  // loop over all the returned track ids
656  int nTracks = 0; // number of tracks that contributed energy to the hit
657  int loneID = -1; // Which ID it was if there was only one
658  for(unsigned int j = 0; j < trackIDs.size(); ++j){
659  const int trackID = trackIDs[j].trackID;
660  if(trackID != sim::kNoiseId){
661  ++nTracks;
662  loneID = trackID;
663  }
664 
665  // check if the current id matches the input track id
666  if(ids.find(trackID) != ids.end()){
667  ++counts[trackID].nHits[hit->View()];
668  }
669  } // end for j
670  if(nTracks == 1){
671  ++counts[loneID].nOnlyHits[hit->View()];
672  }
673  } // end for i
674 
675  const int kMinHits = 2; // minimum number of hits to be considered visible
676  const int kMinOnlyHit = 1; // minimum number of hits in which it is only particle contributing energy to be visible
677 
678  std::set<int> ret;
679  for(std::set<int>::iterator it = ids.begin(); it != ids.end(); ++it){
680  const Counts c = counts[*it];
681  // Check if this track passes the visible prong requirements
682  if(c.nOnlyHits[geo::kX] >= kMinOnlyHit &&
683  c.nOnlyHits[geo::kY] >= kMinOnlyHit){
684 
685  if(c.nHits[geo::kX] >= kMinHits && c.nHits[geo::kY] >= kMinHits)
686  ret.insert(*it);
687  }
688  }
689  return ret;
690  }
691 
692  //-------------------------------------------------------------------
694  geo::View_t view) const
695  {
696  const unsigned int N = p->NumberTrajectoryPoints();
697 
698  if(N < 2){
699  mf::LogWarning("RecoCheckAna") << "Particle has no trajectory points";
700  return 0;
701  }
702 
704 
705  double dist = 0;
706 
707  for(unsigned int n = 0; n < N-1; ++n){
708  TVector3 p0 = p->Position(n).Vect();
709  TVector3 p1 = p->Position(n+1).Vect();
710  if(p0.Z() < 0 || p1.Z() < 0 ||
711  p0.Z() > geom->DetLength() || p1.Z() > geom->DetLength() ||
712  fabs(p0.X()) > geom->DetHalfWidth() || fabs(p1.X()) > geom->DetHalfWidth() ||
713  fabs(p0.Y()) > geom->DetHalfHeight() || fabs(p1.Y()) > geom->DetHalfHeight()) continue;
714 
715  if(view == geo::kX){
716  p0.SetY(0);
717  p1.SetY(0);
718  }
719  if(view == geo::kY){
720  p0.SetX(0);
721  p1.SetX(0);
722  }
723  dist += (p1-p0).Mag();
724  }
725 
726  return dist;
727  }
728 
729  //.....................................................................
730  double KalmanTrackAna::Containment(TVector3 start) const
731  {
732  double dist = 0;
733  double xmin,xmax,ymin,ymax,zmin,zmax;
734 
735  // define the boundaries of containment
737 
738  double detHalfWidth = geom->DetHalfWidth();
739  double detHalfHeight = geom->DetHalfHeight();
740  double detLength = geom->DetLength();
741 
742  // fd detector edges determined from event display for 28 block detector
743  // get rid of offset between geom->DetHalfWidth()(Height()) and actual detector edges
744  double voffset = 15.0;
745  xmin = -detHalfWidth+voffset;
746  xmax = detHalfWidth-voffset;
747  ymin = -detHalfHeight+voffset;
748  ymax = detHalfHeight-voffset;
749  zmin = 0;
750  zmax = detLength;
751 
752  TVector3 mins(xmin,ymin,zmin);
753  TVector3 maxs(xmax,ymax,zmax);
754 
755  float mindist = 1e10;
756 
757  if(start.X() > xmin && start.X() < xmax &&
758  start.Y() > ymin && start.Y() < ymax &&
759  start.Z() > zmin && start.Z() < zmax){
760  TVector3 mindiffs = mins - start;
761  TVector3 maxdiffs = maxs - start;
762  if(TMath::Abs(mindiffs.X()) < mindist){ mindist = TMath::Abs(mindiffs.X()); }
763  if(TMath::Abs(mindiffs.Y()) < mindist){ mindist = TMath::Abs(mindiffs.Y()); }
764  if(TMath::Abs(mindiffs.Z()) < mindist){ mindist = TMath::Abs(mindiffs.Z()); }
765  if(TMath::Abs(maxdiffs.X()) < mindist){ mindist = TMath::Abs(maxdiffs.X()); }
766  if(TMath::Abs(maxdiffs.Y()) < mindist){ mindist = TMath::Abs(maxdiffs.Y()); }
767  if(TMath::Abs(maxdiffs.Z()) < mindist){ mindist = TMath::Abs(maxdiffs.Z()); }
768  dist = mindist;
769  }
770 
771  return dist;
772  }
773 
774  //.....................................................................
775  void KalmanTrackAna::CheckAssociations(art::Ptr<rb::Cluster> slice,std::vector<art::Ptr<rb::Track> > tracks, bool filtered)
776  {
777  if(filtered && tracks.size() > 0){
778  // mf::LogVerbatim("KalmanTrackAna")<<"Filtered slice has non-zero number of reconstructed tracks. This should never happen.";
779  // Should never create tracks on slices that have been filtered
780  // abort();
781  }
782  if(slice->IsNoise() && tracks.size() > 0){
783  mf::LogVerbatim("KalmanTrackAna")<<"Noise slice has non-zero number of reconstructed tracks. This should never happen.";
784  // Should never create tracks on slices that have been filtered
785  abort();
786  }
787 
788  // check that reco'd track hits all belong to the slice hits
789  art::PtrVector<rb::CellHit> slicehits = slice->AllCells();
790  for(unsigned int itrack = 0; itrack < tracks.size(); ++itrack){
791  // get all of the track hits for this track
792  art::PtrVector<rb::CellHit> trackhits = tracks[itrack]->AllCells();
793  // loop over the track hits
794  for(unsigned int itrhit = 0; itrhit < trackhits.size(); ++itrhit){
795  bool foundhit = false;
796  // loop over all the slice hits
797  for(unsigned int islhit = 0; islhit < slicehits.size(); ++islhit){
798  if(trackhits[itrhit] == slicehits[islhit]){
799  foundhit = true;
800  break;
801  }
802  }
803  // after looking at all the slice hits, did we find this track hit. If not, yell.
804  if(!foundhit){
805  mf::LogVerbatim("KalmanTrackAna")<<"Track made from hits on a slice other than the slice the track is associated to. This should never happen.";
806  // Should never create tracks from hits on slices other than the slice its associated with
807  abort();
808  }
809  }
810  }
811 
812  }
813 
814 } // end namespace trk
815 
double E(const int i=0) const
Definition: MCParticle.h:232
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
void reconfigure(const fhicl::ParameterSet &pset)
TH2D * fXPlaneVCell
Plane vs Cell detector map for x view.
double Containment(TVector3 start) const
TH1D * fXPosEndCAF
x end position of tracks in CAF
TH1D * fYHitsPerCell
number of hits found in each y cell
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
TH1D * fYPosEnd
y end position of tracks
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:217
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
TH1D * fXRecoPos
reconstructed x position from yz track hit
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
std::map< std::string, double > xmax
set< int >::iterator it
art::Ptr< rb::CellHit > XCell(unsigned int xIdx) const
Get the ith cell in the x-view.
Definition: Cluster.cxx:157
TH2D * fXvYEnd
Position map for track end position.
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
TH2D * fYPlaneVCell
Plane vs Cell detector map for y view.
KalmanTrackAna(fhicl::ParameterSet const &pset)
const sim::Particle * TrackIDToMotherParticle(int const &id) const
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
static const int kNoiseId
flag for noise id
Definition: Simulation.h:13
recovalid::CAFCutter fCAFCutter
the CAFCutter helper class
TH1D * fTrackLength
Length of reconstructed 3d tracks.
TH1D * fLengthCAF
Length of reconstructed tracks stored in CAF.
std::set< int > FindVisibleProngs(const art::PtrVector< rb::CellHit > &allhits, std::set< int > ids) const
const Var weight
X or Y views.
Definition: PlaneGeo.h:30
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
TH1D * fZPosStartCAF
z start position of tracks in CAF
geo::View_t View() const
Definition: CellHit.h:41
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
TH1D * fCosZ
Cosine theta of 3d tracks with respect to the Z axis.
TH1D * fTracksPerEvent
number of 3d tracks found per event
TH2D * fPlaneVCell
Plane vs Cell map for whole detector.
Vertical planes which measure X.
Definition: PlaneGeo.h:28
float Z() const
Definition: RecoHit.h:38
Definition: event.h:19
double DetLength() const
std::string fSliceLabel
Where to find Slices to reconstruct.
TH1D * f2DTracksPerSlice
number of 2d tracks per non noise slice
void abs(TH1 *hist)
TH1D * fCosY2d
Cosine theta of 2D tracks with respect ot the Y axis.
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
TH1D * fYPosStart
y start position of tracks
const PlaneGeo * Plane(unsigned int i) const
DEFINE_ART_MODULE(TestTMapFile)
TH1D * fXHitsPerCell
number of hits found in each x cell
TH1D * fCosZ2d
Cosine theta of 2d tracks with respect to the Z axis.
virtual TVector3 Start() const
Definition: Prong.h:73
TH1D * fXPosStart
x start position of tracks
int NumberDaughters() const
Definition: MCParticle.h:216
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
TH1D * fXPosStartCAF
x start position of tracks in CAF
double HitCollectionEfficiency(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, const std::vector< const rb::CellHit * > &allhits, const geo::View_t &view, std::map< int, double > *effMap=0, bool energyEff=false, double *desiredEnergy=0, double *totalEnergy=0, int *desiredHits=0, int *totalHits=0) const
Returns the fraction of all energy in an event from a specific set of Geant4 track IDs that are repre...
int TrackId() const
Definition: MCParticle.h:209
int Daughter(const int i) const
Definition: MCParticle.cxx:112
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
TH1D * fYRecoPos
reconstructed y position from xz track hit
TH1D * f2DTracksPerEvent
number of 2d tracks per event
double dist
Definition: runWimpSim.h:113
Double_t ymax
Definition: plot.C:25
TH1D * fHitsCAF
Number of hits in each track stored in CAF.
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
unsigned short Cell() const
Definition: CellHit.h:40
TH1D * fZPosStart
z start position of tracks
Track finder for cosmic rays.
TString part[npart]
Definition: Style.C:32
TH1D * fXPos
position of x track hits
TH1D * fYPosStartCAF
y start position of tracks in CAF
bool passCuts(int cut, const caf::StandardRecord *sr)
Definition: CAFCutter.cxx:37
void hits()
Definition: readHits.C:15
length
Definition: demo0.py:21
TH1D * fCosY
Cosine theta of 3D tracks with respect ot the Y axis.
const int xbins
Definition: MakePlots.C:82
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
TH1D * fCosX
Cosine theta of 3d tracks with respect to the X axis.
void CheckAssociations(art::Ptr< rb::Cluster > slice, std::vector< art::Ptr< rb::Track > > tracks, bool filtered)
art::Ptr< rb::CellHit > YCell(unsigned int yIdx) const
Get the ith cell in the y-view.
Definition: Cluster.cxx:165
TH1D * fXPosEnd
x end position of tracks
int fCutType
what cuts to apply (see CAFCutter.h)
TH1D * fZPosEndCAF
z end position of tracks in CAF
TH1D * fCosX2d
Cosine theta of 2d tracks with respect to the X axis.
T get(std::string const &key) const
Definition: ParameterSet.h:231
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
int evt
std::string fTrackLabel
Where to find Tracks to reconstruct.
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
TH1D * fTrackLength2d
Length of reconstructed 3d tracks.
const art::Ptr< simb::MCTruth > & ParticleToMCTruth(const sim::Particle *p) const
TH1D * fHitsPerPlane
number of hits found in each plane
TH1D * fHitsPerTrack2d
Number of hits in each reconstructed 3d track.
TH1D * fCosYCAF
Cosine theta of tracks with respect ot the Y axis for CAF comparison.
TH2D * fZvX
Z vs X of track hit positions.
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
const double j
Definition: BetheBloch.cxx:29
TH2D * fXvYReco
X vs Y positions of track hits.
TH1D * fYPosEndCAF
y end position of tracks in CAF
double DetHalfHeight() const
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
TH2D * fZvXReco
Z vs X of reconstructed track hit positions.
bool fApplyCAFCuts
should we apply the CAF cuts?
void analyze(art::Event const &evt)
size_type size() const
Definition: PtrVector.h:308
double HitCollectionPurity(const std::set< int > &trackIDs, const std::vector< rb::WeightedHit > &whits, std::map< int, double > *purMap=0, std::map< int, int > *parents=0, bool energyPur=false) const
Returns the fraction of hits in a collection that come from the specified Geant4 track ids...
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
const int ybins
std::string fCAFLabel
label for the process that made the standard records
TH1D * fZPosEnd
z end position of tracks
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
double DetHalfWidth() const
TH1D * fViewCAF
View stored in CAF.
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
TH1D * fYHitsPerPlane
number of hits found in each y plane
void CheckRecoTracks(std::vector< art::Ptr< rb::Track > > const &trkcol, art::PtrVector< rb::CellHit > const &allhits)
Definition: structs.h:12
float X() const
Definition: RecoHit.h:36
TH1D * fXHitsPerPlane
number of hits found in each x plane
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:133
float Mag() const
const std::set< int > GetTrackIDList() const
Get all G4 track ids present in the event.
Definition: BackTracker.h:750
assert(nhit_max >=nhit_nbins)
Double_t ymin
Definition: plot.C:24
double TotalLengthInDetector(const sim::Particle *p, geo::View_t view) const
TH2D * fXvYStart
Position map for track start position.
Helper class for Reco Validation modules.
Definition: CAFCutter.h:58
bool PassMuonCuts(int trackID, art::PtrVector< rb::CellHit > const &hits) const
Tool for NumuEAna that allows one to see if primary muon (or any track ID you feed the function) cont...
TH2D * fZvY
Z vs Y of track hit positions.
float Y() const
Definition: RecoHit.h:37
bool NeutrinoSet() const
Definition: MCTruth.h:77
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Int_t trackID
Definition: plot.C:84
TH1D * fHitsPerTrack
Number of hits in each reconstructed 3d track.
Float_t e
Definition: plot.C:35
Event generator information.
Definition: MCNeutrino.h:18
int Mode() const
Definition: MCNeutrino.h:149
TH2D * fZvYReco
Z vs Y of reconstructed track hit positions.
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
Definition: Prong.h:71
TH1D * fCosXCAF
Cosine theta of tracks with respect to the X axis for CAF comparison.
TH1D * fTracksPerSlice
number of 3d tracks found per non noise slice
const Binning zbins
TH1D * fYPos
position of y track hits
TH1D * fCosZCAF
Cosine theta of tracks with respect to the Z axis for CAF comparison.