CosmicTrackAna_module.cc
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////
2 /// \brief Cosmic track analysis module CosmicTrackAna.
3 /// Makes some validation plots.
4 /// \author Gavin S. Davies - gsdavies@iastate.edu
5 /// \date
6 //////////////////////////////////////////////////////////////////////////
7 #include <string>
8 #include <vector>
9 #include <limits>
10 #include <algorithm>
11 #include <initializer_list>
12 
13 // ROOT includes
14 #include "TH1D.h"
15 #include "TH2D.h"
16 #include "TTree.h"
17 #include "TVector3.h"
18 #include "TSpline.h"
19 #include "TStopwatch.h"
20 
21 // NOvA includes
22 #include "Calibrator/Calibrator.h"
24 #include "GeometryObjects/Geo.h"
26 #include "Geometry/Geometry.h"
27 #include "Geometry/LiveGeometry.h"
28 #include "MCCheater/BackTracker.h"
29 #include "NovaDAQConventions/DAQConventions.h"
30 #include "RecoBase/CellHit.h"
31 #include "RecoBase/RecoHit.h"
32 #include "RecoBase/Track.h"
33 #include "Simulation/Particle.h"
34 #include "Simulation/FLSHit.h"
37 #include "Utilities/func/MathUtil.h"
42 
43 // Framework includes
50 #include "art_root_io/TFileService.h"
51 #include "art_root_io/TFileDirectory.h"
54 #include "fhiclcpp/ParameterSet.h"
57 
59  kDCosX = 0, ///< direction cosine in x
60  kDCosY, ///< direction cosine in y
61  kDCosZ, ///< direction cosine in z
62  kEndX, ///< track end position in x
63  kEndY, ///< track end position in y
64  kEndZ, ///< track end position in z
65  kStartX, ///< track start position in x
66  kStartY, ///< track start position in y
67  kStartZ, ///< track start position in z
68  kDiffX, ///< difference in track start and end position in x
69  kDiffY, ///< difference in track start and end position in y
70  kDiffZ, ///< difference in track start and end position in z
71  kFractionHitsUsed, ///< Fraction of hits in event used in track
72  kFractionSigUsed, ///< Fraction of hits in event used in track
73  kTrkCells, ///< number of cells in a track
74  kTracksFound, ///< number of tracks found
75  kResidualsX, ///< residuals for the tracks
76  kResidualsY, ///< residuals for the tracks
77  kCellHitSig, ///< signal size of calibrated cell hits
78  kTrkSig, ///< track signal (Q)
79  kVtxDistToEdge, ///< distance in cm to closest edge of detector at track vertex
80  kEndDistToEdge, ///< distance in cm to closest edge of detector at track end
81  kHitXPosition, ///< position of reco hits in x
82  kHitYPosition, ///< position of reco hits in y
83  kHitZPosition, ///< position of reco hits in z
84  kTrackLength, ///< length of track
85  kCellHitX, ///<
86  kCellHitY, ///<
87  kPlaneHitX, ///<
88  kPlaneHitY, ///<
89  kCosCosmic, ///< cosine theta of cosmic track with respect to zenith
90  kAzimuth, ///< azimuthal angle of the cosmic ray
91  kTNSX, ///< distribution of x/y cellhit times in nano seconds
92  kTNSY, ///< distribution of x/y cellhit times in nano seconds
93  kTrackRangeMomentum, ///< momentum of track from range
94  kCellsPerPlaneX, ///< average number of cells hit per plane in x
95  kCellsPerPlaneY, ///< average number of cells hit per plane in y
96  kStepRelativeToMedian, ///< size of steps along trajectory normalized by the median step
98 };
99 
101  kDCosXDiff = 0, ///< Fraction difference in direction cosine in x
102  kDCosYDiff, ///< Fraction difference in direction cosine in y
103  kDCosZDiff, ///< Fraction difference in direction cosine in z
105  kPrimaryTrueEnergy, ///< True energy of particle zero (GeV)
106  kDeltaTrackLength, ///< difference between true and reco length in cm
107  kDeltaRangeMomentum, ///< difference in reco momentum from range and true momentum
108  kDeltaRangeMomentumUncor, ///< difference in reco momentum from range uncorrected for track length
109  kDeltaZEnd, ///< difference in z end point of reconstructed track vs simulated
110  kDeltaEnd, ///< difference in end point of reconstructed track vs simulated
111  kDeltaBeg, ///< difference in start point of reconstructed track vs simulated
112  kTrackEfficiency, ///< efficiency in selecting hits for the track
113  kTrackPurity, ///< purity in selecting hits for the track
114  kEndTrueX, ///< truth end position in x
115  kEndTrueY, ///< truth end position in y
116  kEndTrueZ, ///< truth end position in z
117  kStartTrueX, ///< truth start position in x
118  kStartTrueY, ///< truth start position in y
119  kStartTrueZ, ///< truth start position in z
120  kCosCosmicTrue, ///< cosine theta of cosmic track with respect to zenith
121  kAzimuthTrue, ///< azimuthal angle of the cosmic ray
123 };
124 
126  kFracHitsUsedVsNum = 0, ///< Fraction hits used vs the number of hits that could be used
127  kEndXY, ///< track end position in xy
128  kEndZX, ///< track end position in zx
129  kEndZY, ///< track end position in zy
130  kStartXY, ///< track start position in xy
131  kStartZX, ///< track start position in zx
132  kStartZY, ///< track start position in zy
133  kDirXY, ///< direction in x versus y
134  kResidVsCosX, ///< residual for each hit vs cos(theta)
136  kResidVsAzX, ///< residual for each hit vs azimuth
138  kHitXYPosition, ///< reco hit positions in xy
139  kHitXZPosition, ///< reco hit positions in xz
140  kHitYZPosition, ///< reco hit positions in yz
143  kCellsPerPlaneXY, ///< cells hit per plane in Y vs X
144  kPEPerCMVsWX, ///< PE/cm deposited as a function of W in X view
145  kPEPerCMVsWY, ///< PE/cm deposited as a function of W in Y view
146  kStepRelativeToMedianVsAz, ///< largest trajectory step relative to the median vs azimuth
148 };
149 
151  kDCosXRecoVsTrue = 0, ///< Direction cosine in x
152  kDCosYRecoVsTrue, ///< Direction cosine in y
153  kDCosZRecoVsTrue, ///< Direction cosine in z
154  kDCosXDiffVsZLen, ///< Fraction difference in direction cosine in x vs z length
155  kDCosYDiffVsZLen, ///< Fraction difference in direction cosine in y vs z length
156  kDCosZDiffVsZLen, ///< Fraction difference in direction cosine in z vs z length
158  kTracksFoundVsLeft, ///< number of tracks found in event vs particles leaving tracks
159  kDeltaZEndVsDCosZ, ///< difference in z endpoint of recoed track & simulated vs dcosZ
160  kDeltaZBegVsDCosZ, ///< difference in z start point of recoed track & simulated vs dcosZ
161  kDeltaZEndVsDistToEdge, ///< difference in z endpoint of recoed track & simulated vs dcosZ
162  kDiffEndVsDCosZ, ///< difference in end point of recoed track & simulated vs dcosZ
163  kDiffBegVsDCosZ, ///< difference in start point of recoed track & simulated vs dcosZ
164  kDiffEndVsAzimuth, ///< difference in end point of recoed track & simulated vs azimuth
165  kDiffBegVsAzimuth, ///< difference in start point of recoed track & simulated vs azimuth
166  kDiffEndVsCosTheta, ///< difference in end point of recoed track & simulated vs cosTheta
167  kDiffBegVsCosTheta, ///< difference in start point of recoed track & simulated vs cosTheta
168  kDeltaRangeMomVsLength, ///< difference in reco and true range momentum as a function of reco length
169  kDeltaRangeMomUncorVsLength, ///< difference in reco and true range momentum as a function of reco length
170  kDeltaRangeMomUncorVsCosTheta,///< difference in reco and true range momentum as a function of cos theta
171  kDeltaRangeMomUncorVsAzimuth, ///< difference in reco and true range momentum as a function of azimuth
172  kDeltaLengthVsLength, ///< difference in recon and true length as a function of reco length
173  kdEdxVsRange, ///< dE/dx vs remaining length in a track
174  kDeltaLengthVsPlaneAsymmetry, ///< fractional length vs asymmetry in number of planes in each view
175  kDeltaLengthVsDeltaStartPlane,///< fractional length vs difference in start plane in each view
176  kDeltaLengthVsDeltaEndPlane, ///< fractional length vs difference in end plane in each view
177  kDeltaLengthVsCellsPerPlane, ///< fractional length vs maximum cells per plane
179 };
180 
181 
182 namespace trk {
183  // A module to analyse cosmic tracks produced by CosmicTrack module
185  public:
186  explicit CosmicTrackAna(fhicl::ParameterSet const& pset);
187  virtual ~CosmicTrackAna();
188 
189  void analyze(art::Event const& evt);
190  void reconfigure(const fhicl::ParameterSet& pset);
191 
192  virtual void beginSubRun(const art::SubRun& sr);
193 
194  void beginJob();
195  void endJob();
196 
197  private:
198 
199  void Make1DRecoHistograms(std::vector<TH1D *> &hists,
200  TString const sample);
201  void Make1DTruthHistograms(std::vector<TH1D *> &hists,
202  TString const sample);
203  void Make2DRecoHistograms(std::vector<TH2D *> &hists,
204  TString const sample);
205  void Make2DTruthHistograms(std::vector<TH2D *> &hists,
206  TString const sample);
207 
208  std::pair<CosmicTrackInfo, CosmicSignalHitFraction> FillRecoInfo(rb::Track const& track,
209  rb::Cluster const& slice);
211  rb::Cluster const& slice,
212  std::map<geo::CellUniqueId, float> & trueCellPathLength);
213 
214  void FillEventMCTruthHistograms(size_t trkMax);
215 
217  CosmicTrackInfo const& ri);
218  void FillTrackHistograms(CosmicTrackInfo const& ri,
219  CosmicSignalHitFraction const& shf,
220  rb::Track const& track,
221  std::map<geo::CellUniqueId, float> const& trueCellPathLength);
222 
224 
225  CosmicTrackSelection fSelection; ///< object encapsulating selection criteria
226  CosmicTrackUtilities fUtilities; ///< object encapsulating utilitie methods
227  std::string fTrackLabel; ///< Where to find Tracks to reconstruct
228  std::string fExposureLabel; ///< Where to find CosmicExposure info
229  double fMinSig; ///< Minimum signal of cell hit to use in a
230  ///< track - using PE values
231  bool fFindPurityAndEfficiency; ///< Set to true in fcl file to calculate them,
232  ///< false by default
233  ///< because the calculation is so slow
234  //int fTicks; ///< # of clock ticks to read out one event
235  //int fTickWindow; ///< Size of window for ticks of TDC to slice event in time
236  double fLiveTime; ///< Livetime in sec according to file
237  unsigned int fNevents; ///< total number of events counter
238  int fNTracks; ///< total number of reconstructed tracks counter
239  int fNTracksMC; ///< total number of MC particles
240  int fRun; ///< Run number
241  int fSubRun; ///< Subrun number
242  double fCellHalfDepthZ; ///< extent of a cell in the z direction
243  double fDetHalfWidth; ///< detector half width, cached
244  double fDetHalfHeight; ///< detector half height, cached
245  double fDetLength; ///< detector length, cached
246 
247  TH1D* fCosmicRate; ///< histogram holding single value for cosmic rate in Hz at end of job
248  TH1D* fExposure; ///< histogram holding single value for exposure in s at end of job
249 
250  std::vector<TH1D *> fOneDHistsRecoAll; ///< one dimensional histograms for all tracks
251  std::vector<TH1D *> fOneDHistsRecoStop; ///< one dimensional histograms for stopping tracks
252  std::vector<TH1D *> fOneDHistsTruthAll; ///< one dimensional histograms for all tracks
253  std::vector<TH1D *> fOneDHistsTruthStop; ///< one dimensional histograms for stopping tracks
254  std::vector<TH2D *> fTwoDHistsRecoAll; ///< two dimensional histograms for all tracks
255  std::vector<TH2D *> fTwoDHistsRecoStop; ///< two dimensional histograms for stopping tracks
256  std::vector<TH2D *> fTwoDHistsTruthAll; ///< two dimensional histograms for all tracks
257  std::vector<TH2D *> fTwoDHistsTruthStop; ///< two dimensional histograms for stopping tracks
258 
259  CosmicTrackTree fTrigger; ///< event information
260  CosmicTrackInfo fTrackTruth; ///< single track mc truth information
261  CosmicTrackInfo fTrackReco; ///< single track reco information
262  CosmicSignalHitFraction fSigHitFrac; ///< fraction of signal and hits used
263  CosmicHitInfo fCosmicHit; ///< information about hits
264  TTree* fTree; ///< initializer of TTree for this module
265  TTree* fHitTree; ///< initializer of Hit TTree for this module
270 
271  };
272 }
273 
274 //static const double ns = 1.0E-9;
275 
276 namespace trk{
277 
278  bool sortTVector3ByZ(TVector3 const& a, TVector3 const& b)
279  {
280  return a.Z() < b.Z();
281  }
282 
283  bool sortTVector3ByY(TVector3 const& a, TVector3 const& b)
284  {
285  return a.Y() > b.Y();
286  }
287 
288  bool sortFLSHits(sim::FLSHit const& a, sim::FLSHit const& b)
289  {
290  if(a.GetPlaneID() != b.GetPlaneID())
291  return a.GetPlaneID() < b.GetPlaneID();
292 
293  return a.GetCellID() < b.GetCellID();
294  }
295 
296  //......................................................................
298  : EDAnalyzer(pset)
299  , fMinSig (0)
300  , fLiveTime (0)
301  , fNevents (0)
302  , fNTracks (0)
303  , fNTracksMC (0)
304  , fRun (-999)
305  , fSubRun (-999)
306  , fCellHalfDepthZ (-999.)
307  , fDetHalfWidth (-999.)
308  , fDetHalfHeight (-999.)
309  , fDetLength (-999.)
310  {
311  reconfigure(pset);
312  }
313 
314  //......................................................................
316  {
317  }
318 
319  //......................................................................
321  {
322 
323  MF_LOG_DEBUG("CosmicTrackAna") << "cell half depth: " << fGeom->Plane(0)->Cell(0)->HalfD() << "\n"
324  << "cell half width: " << fGeom->Plane(0)->Cell(0)->HalfW() << "\n"
325  << "cell half leng: " << fGeom->Plane(0)->Cell(0)->HalfL();
326 
327  // will always want the reco histograms, not sure we will want the truth
328  // histograms until we determine we have MC
329  this->Make1DRecoHistograms(fOneDHistsRecoAll, "AllReco");
330  this->Make1DRecoHistograms(fOneDHistsRecoStop, "StopReco");
331  this->Make2DRecoHistograms(fTwoDHistsRecoAll, "AllReco");
332  this->Make2DRecoHistograms(fTwoDHistsRecoStop, "StopReco");
333 
334  fCosmicRate = fTFS->make<TH1D>("cosmicrate", ";Cosmic Rate (kHz)", 500, 0., 100);
335  fExposure = fTFS->make<TH1D>("exposure", ";Exposure (s)", 1000, 0., 10);
336 
337 
338  std::string trigst("run/I:subrun/I:event/I:nevents/I:ntracks/I");
339  std::string shfst ("signalFraction/F:hitFraction/F:signal/F:totcell/F:hitsX/F:hitsY/F:signalX/F:signalY/F");
340  std::string ctist ("efficiency/F:purity/F:pdg/I");
341  ctist.append(":length/F:momentum/F:momentumCor/F:momentumBPF/F:calorE/F");
342  ctist.append(":dcosX/F:dcosY/F:dcosZ/F");
343  ctist.append(":vtxX/F:vtxY/F:vtxZ/F");
344  ctist.append(":endX/F:endY/F:endZ/F");
345  ctist.append(":deltaVtxPos/F:deltaEndPos/F");
346  ctist.append(":cosTheta/F:azimuth/F");
347  ctist.append(":vtxToEdge/F:endToEdge/F");
348  ctist.append(":cellsPerPlaneX/F:cellsPerPlaneY/F");
349  ctist.append(":planesX/F:planesY/F:xStartPlane/F:yStartPlane/F:xEndPlane/F:yEndPlane/F");
350  ctist.append(":avgResidual/F:maxResidual/F:triCellFrac/F");
351  ctist.append(":stop/O:goodVtx/O:isCalibTrk/O");
352  fTree = fTFS->make<TTree>("Cosmics","CosmicTrackAna Tree");
353  fTree->Branch("trigger", &fTrigger, trigst.c_str());
354  fTree->Branch("truth", &fTrackTruth, ctist.c_str());
355  fTree->Branch("reco", &fTrackReco, ctist.c_str());
356  fTree->Branch("signalHitFraction", &fSigHitFrac, shfst.c_str());
357 
358  std::string chist ("pathLength/F:calibPathLength/F:truePathLength/F:triCellPathLength/F:PE/F:PECorr/F:w/F:distFromEnd/F:cell/I:plane/I:module/I:view/I");
359  fHitTree = fTFS->make<TTree>("Hits","CosmicHit Tree");
360  fHitTree->Branch("hits", &fCosmicHit, chist.c_str());
361 
362  }
363 
364  //......................................................................
365  void CosmicTrackAna::Make1DRecoHistograms(std::vector<TH1D *> &hists,
366  TString const sample)
367  {
368 
369  hists.resize(kMax1DRecoHist);
370 
371  double maxlength = 6600.0;
372  double xmax = 1200.0;
373  double ymax = 1200.0;
374  double zmax = 6500.0;
375  int xbins = 240;
376  int ybins = 240;
377  int zbins = 650;
378  double ncell = 384;
379  double nplane = 960;
380 
381  hists[kDCosX] = fTFS->make<TH1D>("dirCosX"+sample, ";Track direction cos_{X};Number of tracks", 101, -1., 1.);
382  hists[kDCosY] = fTFS->make<TH1D>("dirCosY"+sample, ";Track direction cos_{Y};Number of tracks", 101, -1., 1.);
383  hists[kDCosZ] = fTFS->make<TH1D>("dirCosZ"+sample, ";Track direction cos_{Z};Number of tracks", 101, -1., 1.);
384  hists[kEndX] = fTFS->make<TH1D>("endX"+sample, ";Track end X (cm);Number of tracks", xbins, -xmax, xmax);
385  hists[kEndY] = fTFS->make<TH1D>("endY"+sample, ";Track end Y (cm);Number of tracks", ybins, -ymax, ymax);
386  hists[kEndZ] = fTFS->make<TH1D>("endZ"+sample, ";Track end Z (cm);Number of tracks", zbins, 0., zmax);
387  hists[kStartX] = fTFS->make<TH1D>("startX"+sample, ";Track start X (cm);Number of tracks", xbins, -xmax, xmax);
388  hists[kStartY] = fTFS->make<TH1D>("startY"+sample, ";Track start Y (cm);Number of tracks", ybins, -ymax, ymax);
389  hists[kStartZ] = fTFS->make<TH1D>("startZ"+sample, ";Track start Z (cm);Number of tracks", zbins, 0., zmax);
390  hists[kDiffX] = fTFS->make<TH1D>("lenX"+sample, ";Track end minus start X (cm);Number of tracks", 100, 0., 1000.);
391  hists[kDiffY] = fTFS->make<TH1D>("lenY"+sample, ";Track end minus start Y (cm);Number of tracks", 100, 0., 1000.);
392  hists[kDiffZ] = fTFS->make<TH1D>("lenZ"+sample, ";Track end minus start Z (cm);Number of tracks", 150, 0., 1500.);
393  hists[kFractionHitsUsed] = fTFS->make<TH1D>("fractionHitsUsed"+sample, ";Used / total hits;Number of spills", 100, 0., 1.);
394  hists[kFractionSigUsed] = fTFS->make<TH1D>("fractionSigUsed"+sample, ";Used / total non-noise hits;Number of spills", 100, 0., 1.);
395  hists[kTracksFound] = fTFS->make<TH1D>("numTracks"+sample, ";Number of reconstructed tracks;", 25, 0, 125);
396  hists[kCellHitSig] = fTFS->make<TH1D>("cellHitSig"+sample, ";Hit PE;Number of CellHits", 100, 0., 100.);
397  hists[kTrkSig] = fTFS->make<TH1D>("trksig"+sample, ";Track Signal (Q);counts", 100,0.,100000.);
398  hists[kTrkCells] = fTFS->make<TH1D>("numHits"+sample, ";Number of CellHits in track;Number of tracks", 500,0,500);
399  hists[kEndDistToEdge] = fTFS->make<TH1D>("endDistToEdge"+sample, ";Track end distance to edge (cm);Number of tracks", 300, -1500., 1500.);
400  hists[kVtxDistToEdge] = fTFS->make<TH1D>("vtxDistToEdge"+sample, ";Track vertex distance to edge (cm);Number of tracks", 300, -1500., 1500.);
401  hists[kTrackLength] = fTFS->make<TH1D>("trackLength"+sample, ";Track Length (cm);Number of tracks", 264, 0., 2.*maxlength);
402  hists[kCosCosmic] = fTFS->make<TH1D>("coszenith"+sample, ";cos#theta_{zenith};Number of tracks", 100, 0., 1.);
403  hists[kAzimuth] = fTFS->make<TH1D>("azimuth"+sample, ";#phi (azimuth angle);Number of tracks", 36, 0., 360.);
404  hists[kHitXPosition] = fTFS->make<TH1D>("hitXPosition"+sample, ";Hit X (cm);Number of RecoHits", xbins, -xmax, xmax);
405  hists[kHitYPosition] = fTFS->make<TH1D>("hitYPosition"+sample, ";Hit Y (cm);Number of RecoHits", ybins, -ymax, ymax);
406  hists[kHitZPosition] = fTFS->make<TH1D>("hitZPosition"+sample, ";Hit Z (cm);Number of RecoHits", zbins, 0., zmax);
407  hists[kCellHitX] = fTFS->make<TH1D>("CellX"+sample, ";Track hit cell # in X view;Number of CellHits", ncell, -0.5, ncell-0.5);
408  hists[kCellHitY] = fTFS->make<TH1D>("CellY"+sample, ";Track hit cell # in Y view;Number of CellHits", ncell, -0.5, ncell-0.5);
409  hists[kPlaneHitX] = fTFS->make<TH1D>("PlaneX"+sample, ";Track hit plane # in X view;Number of CellHits", nplane, 0, nplane-0.5);
410  hists[kPlaneHitY] = fTFS->make<TH1D>("PlaneY"+sample, ";Track hit plane # in Y view;Number of CellHits", nplane, 0, nplane-0.5);
411  hists[kResidualsX] = fTFS->make<TH1D>("residualsX"+sample, ";Fit residual (cm) in X view;Number of hits", 400, 0., 1000.);
412  hists[kResidualsY] = fTFS->make<TH1D>("residualsY"+sample, ";Fit residual (cm) in Y view;Number of hits", 400, 0., 1000.);
413  hists[kTNSX] = fTFS->make<TH1D>("HitTNSX"+sample, ";Hit time (ns) in X view;Number of hits", 70000, -1.e5, 6.e5);
414  hists[kTNSY] = fTFS->make<TH1D>("HitTNSY"+sample, ";Hit time (ns) in Y view;Number of hits", 70000, -1.e5, 6.e5);
415  hists[kTrackRangeMomentum]= fTFS->make<TH1D>("rangeMomentum"+sample, ";p_{range} (GeV);Number of tracks", 100, 0., 100.);
416  hists[kCellsPerPlaneX] = fTFS->make<TH1D>("cellsPerPlaneX"+sample, ";<Cells/Plane>_{X view};Number of tracks", 200, 0., 50. );
417  hists[kCellsPerPlaneY] = fTFS->make<TH1D>("cellsPerPlaneY"+sample, ";<Cells/Plane>_{X view};Number of tracks", 200, 0., 50. );
418  hists[kStepRelativeToMedian] = fTFS->make<TH1D>("stepsRelativeToMedian"+sample, ";Step Size/Median Step;Number of tracks", 200, 0., 50.);
419  for( auto hist : hists) hist->Sumw2();
420 
421  return;
422  }
423 
424  //......................................................................
425  void CosmicTrackAna::Make1DTruthHistograms(std::vector<TH1D *> &hists,
426  TString const sample)
427  {
428 
429  hists.resize(kMax1DTruthHist);
430 
431  double maxlength = 6600.0;
432 
433  hists[kDCosXDiff] = fTFS->make<TH1D>("dcosXDiff"+sample, ";Fractional error in X direction cosine;Number of tracks", 200, -1., 3.);
434  hists[kDCosYDiff] = fTFS->make<TH1D>("dcosYDiff"+sample, ";Fractional error in Y direction cosine;Number of tracks", 200, -1., 3.);
435  hists[kDCosZDiff] = fTFS->make<TH1D>("dcosZDiff"+sample, ";Fractional error in Z direction cosine;Number of tracks", 200, -1., 3.);
436  hists[kPrimaryTrueEnergy] = fTFS->make<TH1D>("primaryTrueEnergy"+sample, ";True primary energy (GeV);Number of tracks", 200, 0., 50.);
437  hists[kRecoEff] = fTFS->make<TH1D>("recoeff"+sample, ";Fraction of tracks reconstructed;Number of spills", 100, 0., 1.);
438  hists[kTrackPurity] = fTFS->make<TH1D>("trackPurity"+sample, ";Tracking hit purity;Number of tracks", 100, 0., 2.);
439  hists[kTrackEfficiency] = fTFS->make<TH1D>("trackEfficiency"+sample, ";Tracking hit efficiency;Number of tracks", 100, 0., 2.);
440  hists[kDeltaTrackLength] = fTFS->make<TH1D>("errTrackLength"+sample, ";Error in track length (cm);Number of tracks", 264, -maxlength, maxlength);
441  hists[kDeltaRangeMomentum]= fTFS->make<TH1D>("errRangeMom"+sample, ";Error in #Delta p_{range} (GeV/c);Number of tracks", 400, -2, 2.);
442  hists[kDeltaRangeMomentumUncor]= fTFS->make<TH1D>("deltaRangeMomUncor"+sample, ";Error in #Delta p_{range} (uncorrected) (GeV/c);Number of tracks", 100, -2, 2.);
443  hists[kDeltaZEnd] = fTFS->make<TH1D>("errZEnd"+sample, ";Error in Z end position (cm);Number of tracks", 200, -100., 100.);
444  hists[kDeltaEnd] = fTFS->make<TH1D>("errEnd"+sample, ";Error in end position (cm);Number of tracks", 200, -100., 100.);
445  hists[kDeltaBeg] = fTFS->make<TH1D>("errStart"+sample, ";Error in start position (cm);Number of tracks", 200, -100., 100.);
446  hists[kEndTrueX] = fTFS->make<TH1D>("trueEndX"+sample, ";True end position X (cm);Number of true tracks", 80, -800., 800.);
447  hists[kEndTrueY] = fTFS->make<TH1D>("trueEndY"+sample, ";True end position Y (cm);Number of true tracks", 80, -800., 800.);
448  hists[kEndTrueZ] = fTFS->make<TH1D>("trueEndZ"+sample, ";True end position Z (cm);Number of true tracks", 300, 0., 6000.);
449  hists[kStartTrueX] = fTFS->make<TH1D>("trueStartX"+sample, ";True start position X (cm);Number of true tracks", 80, -800., 800.);
450  hists[kStartTrueY] = fTFS->make<TH1D>("trueStartY"+sample, ";True start position Y (cm);Number of true tracks", 80, -800., 800.);
451  hists[kStartTrueZ] = fTFS->make<TH1D>("trueStartZ"+sample, ";True start position Z (cm);Number of true tracks", 300, 0., 6000.);
452  hists[kCosCosmicTrue] = fTFS->make<TH1D>("coscosmic"+sample, ";True cos#theta;Number of true tracks", 100, 0., 1.);
453  hists[kAzimuthTrue] = fTFS->make<TH1D>("azimuth"+sample, ";True #phi;Number of true tracks", 36, 0., 360.);
454 
455  for( auto hist : hists) hist->Sumw2();
456 
457  return;
458  }
459 
460  //......................................................................
461  void CosmicTrackAna::Make2DRecoHistograms(std::vector<TH2D *> &hists,
462  TString const sample)
463  {
464 
465  hists.resize(kMax2DRecoHist);
466 
467  double xmax = 1200.0;
468  double ymax = 1200.0;
469  double zmax = 6500.0;
470  int xbins = 24;
471  int ybins = 24;
472  int zbins = 130;
473  double ncell = 384;
474  double nplane = 960;
475 
476  hists[kEndXY] = fTFS->make<TH2D>("endXY"+sample, ";End X (cm); End Y (cm)", xbins, -xmax, xmax, ybins, -ymax, ymax);
477  hists[kEndZY] = fTFS->make<TH2D>("endZY"+sample, ";End Z (cm); End Y (cm)", zbins, 0., zmax, ybins, -ymax, ymax);
478  hists[kEndZX] = fTFS->make<TH2D>("endZX"+sample, ";End Z (cm); End X (cm)", zbins, 0., zmax, xbins, -xmax, xmax);
479  hists[kStartXY] = fTFS->make<TH2D>("startXY"+sample, ";Start X (cm); Start Y (cm)", xbins, -xmax, xmax, ybins, -ymax, ymax);
480  hists[kStartZY] = fTFS->make<TH2D>("startZY"+sample, ";Start Z (cm); Start Y (cm)", zbins, 0., zmax, ybins, -ymax, ymax);
481  hists[kStartZX] = fTFS->make<TH2D>("startZX"+sample, ";Start Z (cm); Start X (cm)", zbins, 0., zmax, xbins, -xmax, xmax);
482  hists[kDirXY] = fTFS->make<TH2D>("dirXY"+sample, ";Prong Direction X;Prong Direction Y", 101, -1., 1., 101, -1., 1.);
483  hists[kFracHitsUsedVsNum] = fTFS->make<TH2D>("fractionHitsUsedVsNum"+sample, ";Total Hits; Fraction Used", 250, 0., 500., 100, 0., 1.);
484  hists[kHitXYPosition] = fTFS->make<TH2D>("hitXYPosition"+sample, ";Hit X (cm);Hit Y (cm)", 20, -200., 200., 50, -250., 250.);
485  hists[kHitXZPosition] = fTFS->make<TH2D>("hitXZPosition"+sample, ";Hit Z (cm);Hit X (cm)", 150, 0., 1500., 40, -200., 200.);
486  hists[kHitYZPosition] = fTFS->make<TH2D>("hitYZPosition"+sample, ";Hit Z (cm);Hit Y (cm)", 150, 0., 1500., 50, -250., 250.);
487  hists[kCellPlaneX] = fTFS->make<TH2D>("ViewX"+sample, ";Plane no. ;Cell no.", nplane, 0., nplane-0.5, ncell, 0., ncell-0.5);
488  hists[kCellPlaneY] = fTFS->make<TH2D>("ViewY"+sample, ";Plane no. ;Cell no.", nplane, 0., nplane-0.5, ncell, 0., ncell-0.5);
489  hists[kResidVsCosX] = fTFS->make<TH2D>("residVsCosX"+sample, ";Residual (cm);cos(#theta)", 200, 0., 1000., 100, -1., 1.);
490  hists[kResidVsCosY] = fTFS->make<TH2D>("residVsCosY"+sample, ";Residual (cm);cos(#theta)", 200, 0., 1000., 100, -1., 1.);
491  hists[kResidVsAzX] = fTFS->make<TH2D>("residVsAzX"+sample, ";Residual (cm);#phi", 200, 0., 1000., 180, 0., 360.);
492  hists[kResidVsAzY] = fTFS->make<TH2D>("residVsAzY"+sample, ";Residual (cm);#phi", 200, 0., 1000., 180, 0., 360.);
493  hists[kCellsPerPlaneXY] = fTFS->make<TH2D>("cellsPerPlaneXY"+sample, ";<Cells/Plane>_{X view};<Cells/Plane>_{Y view}", 200, 0., 50., 200, 0., 50.);
494  hists[kPEPerCMVsWX] = fTFS->make<TH2D>("pePerCMVsWX"+sample, ";W (cm);PE/cm", xbins, -xmax, xmax, 60, 0., 300.);
495  hists[kPEPerCMVsWY] = fTFS->make<TH2D>("pePerCMVsWY"+sample, ";W (cm);PE/cm", ybins, -ymax, ymax, 60, 0., 300.);
496  hists[kStepRelativeToMedianVsAz] = fTFS->make<TH2D>("stepRelativeToMedianVsAz"+sample, ";#phi;Largest Step/Median Step", 120, 0., 360., 100, 0., 10.);
497 
498  for( auto hist : hists) hist->Sumw2();
499 
500  return;
501  }
502 
503  //......................................................................
504  void CosmicTrackAna::Make2DTruthHistograms(std::vector<TH2D *> &hists,
505  TString const sample)
506  {
507 
508  hists.resize(kMax2DTruthHist);
509 
510  double maxlength = 6600.0;
511 
512  hists[kDCosXRecoVsTrue] = fTFS->make<TH2D>("dcosXRecoVsTrue"+sample, ";dcos X_{true};dcos X_{reco}", 250, -1., 1., 250, -1., 1.);
513  hists[kDCosYRecoVsTrue] = fTFS->make<TH2D>("dcosYRecoVsTrue"+sample, ";dcos Y_{true};dcos Y_{reco}", 250, -1., 1., 250, -1., 1.);
514  hists[kDCosZRecoVsTrue] = fTFS->make<TH2D>("dcosZRecoVsTrue"+sample, ";dcos Z_{true};dcos Z_{reco}", 250, -1., 1., 250, -1., 1.);
515  hists[kDCosXDiffVsZLen] = fTFS->make<TH2D>("dcosXDiffVsZLen"+sample, ";Z length (cm);dcos_{X} Fractional Difference", 1320, 0., maxlength, 200, -1., 5.);
516  hists[kDCosYDiffVsZLen] = fTFS->make<TH2D>("dcosYDiffVsZLen"+sample, ";Z length (cm);dcos_{Y} Fractional Difference", 1320, 0., maxlength, 200, -1., 3.);
517  hists[kDCosZDiffVsZLen] = fTFS->make<TH2D>("dcosZDiffVsZLen"+sample, ";Z length (cm);dcos_{Z} Fractional Difference", 1320, 0., maxlength, 200, -1., 3.);
518  hists[kTracksFoundVsLeft] = fTFS->make<TH2D>("tracksFoundVsLeft"+sample, ";True tracks;Reco tracks", 25, 0, 125, 25, 0, 125);
519  hists[kTrackFindEff] = fTFS->make<TH2D>("trackfindeff"+sample, ";True tracks;Reco Efficiency", 20, 0, 10, 100, 0, 1);
520  hists[kDeltaZEndVsDCosZ] = fTFS->make<TH2D>("deltaZEndVsDCosZ"+sample, ";#Delta Z_{end} (cm);dZ/dS;", 200, -100., 100., 20, -1., 1.);
521  hists[kDeltaZBegVsDCosZ] = fTFS->make<TH2D>("deltaZBegVsDCosZ"+sample, ";#Delta Z_{start} (cm);dZ/dS;", 200, -100., 100., 20, -1., 1.);
522  hists[kDiffEndVsDCosZ] = fTFS->make<TH2D>("diffEndVsDCosZ"+sample, ";dZ/dS;#Delta_{end} (cm)", 20, -1., 1., 200, -100., 100.);
523  hists[kDiffBegVsDCosZ] = fTFS->make<TH2D>("diffBegVsDCosZ"+sample, ";dZ/dS;#Delta_{start} (cm)", 20, -1., 1., 200, -100., 100.);
524  hists[kDiffEndVsAzimuth] = fTFS->make<TH2D>("diffEndVsAzimuth"+sample, ";#phi;#Delta_{end} (cm);", 120, 0., 360., 200, -100., 100.);
525  hists[kDiffBegVsAzimuth] = fTFS->make<TH2D>("diffBegVsAzimuth"+sample, ";#phi;#Delta_{start} (cm);", 120, 0., 360., 200, -100., 100.);
526  hists[kDiffEndVsCosTheta] = fTFS->make<TH2D>("diffEndVsCosTheta"+sample, ";cos(#theta);#Delta_{end} (cm);", 20, -1., 1., 200, -100., 100.);
527  hists[kDiffBegVsCosTheta] = fTFS->make<TH2D>("diffBegVsCosTheta"+sample, ";cos(#theta);#Delta_{start} (cm);", 20, -1., 1., 200, -100., 100.);
528  hists[kDeltaZEndVsDistToEdge] = fTFS->make<TH2D>("deltaZEndVsDistToEdge"+sample, ";#Delta Z_{end} (cm);Dist To Edge (cm);", 200, -100., 100., 200, -100., 100.);
529  hists[kDeltaRangeMomVsLength] = fTFS->make<TH2D>("deltaMomVsLength"+sample, ";l_{reco} (cm);#Deltap/p", 100, 0., 2500., 2000, -1., 1.);
530  hists[kDeltaRangeMomUncorVsLength] = fTFS->make<TH2D>("deltaMomUncorVsLength"+sample, ";l_{reco} (cm);#Deltap/p", 100, 0., 2500., 2000, -1., 1.);
531  hists[kDeltaRangeMomUncorVsCosTheta] = fTFS->make<TH2D>("deltaMomUncorVsCosTheta"+sample, ";cos#theta_{reco};#Deltap/p", 100, 0., 1., 2000, -1., 1.);
532  hists[kDeltaRangeMomUncorVsAzimuth] = fTFS->make<TH2D>("deltaMomUncorVsAzimuth"+sample, ";#phi_{reco} (deg);#Deltap/p", 100, 0., 360., 2000, -1., 1.);
533  hists[kDeltaLengthVsLength] = fTFS->make<TH2D>("deltaLengthVsLength"+sample, ";l_{reco} (cm);#Deltal/l", 100, 0., 2500., 2000, -1., 1.);
534  hists[kdEdxVsRange] = fTFS->make<TH2D>("dEdxVsRange"+sample, ";Range (cm);dE/dx (MeV)", 5000, 0., 5000., 2000, 0., 20.);
535  hists[kDeltaLengthVsPlaneAsymmetry] = fTFS->make<TH2D>("deltaLengthVsPlaneAsymmetry"+sample, ";Plane Asymmetry;#Deltal/l", 200, -1., 1., 2000, -1., 1.);
536  hists[kDeltaLengthVsDeltaStartPlane] = fTFS->make<TH2D>("deltaLengthVsDeltaStartPlane"+sample, ";#Delta Start Plane;#Deltal/l", 50, 0., 50., 2000, -1., 1.);
537  hists[kDeltaLengthVsDeltaEndPlane] = fTFS->make<TH2D>("deltaLengthVsDeltaEndPlane"+sample, ";#Delta Start Plane;#Deltal/l", 50, 0., 50., 2000, -1., 1.);
538  hists[kDeltaLengthVsCellsPerPlane] = fTFS->make<TH2D>("deltaLengthVsCellsPerPlane"+sample, ";Maximum Cells/Plane;#Deltal/l", 50, 0., 50., 2000, -1., 1.);
539 
540  for( auto hist : hists) hist->Sumw2();
541 
542  return;
543  }
544 
545  //......................................................................
547  {
549  sr.getByLabel(fExposureLabel, ce);
550  if(ce.failedToGet()){
551  mf::LogError("CosmicTrackAna") << "No Cosmic Exposure Info";
552  return;
553  }
554  fLiveTime += ce->totlivetime;
555 
556  MF_LOG_INFO("CosmicTrackAna") << "The live time is now " << fLiveTime;
557 
558  return;
559  }
560 
561  //......................................................................
562  ///
563  /// The analysis method is intended to perform checks of the
564  /// performance of the reconstruction for the purposes of
565  /// optimization, testing, and validation.
566  ///
567  /// \param evt - Read-only access to the event data
568  ///
569  // checks to make include -
570  // 1. are the direction cosines of the track correct
571  // 2. fraction of hits used in tracks
572  // 3. number of particles leaving tracks vs number tracks found
573  // 4. fraction of hits in a track that belong to a given particle
574  // 5. cosmic flux rate and angular distribtion stability
575  // 6. cosmic track efficiency
577  {
578  if(fDetLength < 0){
582  fCellHalfDepthZ = fGeom->Plane(0)->Cell(0)->HalfD();
583  }
584 
585  // make the histograms to record Truth related information
586  // if they havent already been made and if there is truth
587  // information in the event
588  if( fBT->HaveTruthInfo() && fOneDHistsTruthAll.size() < 1){
589  this->Make1DTruthHistograms(fOneDHistsTruthAll, "AllTruth");
590  this->Make1DTruthHistograms(fOneDHistsTruthStop, "StopTruth");
591  this->Make2DTruthHistograms(fTwoDHistsTruthAll, "AllTruth");
592  this->Make2DTruthHistograms(fTwoDHistsTruthStop, "StopTruth");
593  }
594 
595 
596  ++fNevents; //Updating counter of number of events in sub-run
597  if(fRun < 0) fRun = evt.run();
598  if(fSubRun < 0) fSubRun = evt.subRun();
599 
600  fTrigger.run = evt.run();
601  fTrigger.subrun = evt.subRun();
602  fTrigger.event = evt.id().event();
604 
605  // now get the tracks that were made by CosmicTrack
606  // and test that there are tracks in the event
608  evt.getByLabel(fTrackLabel,trks);
609  if (trks->size() < 1) {
610  mf::LogWarning("CTAWarn") << "There are " << trks->size()
611  << " Tracks in this event, skip";
612  return;
613  }
614  MF_LOG_DEBUG("CosmicTrackAna") << "There are " << trks->size() << " Tracks in this event";
615 
616  const unsigned int trkMax = trks->size();
617 
618  if( fBT->HaveTruthInfo() )
619  this->FillEventMCTruthHistograms(trkMax);
620 
621  // get a FindMany to grab the slice associated with each track
622  art::FindOne<rb::Cluster> foc(trks, evt, fTrackLabel);
623 
624  unsigned int numtracks = 0;
625  unsigned int numstop = 0;
626 
627  // TStopwatch ts;
628  // double mcTime = 0.;
629  // double recoTime = 0.;
630 
631  // Loop over tracks
632  for(unsigned int trkIdx = 0; trkIdx < trkMax; ++trkIdx){
633 
634  // find the slice associated to the track
635  cet::maybe_ref<rb::Cluster const> sliceRef(foc.at(trkIdx));
636  if( !sliceRef.isValid() ) continue;
637  rb::Cluster const& slice(sliceRef.ref());
638 
639  rb::Track const& track = trks->at(trkIdx);
640  if(!track.Is3D()) continue;
641 
642  MF_LOG_DEBUG("CosmicTrackAna") << "-----------------------------------------Track starts at "
643  << track.Start().X() << " "
644  << track.Start().Y() << " "
645  << track.Start().Z() << "---------------------------------";
646 
647  // fill the histogram for the step sizes compared to median for a track
648  if(std::abs(track.Stop().Z() - track.Start().Z()) >= 50. &&
649  -track.Dir().Y() < 0.95*track.Dir().Mag() ){
650 
651  std::vector<float> steps;
652  float median = -999.;
653  float largestStep = -999.;
654  fSelection.GoodSteps(track, steps, median);
655 
656  if(median > 0.){
657  for(auto const& step : steps){
659  largestStep = std::max(largestStep, step/median);
660  }
661  fTwoDHistsRecoAll[kStepRelativeToMedianVsAz]->Fill(fUtilities.Azimuth(track.Dir().X(), track.Dir().Z(), track.Dir().Mag()),
662  largestStep);
663  }
664  }// end filling of step size histogram
665 
666  // check the quality of the reconstruction
667  if(fSelection.GoodReco(track) ){
668 
669  //count number of tracks that cross fTrkLengthCut
670  ++numtracks;
671 
672  std::pair<CosmicTrackInfo,
673  CosmicSignalHitFraction> cticshf = this->FillRecoInfo(track, slice);
674 
675  fTrackReco = cticshf.first;
676  fSigHitFrac = cticshf.second;
677 
678  std::map<geo::CellUniqueId, float> trueCellPathLength;
679  fTrackTruth = this->FillTrueInfo(track, slice, trueCellPathLength);
680 
681  if(fTrackReco.stop) ++numstop;
682 
683  fTree->Fill();
684 
685  // ts.Start();
686  if( fBT->HaveTruthInfo() && std::abs(fTrackTruth.pdg) == 13 )
688 
689  // ts.Stop();
690  // mcTime += ts.RealTime();
691  // ts.Reset();
692 
693  MF_LOG_DEBUG("CosmicTrackAna") << "----------------Start find track path length in cell--------------";
694  this->FillTrackHistograms(fTrackReco, fSigHitFrac, track, trueCellPathLength);
695  MF_LOG_DEBUG("CosmicTrackAna") << "----------------End find track path length in cell--------------";
696 
697  }// end if passes track length in z cut
698 
699  }//end loop over tracks
700 
701  // MF_LOG_VERBATIM("CosmicTrackAna") << "mcTime: " << mcTime
702  // << " recoTime: " << recoTime;
703 
704 
705  fOneDHistsRecoAll[kTracksFound] ->Fill(numtracks);
706  fOneDHistsRecoStop[kTracksFound]->Fill(numstop);
707 
708  fNTracks +=numtracks;
709 
710  fTrigger.ntracks = numtracks;
711 
712  return;
713  }
714 
715  //......................................................................
717  {
718 
719  unsigned int numtracksMC = 0;
720 
721  std::vector<sim::Particle> detParticles;
722 
724 
725  // count how many tracks were left by primary particles in the event
726  for(int p = 0; p < pnav.NumberOfPrimaries(); ++p){
727  if( fBT->ParticleToFLSHit(pnav.Primary(p)->TrackId()).size() > 0){
728  if(std::abs(pnav.Primary(p)->PdgCode()) == 13){
729  ++fNTracksMC;
730  ++numtracksMC;
731  detParticles.push_back(*(pnav.Primary(p)));
732  }
733  }
734  }
735 
736  // no point filling histograms for MC if particles don't leave hits
737  if (detParticles.size() < 1 ) {
738  mf::LogWarning("CTAWarn") << "Failed to find any Mother particles leaving"
739  << " hits in the detector, bail";
740  return;
741  }
742  else {
743  MF_LOG_DEBUG("CosmicTrackAna") << "There are " << detParticles.size()
744  << " Mother particles leaving hits in the detector";
745 
746  for(auto detPart : detParticles){
747  auto lengthEndPoints = fUtilities.TrueLengthEndPoints(detPart);
748 
749  fOneDHistsTruthAll[kPrimaryTrueEnergy]->Fill(detPart.E());
750  fOneDHistsTruthAll[kAzimuthTrue] ->Fill(fUtilities.Azimuth(detPart.Px(), detPart.Pz(), detPart.P()));
751  fOneDHistsTruthAll[kCosCosmicTrue] ->Fill(fUtilities.CosTheta(detPart.Py(), detPart.P()));
752  fOneDHistsTruthAll[kStartTrueX] ->Fill(lengthEndPoints.second.first.X());
753  fOneDHistsTruthAll[kStartTrueY] ->Fill(lengthEndPoints.second.first.Y());
754  fOneDHistsTruthAll[kStartTrueZ] ->Fill(lengthEndPoints.second.first.Z());
755  fOneDHistsTruthAll[kEndTrueX] ->Fill(lengthEndPoints.second.second.X());
756  fOneDHistsTruthAll[kEndTrueY] ->Fill(lengthEndPoints.second.second.Y());
757  fOneDHistsTruthAll[kEndTrueZ] ->Fill(lengthEndPoints.second.second.Z());
758  if( fSelection.Stopper( detPart.EndPosition().Vect() ) ){
759  fOneDHistsTruthStop[kPrimaryTrueEnergy]->Fill(detPart.E());
760  fOneDHistsTruthStop[kAzimuthTrue] ->Fill(fUtilities.Azimuth(detPart.Px(), detPart.Pz(), detPart.P()));
761  fOneDHistsTruthStop[kCosCosmicTrue] ->Fill(fUtilities.CosTheta(detPart.Py(), detPart.P()));
762  fOneDHistsTruthStop[kStartTrueX] ->Fill(lengthEndPoints.second.first.X());
763  fOneDHistsTruthStop[kStartTrueY] ->Fill(lengthEndPoints.second.first.Y());
764  fOneDHistsTruthStop[kStartTrueZ] ->Fill(lengthEndPoints.second.first.Z());
765  fOneDHistsTruthStop[kEndTrueX] ->Fill(lengthEndPoints.second.second.X());
766  fOneDHistsTruthStop[kEndTrueY] ->Fill(lengthEndPoints.second.second.Y());
767  fOneDHistsTruthStop[kEndTrueZ] ->Fill(lengthEndPoints.second.second.Z());
768  }
769  }// end loop over particles (muons)
770 
771  double ratio = 0.;
772  fTwoDHistsTruthAll[kTracksFoundVsLeft]->Fill(numtracksMC*1., trkMax*1.); // fill only for MC
773 
774  if(trkMax > numtracksMC)
775  MF_LOG_VERBATIM("CosmicTrackAna") << "reconstructed more tracks than generated: "
776  << trkMax << " vs " << numtracksMC;
777 
778  ratio = (trkMax*1.)/(numtracksMC*1.);
779  fTwoDHistsTruthAll[kTrackFindEff]->Fill(numtracksMC*1.,ratio);
780  fOneDHistsTruthAll[kRecoEff]->Fill(ratio);
781 
782  }// end if particles were found
783 
784  return;
785  }
786 
787  //......................................................................
789  rb::Cluster const& slice,
790  std::map<geo::CellUniqueId, float>& trueCellPathLength)
791  {
792 
793  CosmicTrackInfo ti;
794 
795  ti.purity = -999.;
796  ti.efficiency = -999.;
797  ti.length = -999.;
798  ti.momentum = -999.;
799  ti.momentumCor = -999.;
800  ti.dcosX = -999.;
801  ti.dcosY = -999.;
802  ti.dcosZ = -999.;
803  ti.vtxX = -999.;
804  ti.vtxY = -999.;
805  ti.vtxZ = -999.;
806  ti.endX = -999.;
807  ti.endY = -999.;
808  ti.endZ = -999.;
809  ti.deltaVtxPos = -999.;
810  ti.deltaEndPos = -999.;
811  ti.cosTheta = -999.;
812  ti.azimuth = -999.;
813  ti.stop = false;
814  ti.goodVtx = false;
815  ti.vtxToEdge = -999.;
816  ti.endToEdge = -999.;
817  ti.cellsPerPlaneX = -999.;
818  ti.cellsPerPlaneY = -999.;
819 
820  if( !fBT->HaveTruthInfo() ) return ti;
821 
822  // get the truth information for this track
823  std::vector< const sim::Particle* > parts = fBT->HitsToParticle(track.AllCells());
824 
825  if(parts.size() > 0){
826 
827  // get the particle for the first IDE as it represents the track ID with the most
828  // energy in the collection of hits
829  auto part = parts.front();
830 
831  // if(track.Dir().Z() > 0)
832  // MF_LOG_VERBATIM("CosmicTrackAna") << "simulation z end point: " << simEnd.Z()
833  // << "\ntrack z end point: " << track.Stop().Z();
834 
835  this->FilldEdxHistogram(*part);
836  auto lengthEndPoints = fUtilities.TrueLengthEndPoints(*part);
837 
838  // get the FLSHits for this track, then make a map
839  // of the cell unique id to the path length of the particle in the cell
840  std::vector<sim::FLSHit> flsHits = fBT->ParticleToFLSHit(part->TrackId());
841  // for( auto flsHit : flsHits)
842  // trueCellPathLength[flsHit.GetCellUniqueId()] = flsHit.GetTotalPathLength();
843  MF_LOG_DEBUG("CosmicTrackAna") << "------------Fill true cell path lengths-------------";
844  trueCellPathLength = fUtilities.TruePathLengthInCell(*part, flsHits);
845  MF_LOG_DEBUG("CosmicTrackAna") << "------------End Fill true cell path lengths-------------";
846 
847  // calculating purity and efficiency take forever.
848  // control whether it is done with a flag
850  // check the efficiency and purity of the track hits for
851  // this MC particle
852  std::set<int> trackID;
853  trackID.insert(part->TrackId());
854  ti.purity = fBT->HitCollectionPurity(trackID, track.AllCells());
855  ti.efficiency = fBT->HitCollectionEfficiency(trackID, track.AllCells(), slice.AllCells(), geo::kXorY);
856  }
857 
858  ti.pdg = part->PdgCode();
859  ti.length = lengthEndPoints.first.first;
860  ti.momentum = lengthEndPoints.first.second;
861  ti.momentumCor = 0.;
862  ti.dcosX = part->Px()/part->P();
863  ti.dcosY = part->Py()/part->P();
864  ti.dcosZ = part->Pz()/part->P();
865  ti.vtxX = lengthEndPoints.second.first.X();
866  ti.vtxY = lengthEndPoints.second.first.Y();
867  ti.vtxZ = lengthEndPoints.second.first.Z();
868  ti.endX = lengthEndPoints.second.second.X();
869  ti.endY = lengthEndPoints.second.second.Y();
870  ti.endZ = lengthEndPoints.second.second.Z();
871  ti.cosTheta = fUtilities.CosTheta(part->Py(), part->P());
872  ti.azimuth = fUtilities.Azimuth(part->Px(), part->Pz(), part->P());
873  ti.stop = fSelection.Stopper(lengthEndPoints.second.second);
874  ti.goodVtx = (std::abs(fSelection.DistToEdge(track.Start())) < 10.) ? true : false;
875  ti.vtxToEdge = fSelection.DistToEdge(lengthEndPoints.second.first);
876  ti.endToEdge = fSelection.DistToEdge(lengthEndPoints.second.second);
877  ti.cellsPerPlaneX = 0.;
878  ti.cellsPerPlaneY = 0.;
879 
880  // fill information about the difference between the true and reconstructed end points
881  TVector3 simEnd(ti.endX, ti.endY, ti.endZ);
882  TVector3 simBeg(ti.vtxX, ti.vtxY, ti.vtxZ);
883  TVector3 recEnd(track.Stop().X(), track.Stop().Y(), track.Stop().Z());
884  TVector3 recBeg(track.Start().X(), track.Start().Y(), track.Start().Z());
885  TVector3 diffEnd = recEnd - simEnd;
886  TVector3 diffBeg = recBeg - simBeg;
887 
888  float diffEndSign = 1.;
889  float recoDCosZ = track.Dir().Z()/track.Dir().Mag();
890  if( (recoDCosZ > 0. && recEnd.Z() > simEnd.Z()) ||
891  (recoDCosZ < 0. && recEnd.Z() < simEnd.Z()) ) diffEndSign = -1.;
892 
893  float diffBegSign = 1.;
894  if( (recoDCosZ > 0. && recBeg.Z() < simBeg.Z()) ||
895  (recoDCosZ < 0. && recBeg.Z() > simBeg.Z()) ) diffBegSign = -1.;
896 
897  ti.deltaVtxPos = diffBegSign*diffBeg.Mag();
898  ti.deltaEndPos = diffEndSign*diffEnd.Mag();
899 
900  }
901 
902  return ti;
903 
904  }
905 
906  //......................................................................
908  CosmicTrackInfo const& ri)
909  {
910 
911 
912  float rangeMomUncor = ri.momentum/fUtilities.LengthCorrection(ri.length);
913  float endDiffSign = ti.deltaEndPos/std::abs(ti.deltaEndPos);
914  float begDiffSign = ti.deltaVtxPos/std::abs(ti.deltaVtxPos);
915 
916  fOneDHistsTruthAll[kDeltaZEnd] ->Fill(endDiffSign*std::abs(ri.endZ - ti.endZ));
919  fOneDHistsTruthAll[kDCosXDiff] ->Fill((ti.dcosX - ri.dcosX)/ti.dcosX);
920  fOneDHistsTruthAll[kDCosYDiff] ->Fill((ti.dcosY - ri.dcosY)/ti.dcosX);
921  fOneDHistsTruthAll[kDCosZDiff] ->Fill((ti.dcosZ - ri.dcosZ)/ti.dcosX);
930  fTwoDHistsTruthAll[kDCosXDiffVsZLen] ->Fill(std::abs(ri.vtxZ - ri.endZ), (ti.dcosX - ri.dcosZ)/ti.dcosX);
931  fTwoDHistsTruthAll[kDCosYDiffVsZLen] ->Fill(std::abs(ri.vtxZ - ri.endZ), (ti.dcosY - ri.dcosY)/ti.dcosY);
932  fTwoDHistsTruthAll[kDCosZDiffVsZLen] ->Fill(std::abs(ri.vtxZ - ri.endZ), (ti.dcosZ - ri.dcosZ)/ti.dcosZ);
933  fTwoDHistsTruthAll[kDeltaZEndVsDCosZ] ->Fill(endDiffSign * std::abs(ri.endZ - ti.endZ), ri.dcosZ);
934  fTwoDHistsTruthAll[kDeltaZEndVsDistToEdge] ->Fill(endDiffSign * std::abs(ri.endZ - ti.endZ), ti.endToEdge);
935  fTwoDHistsTruthAll[kDeltaZBegVsDCosZ] ->Fill(begDiffSign * std::abs(ri.vtxZ - ti.vtxZ), ri.dcosZ);
949 
950  // check and see how well the reconstruction works for stoppers selected based on the reconstruction
951  if( fSelection.CleanSample(ri) ){
952  fOneDHistsTruthStop[kDCosXDiff] ->Fill((ti.dcosX - ri.dcosX)/ti.dcosX);
953  fOneDHistsTruthStop[kDCosYDiff] ->Fill((ti.dcosY - ri.dcosY)/ti.dcosY);
954  fOneDHistsTruthStop[kDCosZDiff] ->Fill((ti.dcosZ - ri.dcosZ)/ti.dcosZ);
955  fOneDHistsTruthStop[kDeltaZEnd] ->Fill(endDiffSign*std::abs(ri.endZ - ti.endZ));
966  fTwoDHistsTruthStop[kDCosXDiffVsZLen] ->Fill(std::abs(ri.vtxZ - ri.endZ), (ti.dcosX - ri.dcosX)/ti.dcosX);
967  fTwoDHistsTruthStop[kDCosYDiffVsZLen] ->Fill(std::abs(ri.vtxZ - ri.endZ), (ti.dcosY - ri.dcosY)/ti.dcosY);
968  fTwoDHistsTruthStop[kDCosZDiffVsZLen] ->Fill(std::abs(ri.vtxZ - ri.endZ), (ti.dcosZ - ri.dcosZ)/ti.dcosZ);
969  fTwoDHistsTruthStop[kDeltaZEndVsDCosZ] ->Fill(endDiffSign*std::abs(ri.endZ - ti.endZ), ri.dcosZ);
970  fTwoDHistsTruthStop[kDeltaZEndVsDistToEdge] ->Fill(endDiffSign*std::abs(ri.endZ - ti.endZ), ri.endToEdge);
971  fTwoDHistsTruthStop[kDeltaZBegVsDCosZ] ->Fill(begDiffSign*std::abs(ri.vtxZ - ti.vtxZ), ri.dcosZ);
987  }
988 
989  return;
990  }
991 
992  //......................................................................
994  {
995 
996  // loop over the particle trajectory points and determine the
997  // energy loss for each set
998  TVector3 curPos, prevPos;
999  float totalLength = 0.;
1000 
1001  std::vector< std::pair<float, float> > dEdXVsTotalL;
1002 
1003  for(size_t t = 1; t < part.NumberTrajectoryPoints(); ++t){
1004 
1005  if( std::abs(part.Position(t).X()) < fDetHalfWidth &&
1006  std::abs(part.Position(t).Y()) < fDetHalfHeight &&
1007  part.Position(t).Z() > 0. && part.Position(t).Z() < fDetLength &&
1008  std::abs(part.Position(t-1).X()) < fDetHalfWidth &&
1009  std::abs(part.Position(t-1).Y()) < fDetHalfHeight &&
1010  part.Position(t-1).Z() > 0. && part.Position(t-1).Z() < fDetLength ){
1011 
1012  curPos.SetXYZ(part.Position(t).X(), part.Position(t).Y(), part.Position(t).Z());
1013  prevPos.SetXYZ(part.Position(t-1).X(), part.Position(t-1).Y(), part.Position(t-1).Z());
1014 
1015  float dx = (curPos - prevPos).Mag();
1016  float dE = 1.e3*(part.E(t-1) - part.E(t)); // in MeV
1017 
1018  totalLength += dx;
1019 
1020  if(dE < 0. && std::abs(dE) > 1.e-3)
1021  MF_LOG_WARNING("CosmicTrackAna") << "MC dE is negative!! " << dE
1022  << " " << dx
1023  << " " << part.E(t) << " " << part.E(t-1);
1024  if(dx <= 0.){
1025  MF_LOG_WARNING("CosmicTrackAna") << "MC indicates a bad dx: " << dx;
1026  continue;
1027  }
1028 
1029  dEdXVsTotalL.push_back(std::pair<float, float>(dE/dx, totalLength));
1030  // MF_LOG_VERBATIM("CosmicTrackAna") << dE/dx << " " << totalLength;
1031 
1032  }
1033  }
1034 
1035  // loop over the vector of pairs to plot dE/dx vs range
1036  for(auto const& dEdxTotalL : dEdXVsTotalL){
1037  fTwoDHistsTruthAll[kdEdxVsRange]->Fill(totalLength - dEdxTotalL.second, dEdxTotalL.first);
1038  if( fSelection.Stopper(curPos) )
1039  fTwoDHistsTruthStop[kdEdxVsRange]->Fill(totalLength - dEdxTotalL.second, dEdxTotalL.first);
1040  }
1041 
1042  return;
1043 
1044  }
1045 
1046  //......................................................................
1047  std::pair<CosmicTrackInfo, CosmicSignalHitFraction> CosmicTrackAna::FillRecoInfo(rb::Track const& track,
1048  rb::Cluster const& slice)
1049  {
1050  // determine the number of planes with hits in each view
1051  std::set<unsigned short> xPlanes;
1052  std::set<unsigned short> yPlanes;
1053  for(size_t c = 0; c < track.NCell(); ++c){
1054  if (track.Cell(c)->View() == geo::kX) xPlanes.insert(track.Cell(c)->Plane());
1055  else if(track.Cell(c)->View() == geo::kY) yPlanes.insert(track.Cell(c)->Plane());
1056  }// end loop over cells in track
1057 
1058  // get the residuals along the track
1059  std::vector<float> residuals = fUtilities.TrackResiduals(track);
1060  std::sort(residuals.begin(), residuals.end());
1061  float avResid = 0.;
1062  for(auto resid : residuals) avResid += resid;
1063 
1064  // find the number of tricells in this track
1065  std::map<std::pair<uint32_t, uint32_t>, bool> planeCellMap;
1066  fUtilities.FindTriCells(track, planeCellMap);
1067  float numTriCell = 0.;
1068  for(auto itr : planeCellMap)
1069  if(itr.second) numTriCell += 1.;
1070 
1071  // find the momentum from the BreakPointFitter
1072  bpfit::Path path(*fGeom, track.Trajectory());
1073 
1074  double bpfMom = 0.;
1075  path.MuonParams(0, 0, &bpfMom, 0, 0);
1076 
1077  CosmicTrackInfo ri;
1078  ri.planesX = xPlanes.size();
1079  ri.planesY = yPlanes.size();
1080  ri.length = track.TotalLength();
1081  ri.momentum = fUtilities.RangeMomentum(track);
1083  ri.momentumBPF = bpfMom;
1084  ri.calorE = track.TotalGeV();
1085  ri.vtxX = track.Start().X();
1086  ri.vtxY = track.Start().Y();
1087  ri.vtxZ = track.Start().Z();
1088  ri.endX = track.Stop().X();
1089  ri.endY = track.Stop().Y();
1090  ri.endZ = track.Stop().Z();
1091  ri.dcosX = track.Dir().X()/track.Dir().Mag();
1092  ri.dcosY = track.Dir().Y()/track.Dir().Mag();
1093  ri.dcosZ = track.Dir().Z()/track.Dir().Mag();
1094  ri.cosTheta = fUtilities.CosTheta(track.Dir().Y(), track.Dir().Mag());
1095  ri.azimuth = fUtilities.Azimuth(track.Dir().X(), track.Dir().Z(), track.Dir().Mag());
1096  ri.vtxToEdge = fSelection.DistToEdge(track.Start());
1097  ri.endToEdge = fSelection.DistToEdge(track.Stop());
1098  ri.xStartPlane = 1. * track.MaxPlane(geo::kX);
1099  ri.yStartPlane = 1. * track.MaxPlane(geo::kY);
1100  ri.xEndPlane = 1. * track.MinPlane(geo::kX);
1101  ri.yEndPlane = 1. * track.MinPlane(geo::kY);
1102  ri.cellsPerPlaneX = (ri.planesX > 0.) ? (1. * track.NXCell())/ri.planesX : -999.;
1103  ri.cellsPerPlaneY = (ri.planesY > 0.) ? (1. * track.NYCell())/ri.planesY : -999.;
1104  ri.avgResidual = avResid/track.NCell();
1105  ri.maxResidual = residuals.back();
1106  ri.triCellFrac = 1.*numTriCell/(1.*track.NCell());
1107  ri.stop = fSelection.Stopper(track.Stop());
1108  ri.goodVtx = (std::abs(fSelection.DistToEdge(track.Start())) < 10.) ? true : false;
1109  ri.isCalibTrk = fSelection.IsCalibTrack(track, slice);
1110 
1111  float trkcells = 0.;
1112  float trksig = 0.;
1113  float totcell = 0.;
1114  float totsig = 0.;
1115  float xsig = 0.;
1116  float ysig = 0.;
1117 
1118  // loop over the cell hits from the cluster
1119  for(size_t i = 0; i < slice.NCell(); ++i){
1120  art::Ptr<rb::CellHit> chit = slice.Cell(i);
1121  if(chit->PE() > fMinSig){
1122  fOneDHistsRecoAll[kCellHitSig]->Fill(chit->PE());
1124  totcell += 1.;
1125  totsig += chit->PE();
1126  }
1127  }
1128 
1129  // get the total signal in each view of the track
1130  for(size_t c = 0; c < track.NCell(); ++c){
1131  if(track.Cell(c)->PE() > fMinSig){
1132  trkcells += 1.;
1133  if (track.Cell(c)->View() == geo::kX) xsig += track.Cell(c)->PE()*track.Weight(c);
1134  else if(track.Cell(c)->View() == geo::kY) ysig += track.Cell(c)->PE()*track.Weight(c);
1135  //MF_LOG_VERBATIM("CosmicTrackAna") << xsig << " " << ysig << " " << " " << track.Cell(c)->PE() << " " << track.Weight(c);
1136  }
1137  }
1138 
1139  trksig = xsig + ysig;
1140 
1142  shf.signalFraction = trksig/totsig;
1143  shf.hitFraction = trkcells/totcell;
1144  shf.signal = trksig;
1145  shf.hitsX = track.NXCell();
1146  shf.hitsY = track.NYCell();
1147  shf.totcell = totcell;
1148  shf.signalX = xsig;
1149  shf.signalY = ysig;
1150 
1151  return std::pair<CosmicTrackInfo, CosmicSignalHitFraction>(ri, shf);
1152  }
1153 
1154  //......................................................................
1156  CosmicSignalHitFraction const& shf,
1157  rb::Track const& track,
1158  std::map<geo::CellUniqueId, float> const& trueCellPathLength)
1159  {
1160 
1161  fOneDHistsRecoAll[kTrkSig] ->Fill(shf.signal);
1162  fOneDHistsRecoAll[kTrkCells] ->Fill(shf.hitsX+shf.hitsY);
1163  fOneDHistsRecoAll[kDCosX] ->Fill(ri.dcosX);
1164  fOneDHistsRecoAll[kDCosY] ->Fill(ri.dcosY);
1165  fOneDHistsRecoAll[kDCosZ] ->Fill(ri.dcosZ);
1166  fOneDHistsRecoAll[kAzimuth] ->Fill(ri.azimuth);
1168  fOneDHistsRecoAll[kStartX] ->Fill(ri.vtxX);
1169  fOneDHistsRecoAll[kStartY] ->Fill(ri.vtxY);
1170  fOneDHistsRecoAll[kStartZ] ->Fill(ri.vtxZ);
1171  fOneDHistsRecoAll[kEndX] ->Fill(ri.endX);
1172  fOneDHistsRecoAll[kEndY] ->Fill(ri.endY);
1173  fOneDHistsRecoAll[kEndZ] ->Fill(ri.endZ);
1176  fOneDHistsRecoAll[kDiffX] ->Fill(std::abs(ri.vtxX - ri.endX));
1177  fOneDHistsRecoAll[kDiffY] ->Fill(std::abs(ri.vtxY - ri.endY));
1178  fOneDHistsRecoAll[kDiffZ] ->Fill(std::abs(ri.vtxZ - ri.endZ));
1185  fTwoDHistsRecoAll[kDirXY] ->Fill(ri.dcosX,ri.dcosY);
1186  fTwoDHistsRecoAll[kStartXY] ->Fill(ri.vtxX, ri.vtxY);
1187  fTwoDHistsRecoAll[kStartZX] ->Fill(ri.vtxZ, ri.vtxX);
1188  fTwoDHistsRecoAll[kStartZY] ->Fill(ri.vtxZ, ri.vtxY);
1189  fTwoDHistsRecoAll[kEndXY] ->Fill(ri.endX, ri.endY);
1190  fTwoDHistsRecoAll[kEndZX] ->Fill(ri.endZ, ri.endX);
1191  fTwoDHistsRecoAll[kEndZY] ->Fill(ri.endZ, ri.endY);
1194 
1195  if( fSelection.CleanSample(ri) ){
1196  fOneDHistsRecoStop[kTrkSig] ->Fill(shf.signal);
1197  fOneDHistsRecoStop[kTrkCells] ->Fill(shf.hitsX+shf.hitsY);
1198  fOneDHistsRecoStop[kDCosX] ->Fill(ri.dcosX);
1199  fOneDHistsRecoStop[kDCosY] ->Fill(ri.dcosY);
1200  fOneDHistsRecoStop[kDCosZ] ->Fill(ri.dcosZ);
1201  fOneDHistsRecoStop[kAzimuth] ->Fill(ri.azimuth);
1202  fOneDHistsRecoStop[kCosCosmic] ->Fill(ri.cosTheta);
1203  fOneDHistsRecoStop[kStartX] ->Fill(ri.vtxX);
1204  fOneDHistsRecoStop[kStartY] ->Fill(ri.vtxY);
1205  fOneDHistsRecoStop[kStartZ] ->Fill(ri.vtxZ);
1206  fOneDHistsRecoStop[kEndX] ->Fill(ri.endX);
1207  fOneDHistsRecoStop[kEndY] ->Fill(ri.endY);
1208  fOneDHistsRecoStop[kEndZ] ->Fill(ri.endZ);
1209  fOneDHistsRecoStop[kTrackLength] ->Fill(ri.length);
1211  fOneDHistsRecoStop[kDiffX] ->Fill(std::abs(ri.vtxX - ri.endX));
1212  fOneDHistsRecoStop[kDiffY] ->Fill(std::abs(ri.vtxY - ri.endY));
1213  fOneDHistsRecoStop[kDiffZ] ->Fill(std::abs(ri.vtxZ - ri.endZ));
1220  fTwoDHistsRecoStop[kDirXY] ->Fill(ri.dcosX,ri.dcosY);
1221  fTwoDHistsRecoStop[kStartXY] ->Fill(ri.vtxX, ri.vtxY);
1222  fTwoDHistsRecoStop[kStartZX] ->Fill(ri.vtxZ, ri.vtxX);
1223  fTwoDHistsRecoStop[kStartZY] ->Fill(ri.vtxZ, ri.vtxY);
1224  fTwoDHistsRecoStop[kEndXY] ->Fill(ri.endX, ri.endY);
1225  fTwoDHistsRecoStop[kEndZX] ->Fill(ri.endZ, ri.endX);
1226  fTwoDHistsRecoStop[kEndZY] ->Fill(ri.endZ, ri.endY);
1229  }
1230 
1231  // make a map of the z position of each trajectory point of the
1232  // track to its corresponding trajectory point. the map will
1233  // be used to find the residual of the reconstruction to the
1234  // cell hits.
1235  // TVector3 cellCenter;
1236  // TVector3 cellExtent;
1237  // TVector3 lowerBound;
1238  // TVector3 upperBound;
1239  TVector3 point;
1240  float distFromEnd = 0.;
1241  float lengthInCell = 0.;
1242  float calibLengthInCell = 0.;
1243  float trueLengthInCell = 0.;
1244  // double xyz[3] = {0.};
1245 
1246  auto trueCellItr = trueCellPathLength.begin();
1247  zBoundMap zToPlaneBounds = fUtilities.MakeZBoundaryMap(track.Trajectory());
1248  std::set<double> calibZBoundaries;
1249  calib::FindZBoundaries(calibZBoundaries);
1250  auto calibZToPlaneBounds = calib::MakeZBoundaryMap(calibZBoundaries,
1251  track.Trajectory());
1252 
1253  std::map< std::pair<uint32_t, uint32_t>, bool > triCellMap;
1254  std::pair<uint32_t, uint32_t> pc;
1255  double cellWidth = 0.;
1256  fUtilities.FindTriCells(track, triCellMap);
1257 
1258  auto residuals = fUtilities.TrackResiduals(track);
1259 
1260  // Loop through both views and grab every CellHit in this track
1261  for(size_t c = 0; c < track.NCell(); ++c){
1262  const art::Ptr<rb::CellHit> chit = track.Cell(c);
1263  const geo::View_t view = chit->View();
1264 
1265  pc.first = chit->Plane();
1266  pc.second = chit->Cell();
1267  cellWidth = 2. * fGeom->Plane(pc.first)->Cell(pc.second)->HalfW();
1268 
1269  trueCellItr = trueCellPathLength.find(fGeom->Plane(pc.first)->Cell(pc.second)->Id());
1270  if( trueCellItr != trueCellPathLength.end() ) trueLengthInCell = trueCellItr->second;
1271  else trueLengthInCell = -9999.;
1272 
1273  // MF_LOG_VERBATIM("CosmicTrackAna") << chit->Plane() << " " << chit->Cell() << "cell hit view is " << view;
1274 
1275  if(view == geo::kX){
1276  fOneDHistsRecoAll[kCellHitX] ->Fill(pc.second);
1277  fOneDHistsRecoAll[kPlaneHitX] ->Fill(pc.first);
1278  fTwoDHistsRecoAll[kCellPlaneX]->Fill(pc.first, pc.second);
1279  if( fSelection.CleanSample(ri) ){
1280  fOneDHistsRecoStop[kCellHitX] ->Fill(pc.second);
1281  fOneDHistsRecoStop[kPlaneHitX] ->Fill(pc.first);
1282  fTwoDHistsRecoStop[kCellPlaneX]->Fill(pc.first, pc.second);
1283  }
1284  }
1285  else if(view == geo::kY){
1286  fOneDHistsRecoAll[kCellHitY] ->Fill(pc.second);
1287  fOneDHistsRecoAll[kPlaneHitY] ->Fill(pc.first);
1288  fTwoDHistsRecoAll[kCellPlaneY]->Fill(pc.first, pc.second);
1289  if( fSelection.CleanSample(ri) ){
1290  fOneDHistsRecoStop[kCellHitY] ->Fill(pc.second);
1291  fOneDHistsRecoStop[kPlaneHitY] ->Fill(pc.first);
1292  fTwoDHistsRecoStop[kCellPlaneY]->Fill(pc.first, pc.second);
1293  }
1294  }
1295 
1296  const rb::RecoHit rhit = track.RecoHit(chit);
1297  if(!rhit.IsCalibrated()) continue;
1298 
1299  point = TVector3(rhit.X(), rhit.Y(), rhit.Z());
1300  distFromEnd = track.DistanceFromEnd(rhit.Z());
1301 
1302  if(view == geo::kX){
1303  fOneDHistsRecoAll[kResidualsX] ->Fill(residuals[c]);
1304  fOneDHistsRecoAll[kTNSX] ->Fill(chit->TNS());
1305  fTwoDHistsRecoAll[kResidVsCosX]->Fill(residuals[c], ri.cosTheta);
1306  fTwoDHistsRecoAll[kResidVsAzX] ->Fill(residuals[c], ri.azimuth);
1307  }
1308  else if(view == geo::kY){
1309  fOneDHistsRecoAll[kResidualsY] ->Fill(residuals[c]);
1310  fOneDHistsRecoAll[kTNSY] ->Fill(chit->TNS());
1311  fTwoDHistsRecoAll[kResidVsCosY]->Fill(residuals[c], ri.cosTheta);
1312  fTwoDHistsRecoAll[kResidVsAzY] ->Fill(residuals[c], ri.azimuth);
1313  }
1314  fOneDHistsRecoAll[kHitXPosition] ->Fill(rhit.X());
1315  fOneDHistsRecoAll[kHitYPosition] ->Fill(rhit.Y());
1316  fOneDHistsRecoAll[kHitZPosition] ->Fill(rhit.Z());
1317  fTwoDHistsRecoAll[kHitXYPosition]->Fill(rhit.X(), rhit.Y());
1318  fTwoDHistsRecoAll[kHitXZPosition]->Fill(rhit.Z(), rhit.X());
1319  fTwoDHistsRecoAll[kHitYZPosition]->Fill(rhit.Z(), rhit.Y());
1320  if( fSelection.CleanSample(ri) ){
1321  if(distFromEnd >= 0. && distFromEnd <= 20000.){
1322  // find the pathlength through the cell
1323  try{
1324  lengthInCell = fUtilities.PathLengthInCell(zToPlaneBounds, point, pc);
1325  calibLengthInCell = calib::PathLengthInCell(calibZToPlaneBounds, point, pc);
1326  }
1327  catch(cet::exception &e){
1328  MF_LOG_WARNING("CosmicTrackAna")
1329  << "caught exception\n"
1330  << e
1331  << "\n when trying to calculate tracklength in cell. Set it to a std::numeric_limits<float>::min()";
1332  lengthInCell = std::numeric_limits<float>::min();
1333  }
1334 
1335  fCosmicHit.pathLength = lengthInCell;
1336  fCosmicHit.calibPathLength = calibLengthInCell;
1337  fCosmicHit.truePathLength = trueLengthInCell;
1338  fCosmicHit.triCellPathLength = (triCellMap[pc]) ? fUtilities.TriCellPathLength(track, pc.first, cellWidth, view) : -9999.;
1339  fCosmicHit.PE = chit->PE();
1340  fCosmicHit.PECorr = rhit.PECorr();
1341  fCosmicHit.distFromEnd = distFromEnd;
1342  fCosmicHit.view = view;
1343  fCosmicHit.cell = pc.second;
1344  fCosmicHit.plane = pc.first;
1345  fCosmicHit.module = std::floor(fCosmicHit.cell/32.);
1346  if(view == geo::kX){
1347  fCosmicHit.w = rhit.Y();
1348  fTwoDHistsRecoStop[kPEPerCMVsWX]->Fill(rhit.Y(), rhit.PECorr()/lengthInCell);
1349  }
1350  else if(view == geo::kY){
1351  fCosmicHit.w = rhit.X();
1352  fTwoDHistsRecoStop[kPEPerCMVsWY]->Fill(rhit.X(), rhit.PECorr()/lengthInCell);
1353  }
1354  fHitTree->Fill();
1355  } // end if in the correct window from end point of track
1356  if(view == geo::kX){
1357  fOneDHistsRecoStop[kResidualsX] ->Fill(residuals[c]);
1358  fOneDHistsRecoStop[kTNSX] ->Fill(chit->TNS());
1359  fTwoDHistsRecoStop[kResidVsCosX]->Fill(residuals[c], ri.cosTheta);
1360  fTwoDHistsRecoStop[kResidVsAzX] ->Fill(residuals[c], ri.azimuth);
1361  }
1362  else if(view == geo::kY){
1363  fOneDHistsRecoStop[kResidualsY] ->Fill(residuals[c]);
1364  fOneDHistsRecoStop[kTNSY] ->Fill(chit->TNS());
1365  fTwoDHistsRecoStop[kResidVsCosY]->Fill(residuals[c], ri.cosTheta);
1366  fTwoDHistsRecoStop[kResidVsAzY] ->Fill(residuals[c], ri.azimuth);
1367  }
1368  fOneDHistsRecoStop[kHitXPosition] ->Fill(rhit.X());
1369  fOneDHistsRecoStop[kHitYPosition] ->Fill(rhit.Y());
1370  fOneDHistsRecoStop[kHitZPosition] ->Fill(rhit.Z());
1371  fTwoDHistsRecoStop[kHitXYPosition]->Fill(rhit.X(), rhit.Y());
1372  fTwoDHistsRecoStop[kHitXZPosition]->Fill(rhit.Z(), rhit.X());
1373  fTwoDHistsRecoStop[kHitYZPosition]->Fill(rhit.Z(), rhit.Y());
1374  }
1375 
1376  } // end for cells
1377 
1378 // if(cosCosmic < 0.1){
1379 // MF_LOG_VERBATIM("CosmicTrackAna") << "track params " << track.TotalLength() << " " << cosCosmic << " " << azimuth
1380 // << "\ndirection " << dir.X() << " " << dir.Y() << " " << dir.Z()
1381 // << "\nminPlane: " << track.MinPlane() << " maxPlane: " << track.MaxPlane();
1382 // for(size_t p = 0; p < track.NTrajectoryPoints(); ++p){
1383 // TVector3 point = track.TrajectoryPoint(p);
1384 // MF_LOG_VERBATIM("CosmicTrackAna") << point.X() << " " << point.Y() << " " << point.Z();
1385 // }
1386 // }
1387 
1388  return;
1389  }
1390 
1391  //......................................................................
1393  {
1394 
1395  fTrackLabel = pset.get< std::string >("TrackLabel" );
1396  fExposureLabel = pset.get< std::string >("ExposureLabel" );
1397  fMinSig = pset.get< double >("MinSig" );
1398  fFindPurityAndEfficiency = pset.get< bool >("FindPurityAndEfficiency", false);
1399 
1400  fSelection.reconfigure(pset.get<fhicl::ParameterSet>("SelectionCriteria"));
1401  fUtilities.reconfigure(pset.get<fhicl::ParameterSet>("Utilities" ));
1402 
1403  return;
1404  }
1405 
1406  //......................................................................
1408  {
1409  fExposure->Fill(fLiveTime);
1410 
1411  if(fNevents != 0 && fLiveTime != 0){
1412  float cosmicrate = 1.e-3*fNTracks/fLiveTime;
1413 
1414  MF_LOG_VERBATIM("CosmicTrackAna") << "There are: " << fNTracks
1415  << " cosmic tracks selected in this job.\n"
1416  << "Livetime is " << fLiveTime
1417  << " according to CosmicExposure.\n"
1418  << "Note, livetime is for FULL SUBRUN, "
1419  << "so rate is wrong if haven't run full subrun.\n"
1420  << "The cosmic rate is therefore "
1421  << cosmicrate << " Hz \n";
1422 
1423  fCosmicRate->Fill(cosmicrate);
1424  }
1425 
1426  return;
1427  }
1428 
1429 
1430 
1431 } // end namespace trk
1432 
1433 namespace trk{
1435 }
double E(const int i=0) const
Definition: MCParticle.h:232
difference between true and reco length in cm
std::pair< CosmicTrackInfo, CosmicSignalHitFraction > FillRecoInfo(rb::Track const &track, rb::Cluster const &slice)
truth start position in x
bool fFindPurityAndEfficiency
because the calculation is so slow
float TNS() const
Definition: CellHit.h:46
direction cosine in x
T max(const caf::Proxy< T > &a, T b)
float LengthCorrection(float const &length)
CosmicTrackAna(fhicl::ParameterSet const &pset)
void reconfigure(const fhicl::ParameterSet &pset)
double HalfL() const
Definition: CellGeo.cxx:198
unsigned int NumberTrajectoryPoints() const
Definition: MCParticle.h:217
back track the reconstruction to the simulation
double fMinSig
track - using PE values
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
float TriCellPathLength(rb::Track const &track, uint32_t const &plane, double const &cellWidth, geo::View_t const &view)
unsigned int fNevents
total number of events counter
track start position in x
int PdgCode() const
Definition: MCParticle.h:211
distribution of x/y cellhit times in nano seconds
fractional length vs difference in end plane in each view
difference in end point of recoed track & simulated vs cosTheta
double HalfD() const
Definition: CellGeo.cxx:205
difference in reco and true range momentum as a function of reco length
azimuthal angle of the cosmic ray
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
std::string fExposureLabel
Where to find CosmicExposure info.
void reconfigure(fhicl::ParameterSet const &pset)
distribution of x/y cellhit times in nano seconds
std::map< std::string, double > xmax
PE/cm deposited as a function of W in Y view.
std::vector< TH2D * > fTwoDHistsTruthAll
two dimensional histograms for all tracks
std::string fTrackLabel
Where to find Tracks to reconstruct.
int GetPlaneID() const
Plane ID.
Definition: FLSHit.h:37
track start position in zx
signal size of calibrated cell hits
difference in reco and true range momentum as a function of azimuth
zBoundMap MakeZBoundaryMap(std::set< double > const &planeZBounds, std::vector< TVector3 > const &trajectory)
Return a map of the z position of each trajectory point on a track to the bounding positions of where...
Definition: CalibUtil.cxx:354
std::map< double, zBounds > zBoundMap
difference in end point of recoed track & simulated vs dcosZ
void Make1DRecoHistograms(std::vector< TH1D * > &hists, TString const sample)
difference in track start and end position in y
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:748
Direction cosine in x.
CosmicTrackSelection fSelection
object encapsulating selection criteria
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
art::ServiceHandle< cheat::BackTracker > fBT
int GetCellID() const
Cell ID.
Definition: FLSHit.h:39
double fDetHalfHeight
detector half height, cached
Fraction of hits in event used in track.
std::vector< TH2D * > fTwoDHistsTruthStop
two dimensional histograms for stopping tracks
geo::View_t View() const
Definition: CellHit.h:41
distance in cm to closest edge of detector at track vertex
efficiency in selecting hits for the track
void FillTrackMCTruthHistograms(CosmicTrackInfo const &ti, CosmicTrackInfo const &ri)
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
reco hit positions in xz
cells hit per plane in Y vs X
difference in start point of recoed track & simulated vs cosTheta
std::vector< sim::FLSHit > ParticleToFLSHit(const int &trackID) const
All the FLSHits that were created by the track id trackID, sorted from most to least light...
Fraction difference in direction cosine in y.
difference in recon and true length as a function of reco length
residual for each hit vs azimuth
residuals for the tracks
Vertical planes which measure X.
Definition: PlaneGeo.h:28
_oneDRecoHists
Cosmic track analysis module CosmicTrackAna. Makes some validation plots.
TH1D * fCosmicRate
histogram holding single value for cosmic rate in Hz at end of job
void FindTriCells(rb::Track const &track, std::map< std::pair< uint32_t, uint32_t >, bool > &planeCellMap)
difference in start point of recoed track & simulated vs dcosZ
float Z() const
Definition: RecoHit.h:38
Definition: event.h:19
double HalfW() const
Definition: CellGeo.cxx:191
A collection of associated CellHits.
Definition: Cluster.h:47
track end position in zx
double DetLength() const
PE/cm deposited as a function of W in X view.
CosmicTrackInfo fTrackReco
single track reco information
track end position in xy
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
difference in reco momentum from range and true momentum
largest trajectory step relative to the median vs azimuth
difference in track start and end position in z
Fraction difference in direction cosine in x vs z length.
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
A single unit of energy deposition in the liquid scintillator.
Definition: FLSHit.h:19
number of cells in a track
const PlaneGeo * Plane(unsigned int i) const
void Make1DTruthHistograms(std::vector< TH1D * > &hists, TString const sample)
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
bool Stopper(TVector3 const &finalPoint)
TString hists[nhists]
Definition: bdt_com.C:3
EDAnalyzer(fhicl::ParameterSet const &pset)
Definition: EDAnalyzer.h:27
int fNTracks
total number of reconstructed tracks counter
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
CosmicTrackUtilities fUtilities
object encapsulating utilitie methods
std::vector< TH1D * > fOneDHistsTruthStop
one dimensional histograms for stopping tracks
difference in track start and end position in x
residuals for the tracks
double dE
void Make2DRecoHistograms(std::vector< TH2D * > &hists, TString const sample)
std::vector< TH2D * > fTwoDHistsRecoStop
two dimensional histograms for stopping tracks
float abs(float number)
Definition: d0nt_math.hpp:39
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
CosmicTrackInfo FillTrueInfo(rb::Track const &track, rb::Cluster const &slice, std::map< geo::CellUniqueId, float > &trueCellPathLength)
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
number of tracks found
std::pair< std::pair< float, float >, std::pair< TVector3, TVector3 > > TrueLengthEndPoints(simb::MCParticle const &part)
position of reco hits in z
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
average number of cells hit per plane in y
cosine theta of cosmic track with respect to zenith
void Make2DTruthHistograms(std::vector< TH2D * > &hists, TString const sample)
number of tracks found in event vs particles leaving tracks
const CellUniqueId & Id() const
Definition: CellGeo.h:55
Fraction difference in direction cosine in z.
cosine theta of cosmic track with respect to zenith
Double_t ymax
Definition: plot.C:25
double fCellHalfDepthZ
extent of a cell in the z direction
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
double fDetHalfWidth
detector half width, cached
track end position in y
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
azimuthal angle of the cosmic ray
CosmicTrackTree fTrigger
event information
track signal (Q)
True energy of particle zero (GeV)
std::map< geo::CellUniqueId, float > TruePathLengthInCell(simb::MCParticle const &part, std::vector< sim::FLSHit > const &flsHits)
TString part[npart]
Definition: Style.C:32
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
track end position in zy
double fLiveTime
Livetime in sec according to file.
Fraction difference in direction cosine in x.
double dx[NP][NC]
position of reco hits in y
const sim::Particle * Primary(const int) const
const int xbins
Definition: MakePlots.C:82
std::vector< TH1D * > fOneDHistsTruthAll
one dimensional histograms for all tracks
track start position in z
Track parameters (s, X0, KE, beta, ...) along a particular path in the detector.
Definition: Path.h:18
difference in start point of recoed track & simulated vs azimuth
art::ServiceHandle< geo::LiveGeometry > fLiveGeom
correl_yv Fill(-(dy[iP-1][iC-1]), hyv->GetBinContent(iP, iC))
truth start position in y
const double a
void FilldEdxHistogram(simb::MCParticle const &part)
average number of cells hit per plane in x
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
track start position in zy
Fraction of hits in event used in track.
fractional length vs difference in start plane in each view
direction cosine in z
float PE() const
Definition: CellHit.h:42
void reconfigure(fhicl::ParameterSet const &pset)
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
std::vector< float > TrackResiduals(rb::Track const &track)
float Azimuth(float const &dcosX, float const &dcosZ, float const &magnitude)
reco hit positions in yz
Collect Geo headers and supply basic geometry functions.
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
caf::StandardRecord * sr
SubRunNumber_t subRun() const
fractional length vs asymmetry in number of planes in each view
art::ServiceHandle< art::TFileService > fTFS
difference in z endpoint of recoed track & simulated vs dcosZ
dE/dx vs remaining length in a track
float CosTheta(float const &dcosY, float const &magnitude)
double DetHalfHeight() const
RunNumber_t run() const
std::vector< TH2D * > fTwoDHistsRecoAll
two dimensional histograms for all tracks
double median(TH1D *h)
Definition: absCal.cxx:524
#define MF_LOG_INFO(category)
TTree * fHitTree
initializer of Hit TTree for this module
double fDetLength
detector length, cached
track end position in x
bool sortTVector3ByY(TVector3 const &a, TVector3 const &b)
length of track
direction in x versus y
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...
Fraction difference in direction cosine in z vs z length.
bool sortTVector3ByZ(TVector3 const &a, TVector3 const &b)
void FillEventMCTruthHistograms(size_t trkMax)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
purity in selecting hits for the track
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
const int ybins
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
const std::string path
Definition: plot_BEN.C:43
Fraction difference in direction cosine in y vs z length.
EventNumber_t event() const
Definition: EventID.h:116
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
difference in end point of reconstructed track vs simulated
bool sortFLSHits(sim::FLSHit const &a, sim::FLSHit const &b)
void analyze(art::Event const &evt)
Fraction hits used vs the number of hits that could be used.
reco hit positions in xy
truth start position in z
difference in reco momentum from range uncorrected for track length
double DetHalfWidth() const
int fNTracksMC
total number of MC particles
void FillTrackHistograms(CosmicTrackInfo const &ri, CosmicSignalHitFraction const &shf, rb::Track const &track, std::map< geo::CellUniqueId, float > const &trueCellPathLength)
distance in cm to closest edge of detector at track end
#define MF_LOG_VERBATIM(category)
zBoundMap MakeZBoundaryMap(std::vector< TVector3 > const &trajectory)
difference in z start point of recoed track & simulated vs dcosZ
float X() const
Definition: RecoHit.h:36
track start position in xy
momentum of track from range
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
difference in start point of reconstructed track vs simulated
bool IsCalibTrack(rb::Track const &track, rb::Cluster const &slice)
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
std::vector< TH1D * > fOneDHistsRecoAll
one dimensional histograms for all tracks
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:137
truth end position in x
int fSubRun
Subrun number.
CosmicHitInfo fCosmicHit
information about hits
residual for each hit vs cos(theta)
#define MF_LOG_DEBUG(id)
std::vector< TVector3 > const & Trajectory() const
return a constant reference to the track trajectory points
Definition: Track.h:86
float Mag() const
const hit & b
Definition: hits.cxx:21
fractional length vs maximum cells per plane
CosmicTrackInfo fTrackTruth
single track mc truth information
art::ServiceHandle< geo::Geometry > fGeom
bool GoodReco(rb::Track const &track)
double TotalGeV(EEnergyCalcScheme escheme=kRecomputeEnergy) const
Simple sum of the estimated GeV of all the hits.
Definition: Cluster.cxx:378
direction cosine in y
virtual void beginSubRun(const art::SubRun &sr)
float DistToEdge(TVector3 const &point)
difference in reco and true range momentum as a function of reco length
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
truth end position in z
difference in reco and true range momentum as a function of cos theta
float Y() const
Definition: RecoHit.h:37
std::vector< TH1D * > fOneDHistsRecoStop
one dimensional histograms for stopping tracks
float PECorr() const
Definition: RecoHit.cxx:47
bool CleanSample(CosmicTrackInfo const &ri)
difference in end point of recoed track & simulated vs azimuth
track start position in y
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
#define MF_LOG_WARNING(category)
Int_t trackID
Definition: plot.C:84
void FindZBoundaries(std::set< double > &planeZBounds)
Find the boundaries in the z direction of planes in the detector.
Definition: CalibUtil.cxx:332
track end position in z
TTree * fTree
initializer of TTree for this module
Float_t e
Definition: plot.C:35
position of reco hits in x
TH1D * fExposure
histogram holding single value for exposure in s at end of job
truth end position in y
Direction cosine in z.
virtual double DistanceFromEnd(double z) const
Definition: Track.cxx:255
bool GoodSteps(rb::Track const &track, std::vector< float > &steps, float &median)
float PathLengthInCell(zBoundMap const &zBounds, TVector3 const &recoHitLoc, std::pair< uint32_t, uint32_t > const &pc)
difference in z endpoint of recoed track & simulated vs dcosZ
EventID id() const
double Weight(unsigned int globalIdx) const
Weight assigned to the cell.
Definition: Cluster.cxx:209
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:190
CosmicSignalHitFraction fSigHitFrac
fraction of signal and hits used
float RangeMomentum(rb::Track const &track)
Direction cosine in y.
const Binning zbins
std::vector< const sim::Particle * > HitsToParticle(const std::vector< const rb::CellHit * > &hits) const
Returns vector of sim::Particle objects contributing to the given collection of hits.
float PathLengthInCell(zBoundMap const &zBounds, TVector3 const &recoHitLoc, std::pair< uint32_t, uint32_t > const &pc)
Return the path length of a track in the cell in question.
Definition: CalibUtil.cxx:498
difference in z end point of reconstructed track vs simulated
size of steps along trajectory normalized by the median step
enum BeamMode string