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
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  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  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  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  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  LOG_DEBUG("CosmicTrackAna") << "----------------Start find track path length in cell--------------";
694  this->FillTrackHistograms(fTrackReco, fSigHitFrac, track, trueCellPathLength);
695  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  // 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  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  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  // 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  LOG_DEBUG("CosmicTrackAna") << "------------Fill true cell path lengths-------------";
844  trueCellPathLength = fUtilities.TruePathLengthInCell(*part, flsHits);
845  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  LOG_WARNING("CosmicTrackAna") << "MC dE is negative!! " << dE
1022  << " " << dx
1023  << " " << part.E(t) << " " << part.E(t-1);
1024  if(dx <= 0.){
1025  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  // 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  //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  // 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  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;
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 // 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 // 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  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)
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
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
SubRunNumber_t subRun() const
Definition: Event.h:72
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:744
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.
TH1 * ratio(TH1 *h1, TH1 *h2)
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
bool Stopper(TVector3 const &finalPoint)
TString hists[nhists]
Definition: bdt_com.C:3
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
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
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
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
std::vector< TH2D * > fTwoDHistsRecoAll
two dimensional histograms for all tracks
double median(TH1D *h)
Definition: absCal.cxx:524
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...
#define LOG_WARNING(category)
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.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
reco hit positions in xy
truth start position in z
T * make(ARGS...args) const
difference in reco momentum from range uncorrected for track length
double DetHalfWidth() const
int fNTracksMC
total number of MC particles
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
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
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
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
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:133
truth end position in x
int fSubRun
Subrun number.
CosmicHitInfo fCosmicHit
information about hits
residual for each hit vs cos(theta)
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
#define LOG_VERBATIM(category)
std::vector< TH1D * > fOneDHistsRecoStop
one dimensional histograms for stopping tracks
float PECorr() const
Definition: RecoHit.cxx:47
bool CleanSample(CosmicTrackInfo const &ri)
#define LOG_INFO(stream)
Definition: Messenger.h:144
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
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
RunNumber_t run() const
Definition: Event.h:77
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
Definition: Event.h:56
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:196
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