PCHitsList_module.cc
Go to the documentation of this file.
1 // \file PCHitsList_module.cc
2 // \brief Module to fill vector of PCHits, calibration hit info from cosmic files
3 // \version $Id: PCHitsList_module.cc,v 1.0 2013/01/13 18:20:04 rtoner Exp $
4 // \author rtoner
5 // \date $Date: 2013/01/13 22:06:10 $
6 
7 //C++ and C:
8 #include <cmath>
9 #include <iostream>
10 #include <vector>
11 #include <utility>
12 
13 //Framework:
18 #include "cetlib_except/exception.h"
20 
21 //NovaSoft:
24 #include "DAQChannelMap/DAQChannelMap.h"
25 #include "Geometry/Geometry.h"
26 #include "Geometry/LiveGeometry.h"
28 #include "NovaDAQConventions/DAQConventions.h"
29 #include "Utilities/func/MathUtil.h"
31 #include "MCCheater/BackTracker.h"
32 #include "RecoBase/HitMap.h"
33 #include "RecoBase/Track.h"
34 #include "RecoBase/RecoHit.h"
35 #include "Utilities/AssociationUtil.h"
37 #include "Simulation/FLSHit.h"
38 
39 #include "TVector3.h"
40 
41 namespace caldp { class PCHit; }
42 
43 namespace calib {
44 
45  class PCHitsList : public art::EDProducer {
46 
47  public:
48  explicit PCHitsList(fhicl::ParameterSet const & pset);
49  virtual ~PCHitsList();
50 
51  void produce(art::Event & e) override;
52 
53  void reconfigure(fhicl::ParameterSet const & p);
54 
55  private:
56  bool IsGoodTrack(rb::Track const& track,
57  rb::Cluster const& slice);
58 
59  bool GoodSteps(rb::Track const& track);
60 
61  void GetTrueEnergyPathAndLightForCell(const rb::Track* track,
62  size_t const& ihit,
63  double & mev,
64  double & truePath,
65  double & trueW,
66  double & truePE,
67  double & poissonLambda) const;
68 
69  double TrueW( const sim::FLSHit& fls ) const;
70 
71  void ProcessTrack(const rb::Track* track,
73  std::vector<caldp::PCHit>* pcHitsXY,
74  std::vector<caldp::PCHit>* pcHitsZ,
75  std::vector<caldp::PCHit>* pcHitsAvg) const;
76 
77  void ProcessTrackTrajectory(const rb::Track* track,
79  std::set<double> const& zPlaneBounds,
80  std::vector<caldp::PCHit>* pcHitsTraj) const;
81 
82  void ProcessTrackForBelowThresholdHits(const rb::Track* track,
84  std::vector<caldp::PCHit>* pcHitsXYThresh) const;
85 
86  void ProcessTrackForFLSHits(const rb::Track* track2,
88  std::vector<caldp::PCHit>* pcHitsXY) const;
89 
90  caldp::PCHit CellHitToPCHit(const rb::CellHit* chit,
91  const rb::Track* trk,
92  double const& path,
93  double const& mev,
94  double const& truePE,
95  double const& truePath,
96  double const& trueW,
97  double const& poissonLambda,
98  daqchannelmap::DAQChannelMap* cmap) const;
99 
107 
108  bool fIsBeam;
109  double fPathCut; ///< Only keep hits with pathlength less than this
110  double fMinPath; ///< Only keep hits with pathlength greater than this
111  int fExtentZCut; ///< Only keep tracks crossing this distance in z
112  double fDirZCut; ///< Only keep tracks shallower than this
113  double fCompletenessCut; ///< Only keep tracks with more than this fraction of slice hits included in each view
114  double fMaxVtxDistFromEdge; ///< Only keep tracks with an initial point within this dist (cm) from the edge
115  double fMaxEndDistFromEdge; ///< Only keep tracks whose final point is not more than this distance from the edge
116  double fBadStepSizeLimit; ///< Only keep tracks without steps outside this limit, see GoodSteps comments
117  double fMaxCellsPerPlane; ///< Only keep tracks with fewer <cells/plane> than this limit
118  double fMaxDeltaEndPlane; ///< Only keep tracks where the difference in end plane between the views is less
119  ///< than this number
120  double fMaxPlaneAsymmetry; ///< Only keep tracks where the asymmetry in the planes in the views is less than
121  ///< this number
122 
127  };
128 }
129 
130 namespace calib {
131 
132  //---------------------------------------------------------------------//
133 
134  PCHitsList::PCHitsList(fhicl::ParameterSet const & pset)
135  : EDProducer(pset)
136  {
137  reconfigure(pset);
138 
139  produces< std::vector<caldp::PCHit> >(fQualXYName);
140  produces< std::vector<caldp::PCHit> >(fQualZName);
141  produces< std::vector<caldp::PCHit> >(fQualAvgName);
142  produces< std::vector<caldp::PCHit> >(fQualTrajName);
143  produces< std::vector<caldp::PCHit> >(fQualXYThreshName);
144  produces< std::vector<caldp::PCHit> >(fQualXYFLSName);
145 
146  produces< art::Assns<caldp::PCHit, rb::Track> >(fQualXYName);
147  produces< art::Assns<caldp::PCHit, rb::Track> >(fQualZName);
148  produces< art::Assns<caldp::PCHit, rb::Track> >(fQualAvgName);
149  produces< art::Assns<caldp::PCHit, rb::Track> >(fQualTrajName);
150  produces< art::Assns<caldp::PCHit, rb::Track> >(fQualXYThreshName);
151  produces< art::Assns<caldp::PCHit, rb::Track> >(fQualXYFLSName);
152  }
153 
154  //---------------------------------------------------------------------//
155 
157  {
158  }
159 
160  //---------------------------------------------------------------------//
161 
163  {
164  fCosmicLabel = pset.get< std::string >("CosmicLabel"); // Label of tracker module
165  fQualXYName = pset.get< std::string >("QualXYName"); // Instance label, "XY" quality hits
166  fQualZName = pset.get< std::string >("QualZName"); // Instance label, "Z" quality hits
167  fQualAvgName = pset.get< std::string >("QualAvgName"); // Instance label, "Avg" quality hits
168  fQualTrajName = pset.get< std::string >("QualTrajName"); // Instance label, "Traj" quality hits
169  fQualXYThreshName = pset.get< std::string >("QualXYThreshName"); // Instance label, "XY" quality hits
170  // falling below threshold
171  fQualXYFLSName = pset.get< std::string >("QualXYFLSName"); // Instance label, "XY" quality hits from FLS
172 
173  fIsBeam = pset.get< bool >("IsBeam");
174  fPathCut = pset.get< double >("PathCut");
175  fMinPath = pset.get< double >("MinPathlength");
176  fExtentZCut = pset.get< int >("ExtentZCut");
177  fDirZCut = pset.get< double >("DirZCut");
178  fCompletenessCut = pset.get< double >("CompletenessCut");
179  fMaxVtxDistFromEdge = pset.get< double >("MaxVtxDistFromEdge", 10.);
180  fMaxEndDistFromEdge = pset.get< double >("MaxEndDistFromEdge", 10.);
181  fBadStepSizeLimit = pset.get< double >("BadStepSizeLimit", 2. );
182  fMaxCellsPerPlane = pset.get< double >("MaxCellsPerPlane", 6. );
183  fMaxDeltaEndPlane = pset.get< double >("MaxDeltaEndPlane", 3. );
184  fMaxPlaneAsymmetry = pset.get< double >("MaxPlaneAsymmetry", 0.1);
185 
186  return;
187  }
188 
189  //---------------------------------------------------------------------//
190 
192  {
193  //Ptrs to vectors to be filled (one instance for each degree of hit quality)
194  std::unique_ptr< std::vector<caldp::PCHit> > pcHitsXY (new std::vector<caldp::PCHit>);
195  std::unique_ptr< std::vector<caldp::PCHit> > pcHitsZ (new std::vector<caldp::PCHit>);
196  std::unique_ptr< std::vector<caldp::PCHit> > pcHitsAvg (new std::vector<caldp::PCHit>);
197  std::unique_ptr< std::vector<caldp::PCHit> > pcHitsTraj (new std::vector<caldp::PCHit>);
198  std::unique_ptr< std::vector<caldp::PCHit> > pcHitsXYThresh(new std::vector<caldp::PCHit>);
199  std::unique_ptr< std::vector<caldp::PCHit> > pcHitsXYFLS (new std::vector<caldp::PCHit>);
200 
201  std::unique_ptr<art::Assns<caldp::PCHit, rb::Track> > hitAssnsXY (new art::Assns<caldp::PCHit, rb::Track>);
202  std::unique_ptr<art::Assns<caldp::PCHit, rb::Track> > hitAssnsZ (new art::Assns<caldp::PCHit, rb::Track>);
203  std::unique_ptr<art::Assns<caldp::PCHit, rb::Track> > hitAssnsAvg (new art::Assns<caldp::PCHit, rb::Track>);
204  std::unique_ptr<art::Assns<caldp::PCHit, rb::Track> > hitAssnsTraj (new art::Assns<caldp::PCHit, rb::Track>);
205  std::unique_ptr<art::Assns<caldp::PCHit, rb::Track> > hitAssnsXYThresh(new art::Assns<caldp::PCHit, rb::Track>);
206  std::unique_ptr<art::Assns<caldp::PCHit, rb::Track> > hitAssnsXYFLS (new art::Assns<caldp::PCHit, rb::Track>);
207 
208  //Get cosmic tracks to use for calibration:
210  e.getByLabel(fCosmicLabel, trackcol);
211 
212  std::vector< art::Ptr<rb::Track> > tracks;
213  art::fill_ptr_vector(tracks, trackcol);
214 
215  art::FindManyP<rb::Cluster> fmpslc(trackcol, e, fCosmicLabel);
216 
218 
219  // prepare to find the path lengths in a cell based on the track trajectory
220  std::set<double> planeZBounds;
221  calib::FindZBoundaries(planeZBounds);
222 
223  for(size_t itrk = 0; itrk < tracks.size(); ++itrk){
224 
225  // Find the slice associated with this track
226  // There should be exactly one
227  const std::vector<art::Ptr<rb::Cluster> > slices = fmpslc.at(itrk);
228  if(slices.size() != 1)
229  throw cet::exception("PCHitsList") << "incorrect number of slices associated to track";
230 
231  const rb::Cluster& slice = *slices[0];
232 
233  // Remember how many hits of each type we had so far, so we know what
234  // index to start associating from when we're done with this track.
235  const size_t xyBeginIdx = pcHitsXY->size();
236  const size_t zBeginIdx = pcHitsZ->size();
237  const size_t avgBeginIdx = pcHitsAvg->size();
238  const size_t trajBeginIdx = pcHitsTraj->size();
239  const size_t xyThreshBeginIdx = pcHitsXYThresh->size();
240  const size_t xyFLSBeginIdx = pcHitsXYFLS->size();
241 
242  const art::Ptr<rb::Track> track = tracks[itrk];
243 
244  MF_LOG_DEBUG("PCHitsList") << "TRACK " << itrk << ": " << track->NCell() << " cells to process.";
245 
246  //-------Track cuts:----------------//
247  if( !this->IsGoodTrack(*track, slice) ) continue;
248 
249  // track passes the cuts, so process it
250  ProcessTrack(track.get(),
251  cmap,
252  pcHitsXY.get(),
253  pcHitsZ.get(),
254  pcHitsAvg.get());
255 
256  ProcessTrackTrajectory(track.get(),
257  cmap,
258  planeZBounds,
259  pcHitsTraj.get());
260 
262  cmap,
263  pcHitsXYThresh.get());
264 
265  ProcessTrackForFLSHits(track.get(),
266  cmap,
267  pcHitsXYFLS.get());
268 
269 
270  // Associate all the PCHits we created for this track. We do this in
271  // bulk because making lots of individual calls as we go is slow.
272  util::CreateAssn(e, *pcHitsXY, track, *hitAssnsXY, xyBeginIdx, UINT_MAX, fQualXYName);
273  util::CreateAssn(e, *pcHitsZ, track, *hitAssnsZ, zBeginIdx, UINT_MAX, fQualZName);
274  util::CreateAssn(e, *pcHitsAvg, track, *hitAssnsAvg, avgBeginIdx, UINT_MAX, fQualAvgName);
275  util::CreateAssn(e, *pcHitsTraj, track, *hitAssnsTraj, trajBeginIdx, UINT_MAX, fQualTrajName);
276  util::CreateAssn(e, *pcHitsXYThresh, track, *hitAssnsXYThresh, xyThreshBeginIdx, UINT_MAX, fQualXYThreshName);
277  util::CreateAssn(e, *pcHitsXYFLS, track, *hitAssnsXYFLS, xyFLSBeginIdx, UINT_MAX, fQualXYFLSName);
278 
279  }//End track loop
280 
281  MF_LOG_DEBUG("PCHitsList")
282  << "SUMMARY: Size of XY vector = " << pcHitsXY->size()
283  << ", Size of Z vector = " << pcHitsZ->size()
284  << ", Size of Avg vector = " << pcHitsAvg->size();
285 
286  e.put(std::move(pcHitsXY), fQualXYName);
287  e.put(std::move(pcHitsZ), fQualZName);
288  e.put(std::move(pcHitsAvg), fQualAvgName);
289  e.put(std::move(pcHitsTraj), fQualTrajName);
290  e.put(std::move(pcHitsXYThresh), fQualXYThreshName);
291  e.put(std::move(pcHitsXYFLS), fQualXYFLSName);
292  e.put(std::move(hitAssnsXY), fQualXYName);
293  e.put(std::move(hitAssnsZ), fQualZName);
294  e.put(std::move(hitAssnsAvg), fQualAvgName);
295  e.put(std::move(hitAssnsTraj), fQualTrajName);
296  e.put(std::move(hitAssnsXYThresh), fQualXYThreshName);
297  e.put(std::move(hitAssnsXYFLS), fQualXYFLSName);
298 
299  MF_LOG_DEBUG("PCHitsList") << "Event finished successfully.";
300 
301  }//end produce
302 
303 
304  //---------------------------------------------------------------------//
306  rb::Cluster const& slice)
307  {
308  MF_LOG_DEBUG("PCHitsList") << "NCellsX=" << track.NXCell() << " NCellsY=" << track.NYCell();
309 
310  //Track must have at least two X cells and 2 Y cells:
311  if(track.NXCell() < 2 || track.NYCell() < 2) {
312  MF_LOG_DEBUG("PCHitsList") << "!---->Failed NCell Cut!";
313  return false;
314  }
315 
316  //Check that direction vector is well-defined:
317  const TVector3 dr = track.Dir().Unit();
318  const double dx = std::abs(dr.X());
319  const double dy = std::abs(dr.Y());
320  const double dz = std::abs(dr.Z());
321 
322  //No undefined directions:
323  if(std::isnan(dx) ||
324  std::isnan(dy) ||
325  std::isnan(dz)) return false;
326 
327  //At least one direction must be non-zero:
328  if(dx == 0 &&
329  dy == 0 &&
330  dz == 0) {
331  MF_LOG_DEBUG("PCHitsList") << "!---->Failed Direction Cut!";
332  return false;
333  }
334 
335  //Must be a unit vector (otherwise, abort event, something went very wrong):
336  if ( std::abs( util::pythag(dx, dy, dz)-1 ) > 1e-4){
337  throw cet::exception("SkipEvent") << "ERROR: Direction was not a unit vector. Skipping Event.\n";
338  }
339 
340  // Additional quality cuts. These remove tracks that often have
341  // badly-reconstructed W positions.
342 
343  // Cut very steep tracks two ways. These can be misreconstructed, or
344  // they can simply make it difficult for the W estimation to do a good
345  // job. OK to use the rb::Track::Dir() method here because the initial
346  // direction of the track is what we want to cut on
347  if(std::abs(track.Start().Z() - track.Stop().Z() ) <= fExtentZCut ||
348  std::abs(track.Dir().Z()) < fDirZCut) return false;
349 
350  // Cut out actual reconstruction failures, where the track doesn't use
351  // a significant amount of the slice hits. If either view is not
352  // complete enough discard the whole thing
353  if(track.NXCell() < fCompletenessCut*slice.NXCell() ||
354  track.NYCell() < fCompletenessCut*slice.NYCell()) return false;
355 
356  // Check that the average number of cells per plane is similar in
357  // each view. If the value is very high in either view, or larger in
358  // one view than the other then you might have some form of FEB flash
359  // or cosmic ray brem that will mess up the calibration
360  std::set<unsigned short> xPlaneList;
361  std::set<unsigned short> yPlaneList;
362  for(size_t c = 0; c < track.NCell(); ++c){
363  if (track.Cell(c)->View() == geo::kX) xPlaneList.insert(track.Cell(c)->Plane());
364  else if(track.Cell(c)->View() == geo::kY) yPlaneList.insert(track.Cell(c)->Plane());
365  }// end loop over cells in track
366 
367  double xCells = track.NXCell();
368  double yCells = track.NYCell();
369  double xPlanes = 1. * xPlaneList.size();
370  double yPlanes = 1. * yPlaneList.size();
371  double cellsPerPlaneX = (xPlanes > 0) ? xCells/xPlanes : 999.;
372  double cellsPerPlaneY = (yPlanes > 0) ? yCells/yPlanes : 999.;
373  if( std::max(cellsPerPlaneX, cellsPerPlaneY) > fMaxCellsPerPlane ) return false;
374 
375  // remove tracks where there is an asymmetry between the number
376  // of planes hit in each view. One class of asymmetry is where
377  // there are extra planes off either end of the track in one view.
378  // Since there is no information in the other view for those planes
379  // you can't be sure you got the length/containment correct
380  double xStart = track.MaxPlane(geo::kX);
381  double yStart = track.MaxPlane(geo::kY);
382  double xEnd = track.MinPlane(geo::kX);
383  double yEnd = track.MinPlane(geo::kY);
384  if(std::abs(xStart - yStart) > fMaxDeltaEndPlane ||
385  std::abs(xEnd - yEnd ) > fMaxDeltaEndPlane ||
386  std::abs(xPlanes - yPlanes) > fMaxPlaneAsymmetry*(xPlanes + yPlanes) ) return false;
387 
388 
389  // Check that the track doesn't take one or two big steps
390  // compared to the rest, that is usually a reconstruction failure.
391  if(!this->GoodSteps(track) ) return false;
392 
393  int junk = -999;
394 
395  // Remove tracks whose vertex position is well away from the edge of the
396  // detector. Muon catcher is not counted as part of the detector proper in
397  // ND so have to check for that as well.
402  return false;
403  }
404 
405  // Remove tracks whose end position is more than fMaxEndDistFromEdge
406  // outside of the detector The minus sign is because LiveGeometry returns a
407  // negative value if the point is outside of the active detector. A track
408  // can stop anywhere inside the detector and still be OK.
409  if( !fIsBeam && (fLiveGeom->DistToEdgeXY(track.Stop()) < -fMaxEndDistFromEdge ||
412  fLiveGeom->DistanceToEdgeInMC(track.Stop(), junk) < -fMaxEndDistFromEdge) ) {
413  return false;
414  }
415 
416  return true;
417  }
418 
419 
420  //......................................................................
421  // tracks that have huge changes in the step size from one trajectory
422  // point to another are probably poorly reconstructed, and should be removed
423  // from the analysis. Usually the biggest steps are at either the start
424  // or end of the track in cases of really bad reconstruction
426  {
427 
428  // first find the median step size
429  size_t size = track.NTrajectoryPoints() - 1;
430 
431  // track with only two distances between trajectory points,
432  // probably not very useful for the analysis so return false
433  if(size < 3) return false;
434 
435  std::vector<float> stepsize(size, 0.);
436  float largestStep = -1.;
437  for(size_t t = 1; t < size+1; ++t){
438  TVector3 const& tp1 = track.TrajectoryPoint(t-1);
439  TVector3 const& tp2 = track.TrajectoryPoint(t);
440  stepsize[t-1] = (tp1 - tp2).Mag();
441  largestStep = std::max(stepsize[t-1], largestStep);
442  }
443 
444  // sort the vector of step sizes and find the median
445  std::sort(stepsize.begin(), stepsize.end());
446 
447  float median = 0.;
448  if(size%2 == 0)
449  median = 0.5*(stepsize[size/2] + stepsize[(size/2) - 1]);
450  else
451  median = stepsize[size/2];
452 
453  // the median represents the value of the step size
454  // distribution where half of the distribution is
455  // below that value and half is above it. If the
456  // largest step is more than twice the median, then
457  // it is very far outside the distribution and this
458  // track was not well reconstructed
459  if(largestStep > median*fBadStepSizeLimit) return false;
460 
461  return true;
462  }
463 
464  //---------------------------------------------------------------------//
465  // all input references are assumed to have been set to zero by the calling method
467  size_t const& ihit,
468  double & mev,
469  double & truePath,
470  double & trueW,
471  double & truePE,
472  double & poissonLambda) const
473  {
474 
475  if (fBT->HaveTruthInfo()){
476  const std::vector<cheat::TrackIDE> ides = fBT->HitToTrackIDE(track->Cell(ihit));
477 
478  for(size_t n = 0; n < ides.size(); ++n){
479  mev += ides[n].energy*1000;
480 
481  // the ides are sorted to have the most contributing track id be first in the list
482  // so grab that track id and get the particle it corresponds to
483  if(n < 1){
484  auto particle = fBT->ParticleNavigator()[ides[n].trackID];
485 
486  std::vector<sim::FLSHit> flsHits = fBT->ParticleToFLSHit(particle->TrackId(), particle->PdgCode());
487  for(auto fls : flsHits){
488  if(fls.GetPlaneID() == track->Cell(ihit)->Plane() &&
489  fls.GetCellID() == track->Cell(ihit)->Cell() ){
490  truePath = fls.GetTotalPathLength();
491  trueW = TrueW(fls);
492  break;
493  }
494  }
495  } // end if on the first ide
496  } // end loop over ides
497 
498  const std::vector<sim::PhotonSignal> pes = fBT->HitToPhotonSignal(track->Cell(ihit));
499  for(size_t n = 0; n < pes.size(); ++n){
500  truePE += pes[n].NPhoton();
501  poissonLambda += pes[n].PoissonLambda();
502  }
503 
504  }
505 
506  return;
507  }
508 
509 
510 
511  //---------------------------------------------------------------------//
512  double PCHitsList::TrueW(const sim::FLSHit& fls) const
513  {
514  double local_enter[3] = {fls.GetEntryX(), fls.GetEntryY(), fls.GetEntryZ()};
515  double local_exit[3] = {fls.GetExitX(), fls.GetExitY(), fls.GetExitZ()};
516  double enter[3], exit[3];
517  int ip = fls.GetPlaneID();
518  int ic = fls.GetCellID();
519 
520  (fGeom->Plane(ip)->Cell(ic))->LocalToWorld(local_enter,enter);
521  (fGeom->Plane(ip)->Cell(ic))->LocalToWorld(local_exit ,exit);
522 
523  if(fGeom->Plane(fls.GetPlaneID())->View() == geo::kX)
524  return ( enter[1] + exit[1] )/2;
525  else if(fGeom->Plane(fls.GetPlaneID())->View() == geo::kY){
526  // TB horizontal modules (Y-view) have readout on opposite end as ND/FD
528  return -( enter[0] + exit[0] )/2;
529  return ( enter[0] + exit[0] )/2;
530  }
531  else
532  mf::LogError("TrueW") << "bad view for fls hit";
533  return 9999.;
534  }
535 
536 
537  //---------------------------------------------------------------------//
540  std::vector<caldp::PCHit>* pcHitsXY,
541  std::vector<caldp::PCHit>* pcHitsZ,
542  std::vector<caldp::PCHit>* pcHitsAvg) const
543  {
544  std::vector<rb::CellHit> chits;
545  for(unsigned int ihit = 0; ihit < track->NCell(); ++ihit){
546  chits.push_back(*track->Cell(ihit).get());
547  }
548 
549  //Calculate path quality:
550  std::vector<EPathQuality> quals;
551  const std::vector<double> paths = BestPathEstimates(track->OfflineChans(), track->PlaneDirMap(), quals);
552 
553  //Check there is a path and quality measure for each cell; if not, something wrong, abort event:
554  if(chits.size() != quals.size() ||
555  chits.size() != paths.size()){
556  throw cet::exception("SkipEvent") << "ERROR: Number of pathlengths or qualities != number of cells. Skipping Event.\n";
557  }
558 
559  //Loop through cells in track:
560  for(size_t ihit = 0; ihit < chits.size(); ++ihit){
561 
562  //Get hit
563  const rb::CellHit& chit = chits[ihit];
564 
565  //Get path:
566  const double path = paths[ihit]; //RBT: Way of ensuring cell and path line up??
567 
568  //Check path not too long:
569  if(path > fPathCut) continue;
570 
571  //Get the correct MeV and true PE in this hit
572  double mev = 0.;
573  double truePE = 0.;
574  double poissonLambda = 0.;
575  double truePath = 0.;
576  double trueW = -9999.;
577  this->GetTrueEnergyPathAndLightForCell(track, ihit, mev, truePath, trueW, truePE, poissonLambda);
578 
579  caldp::PCHit pchit = CellHitToPCHit(&chit, track,
580  path, mev, truePE, truePath, trueW, poissonLambda,
581  cmap);
582 
583  switch(quals[ihit]){
584 
585  case kXYAdjacents:
586  MF_LOG_DEBUG("PCHitsList")
587  << "Quality XY! - hit=" << ihit << " FILLS TO -> plane=" << chit.Plane()
588  << " cell=" << chit.Cell() << " pe=" << chit.PE() << " path=" << path
589  << " W=" << track->W(&chit) << " mev=" << mev;
590  pcHitsXY->push_back(pchit);
591  break;
592 
593  case kZAdjacents:
594  MF_LOG_DEBUG("PCHitsList")
595  << "Quality Z! - hit=" << ihit << " FILLS TO -> plane=" << chit.Plane()
596  << " cell=" << chit.Cell() << " pe=" << chit.PE() << " path=" << path
597  << " W=" << track->W(&chit) << " mev=" << mev;
598  pcHitsZ->push_back(pchit);
599  break;
600 
601  case kAverage:
602  MF_LOG_DEBUG("PCHitsList")
603  << "Quality Avg! - hit=" << ihit << " FILLS TO -> plane=" << chit.Plane()
604  << " cell=" << chit.Cell() << " pe=" << chit.PE() << " path=" << path
605  << " W=" << track->W(&chit) << " mev=" << mev;
606  pcHitsAvg->push_back(pchit);
607  break;
608 
609  case kBadTrack:
610  MF_LOG_DEBUG("PCHitsList") << "!!!!!!BAD hit/track!!!!!! - hit=" << ihit;
611  // Something wrong with this hit/track, so do nothing
612  break;
613 
614  default:
615  throw cet::exception("SkipEvent") << "ERROR: Unknown Path Quality type. Skipping Event.\n";
616  //Unknown path type; something went very wrong; abort this event.
617  }
618 
619  }//End cell loop
620 
621  }//end ProcessTrack
622 
623  //---------------------------------------------------------------------//
626  std::set<double> const& zPlaneBounds,
627  std::vector<caldp::PCHit>* pcHitsTraj) const
628  {
629  auto boundaryMap = calib::MakeZBoundaryMap(zPlaneBounds, track->Trajectory());
630 
631  // don't bother with this track if no z boundaries were found, that implies
632  // that all of the trajectory points are outside of the live geometry
633  if(boundaryMap.empty()) return;
634 
635  for(unsigned int ihit = 0; ihit < track->NCell(); ++ihit){
636  const rb::CellHit& chit = *track->Cell(ihit);
637 
638  std::pair<uint32_t, uint32_t> pc = {chit.Plane(), chit.Cell()};
639 
640  double xyz[3];
641  fGeom->Plane(chit.Plane())->Cell(chit.Cell())->GetCenter(xyz);
642  xyz[1-chit.View()] = track->W(&chit);
643  // TODO - wouldn't this be a better estimate of the hit position?
644  // track->InterpolateXY(xyz[2], xyz[0], xyz[1]);
645  TVector3 point(xyz[0], xyz[1], xyz[2]);
646 
647  MF_LOG_DEBUG("PCHitsList")
648  << ihit << "/" << track->NCell()
649  << " recohit: " << point.X()
650  << " " << point.Y() << " " << point.Z();
651 
652  // don't bother with the point if it is reconstructed outside of the detector
653  // silly but apparently true, IsPointLive cannot handle the muon catcher
654  // IsMuonCatcher returns false if in the FD, so in the FD this check
655  // should be determined by IsPointLive
656  if( !fLiveGeom->IsPointLive(point) && !fLiveGeom->InMuonCatcher(point) ){
657  MF_LOG_DEBUG("PCHitsList")
658  << ihit
659  << " is not valid "
660  << "(x,y,z) = "
661  << point.X() << "," << point.Y() << "," << point.Z();
662  continue;
663  }
664 
665  try{
666  float path = calib::PathLengthInCell(boundaryMap, point, pc);
667 
668  //Check path not too long:
669  if(path > fPathCut || path < fMinPath) continue;
670 
671  //Get the correct MeV and true PE in this hit
672  double mev = 0.;
673  double truePE = 0.;
674  double trueW = -9999.;
675  double truePath = 0.;
676  double poissonLambda = 0.;
677  this->GetTrueEnergyPathAndLightForCell(track, ihit, mev, truePath, trueW, truePE, poissonLambda);
678 
679  pcHitsTraj->push_back(CellHitToPCHit(&chit, track,
680  path, mev, truePE, truePath, trueW, poissonLambda,
681  cmap));
682  }
683  catch(cet::exception &e){
684  MF_LOG_WARNING("PCHitsList")
685  << "caught exception \n" << e
686  << "\n cell " << pc.first << " " << pc.second << " not added to PCHitsList";
687  }
688 
689  }//End cell loop
690 
691  }//end ProcessTrackTrajectory
692 
693  //---------------------------------------------------------------------//
696  std::vector<caldp::PCHit>* pcHitsXYThresh) const
697  {
698  std::map<geo::OfflineChan, double> trueE;
699  std::map<geo::OfflineChan, double> truePE;
700  std::map<geo::OfflineChan, double> truePath;
701  std::map<geo::OfflineChan, double> trueW;
702  std::map<geo::OfflineChan, double> poissonLambda;
703  if(fBT->HaveTruthInfo()){
704  const std::vector<const sim::Particle*> parts = fBT->HitsToParticle(track->AllCells());
705  if(!parts.empty()){
706  const int id = parts[0]->TrackId(); // Biggest contributor
707  std::set<geo::OfflineChan> flsChans;
708  const std::vector<sim::FLSHit> flss = fBT->ParticleToFLSHit(id);
709  for(const sim::FLSHit& fls: flss){
710  const geo::OfflineChan chan(fls.GetPlaneID(), fls.GetCellID());
711  trueE[chan] += 1000*fls.GetEdep();
712  truePath[chan] = fls.GetTotalPathLength();
713  trueW[chan] = TrueW(fls);
714  flsChans.insert(chan);
715  }
716 
717  for(geo::OfflineChan chan: flsChans){
718  // We won't have an actual CellHit for this cell, so have to use
719  // CellTo not HitTo. Could try to do something clever with matching
720  // the times up. But this must be a small effect.
721  const std::vector<sim::PhotonSignal> pes = fBT->CellToPhotonSignal(chan.Plane(), chan.Cell());
722  for(sim::PhotonSignal pe: pes){
723  truePE[chan] += pe.NPhoton();
724  poissonLambda[chan] += pe.PoissonLambda();
725  }
726  }
727  }
728  }
729 
730  // Find cells that would have made it into the sample if they'd been
731  // hit. Downstream analyzers can interpret these as a threshold effect.
732  std::vector<EPathQuality> threshQuals;
733  std::vector<double> threshPaths;
734  std::vector<geo::OfflineChan> threshChans;
735 
736  FindBelowThresholdCalibCandidates(track->OfflineChans(), track->PlaneDirMap(), threshQuals, threshPaths, threshChans);
737 
738  for(unsigned int threshIdx = 0; threshIdx < threshQuals.size(); ++threshIdx){
739  const double path = threshPaths[threshIdx];
740  // Check path not too long. Apply same cut as normal candidates so as
741  // not to bias ourselves.
742  if(path > 10) continue;
743 
744  const int plane = threshChans[threshIdx].Plane();
745  const int cell = threshChans[threshIdx].Cell();
746  const geo::OfflineChan chan(plane, cell);
747  const geo::View_t view = fGeom->Plane(plane)->View();
749  dummy.SetPlane(plane);
750  dummy.SetCell(cell);
751  dummy.SetView(view);
752 
753  caldp::PCHit pchit = CellHitToPCHit(&dummy, track,
754  path, trueE[chan],
755  truePE[chan], truePath[chan], trueW[chan], poissonLambda[chan],
756  cmap);
757 
758  if(threshQuals[threshIdx] == kXYAdjacents)
759  pcHitsXYThresh->push_back(pchit);
760  } // end for threshIdx
761 
762  } // end ProcessTrackForBelowThresholdHits
763 
764  //---------------------------------------------------------------------//
767  std::vector<caldp::PCHit>* pcHitsXY) const
768  {
769  if(!fBT->HaveTruthInfo()) return;
770 
771  rb::HitMap hmap(track->AllCells());
772 
773  // BackTracker won't be able to figure out the true energy, since
774  // there's no CellHit to work back from. Accumulate it here.
775  std::map<geo::OfflineChan, double> extrasTruthMap;
776  std::map<geo::OfflineChan, double> extrasPEMap;
777  std::map<geo::OfflineChan, double> extrasTruePath;
778  std::map<geo::OfflineChan, double> extrasTrueW;
779  std::map<geo::OfflineChan, double> extrasPoissonLambda;
780  std::vector<geo::OfflineChan> extras;
781 
782  const std::vector<const sim::Particle*> parts = fBT->HitsToParticle(track->AllCells());
783  if(!parts.empty()){
784  const int id = parts[0]->TrackId(); // Biggest contributor
785  const std::vector<sim::FLSHit> flss = fBT->ParticleToFLSHit(id);
786  const unsigned int nFLS = flss.size();
787  std::set<geo::OfflineChan> extraset;
788  std::set<geo::OfflineChan> flsCells;
789  for(unsigned int flsIdx = 0; flsIdx < nFLS; ++flsIdx){
790  const sim::FLSHit& fls = flss[flsIdx];
791  const geo::OfflineChan chan(fls.GetPlaneID(), fls.GetCellID());
792  extraset.insert(chan);
793  extrasTruthMap[chan] += 1000*fls.GetEdep();
794  extrasTruePath[chan] = fls.GetTotalPathLength();
795  if(fGeom->Plane(fls.GetPlaneID())->View() == geo::kX)
796  extrasTrueW[chan] = fls.GetYAverage();
797  else if(fGeom->Plane(fls.GetPlaneID())->View() == geo::kY)
798  extrasTrueW[chan] = fls.GetXAverage();
799  flsCells.insert(chan);
800  }
801 
802  for(geo::OfflineChan chan: flsCells){
803  // We might not have an actual CellHit for this cell, so have to use
804  // CellTo not HitTo. Could try to do something clever with matching
805  // the times up. But this must be a small effect.
806  const std::vector<sim::PhotonSignal> pes = fBT->CellToPhotonSignal(chan.Plane(), chan.Cell());
807  for(sim::PhotonSignal pe: pes){
808  extrasPEMap[chan] += pe.NPhoton();
809  extrasPoissonLambda[chan] += pe.PoissonLambda();
810  }
811  }
812 
813  // extraset is a set to get a collection of distinct cells, but
814  // later on we need it to be a vector.
815  extras.insert(extras.end(), extraset.begin(), extraset.end());
816  }
817 
818  const double meant = track->MeanTNS();
819 
820  rb::Track trackNoHits(*track);
821  trackNoHits.Clear();
822 
823  //Calculate path quality
824  std::vector<EPathQuality> quals;
825  const std::vector<double> paths = BestPathEstimates(trackNoHits.OfflineChans(), trackNoHits.PlaneDirMap(), quals, extras);
826 
827  //Loop through cells in track:
828  for(unsigned int ihit = 0; ihit < extras.size(); ++ihit){
829  const geo::OfflineChan chan = extras[ihit];
830 
831  //Get path:
832  const double path = paths[ihit];
833  //Check path not too long:
834  if(path > 10) continue;
835 
836  if(quals[ihit] == kXYAdjacents){
837  //Get the correct MeV in this hit
838  const double mev = extrasTruthMap[chan];
839  const double truePE = extrasPEMap[chan];
840  const double truePath = extrasTruePath[chan];
841  const double trueW = extrasTrueW[chan];
842  const double poissonLambda = extrasPoissonLambda[chan];
843 
844  rb::CellHit chit;
845  chit.SetPlane(chan.Plane());
846  chit.SetCell(chan.Cell());
847  chit.SetView(fGeom->Plane(chan.Plane())->View());
848  if(hmap.CellExists(chan.Plane(), chan.Cell())){
849  chit.SetPE(hmap.Cell(chan.Plane(), chan.Cell())->PE());
850  }
851  else{
852  chit.SetPE(0);
853  }
854  chit.SetTNS(meant, true);
855 
856  caldp::PCHit pchit = CellHitToPCHit(&chit, track,
857  path, mev, truePE, truePath, trueW, poissonLambda,
858  cmap);
859 
860  pcHitsXY->push_back(pchit);
861  }
862  }//End cell loop
863 
864  }//end ProcessTrackForFLSHit
865 
866  //---------------------------------------------------------------------//
868  const rb::Track* trk,
869  double const& path,
870  double const& mev,
871  double const& truePE,
872  double const& truePath,
873  double const& trueW,
874  double const& poissonLambda,
876  {
877  //get more information on hit location
878  uint32_t lchan = cmap->encodeLChan(fGeom->DetId(), chit->Plane(), chit->Cell());
879  uint32_t dchan = cmap->encodeDChan(lchan);
880  uint32_t pixel = cmap->getPixel(dchan);
881  uint32_t apd = cmap->getFEB(dchan);
882  uint32_t dcm = cmap->getDCM(dchan);
883  uint32_t diblock = cmap->getDiBlock(dchan);
884 
885  double w = trk->W(chit);
886 
887  //get distance to readout electronics
888  double xyzd[4];
889  GetXYZD(chit->OfflineChan(), w, xyzd);
890  //make sure the reconstructed hit coordinate is inside the physical
891  //detector if not then change the distance to the electronics to -1.0
892  if(std::abs(xyzd[0]) > fGeom->DetHalfWidth() ||
893  std::abs(xyzd[1]) > fGeom->DetHalfHeight() ||
894  xyzd[2] > fGeom->DetLength() ||
895  xyzd[2] < 0.0) xyzd[3] = -1.0;
896 
897  double flightlen = trk->DistanceFromStart(xyzd[2]);
898 
899  //Define and fill a PCHit for this cell hit
900  caldp::PCHit pchit;
901  pchit.SetPlane ( chit->Plane() ); //Set plane number
902  pchit.SetCell ( chit->Cell() ); //Set cell number
903  pchit.SetPixel ( pixel ); //Set pixel number
904  pchit.SetAPD ( apd ); //Set apd number
905  pchit.SetDCM ( dcm ); //Set dcm number
906  pchit.SetDiblock ( diblock ); //set diblock number
907  pchit.SetView ( chit->View() ); //set hit view
908  pchit.SetPE ( chit->PE() ); //Set uncalibrated PE
909  pchit.SetPath ( path ); //Set path length
910  pchit.SetW ( w ); //Set W (unmeasured coordinate) value
911  pchit.SetTrueMeV ( mev ); //Set true energy of hit
912  pchit.SetTruePE ( truePE ); //Set true PE of hit
913  pchit.SetTruePath ( truePath ); //Set true pathlength through cell
914  pchit.SetTrueW ( trueW );
915  pchit.SetPoissonLambda(poissonLambda );
916  pchit.SetTNS ( chit->TNS() ); //Set time of hit
917  pchit.SetFlightLen ( flightlen ); //Set path from start of track to hit
918  pchit.SetReadDist ( xyzd[3] ); //Set distance from hit to readout electronics
919  pchit.SetGoodTime ( chit->GoodTiming()); //Set if cell had a proper time fit
920 
921  return pchit;
922  } // end CellHitToPCHit
923 
924 } //end namespace
925 
926 //---------------------------------------------------------------------//
927 
float TNS() const
Definition: CellHit.h:46
T max(const caf::Proxy< T > &a, T b)
float GetEntryX() const
Entry point of the particle (position, time and energy)
Definition: FLSHit.h:48
size_t NTrajectoryPoints() const
Definition: Track.h:83
caldp::PCHit CellHitToPCHit(const rb::CellHit *chit, const rb::Track *trk, double const &path, double const &mev, double const &truePE, double const &truePath, double const &trueW, double const &poissonLambda, daqchannelmap::DAQChannelMap *cmap) const
novadaq::cnv::DetId DetId() const
What detector are we in?
std::vector< geo::OfflineChan > OfflineChans() const
Positions of all the CellHits.
Definition: Cluster.cxx:190
back track the reconstruction to the simulation
static bool CreateAssn(art::EDProducer const &prod, art::Event &evt, std::vector< T > &a, art::Ptr< U > b, art::Assns< T, U > &assn, size_t indx=UINT_MAX, std::string const &instance=std::string())
Create a 1 to 1 association between a new product and one already in the event.
std::string fQualXYThreshName
void SetTNS(float aTNS)
Set hit time (ns)
Definition: PCHit.h:95
diblock
print "ROW IS " print row
Definition: geo2elec.py:31
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void SetPE(float aPE)
Set PE value.
Definition: PCHit.h:79
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
void SetTNS(float tns, bool good)
Definition: CellHit.h:56
int fExtentZCut
Only keep tracks crossing this distance in z.
int GetPlaneID() const
Plane ID.
Definition: FLSHit.h:37
void SetAPD(int aAPD)
Set apd value.
Definition: PCHit.h:73
void SetPixel(int aPixel)
Set pixel value.
Definition: PCHit.h:75
void GetTrueEnergyPathAndLightForCell(const rb::Track *track, size_t const &ihit, double &mev, double &truePath, double &trueW, double &truePE, double &poissonLambda) const
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
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Track.cxx:281
std::vector< TrackIDE > HitToTrackIDE(const rb::CellHit &hit, bool useBirksE=false) const
Convenience function. HitsToTrackIDE but for a single hit.
void SetFlightLen(float aLen)
Set path length from start of track.
Definition: PCHit.h:99
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:748
std::string fQualXYFLSName
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
double DistToEdgeXY(TVector3 vertex)
void SetPoissonLambda(float aLambda)
Set number of simulated photons at readout before fluctuations.
Definition: PCHit.h:103
int GetCellID() const
Cell ID.
Definition: FLSHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
float GetTotalPathLength() const
Get path length of all steps in FLSHit.
Definition: FLSHit.cxx:70
void SetTruePE(float aTruePE)
Set True PE value.
Definition: PCHit.h:89
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...
Vertical planes which measure X.
Definition: PlaneGeo.h:28
std::vector< double > BestPathEstimates(std::vector< geo::OfflineChan > const &chans, std::map< unsigned int, TVector3 > const &directions, std::vector< EPathQuality > &quals, std::vector< geo::OfflineChan > const &extras)
Definition: CalibUtil.cxx:98
Definition: event.h:19
double DistanceToEdgeInMC(TVector3 vertex, int &wall)
A collection of associated CellHits.
Definition: Cluster.h:47
double DetLength() const
void ProcessTrack(const rb::Track *track, daqchannelmap::DAQChannelMap *cmap, std::vector< caldp::PCHit > *pcHitsXY, std::vector< caldp::PCHit > *pcHitsZ, std::vector< caldp::PCHit > *pcHitsAvg) const
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void ProcessTrackForFLSHits(const rb::Track *track2, daqchannelmap::DAQChannelMap *cmap, std::vector< caldp::PCHit > *pcHitsXY) const
TString ip
Definition: loadincs.C:5
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
const PlaneGeo * Plane(unsigned int i) const
MaybeLogger_< ELseverityLevel::ELsev_error, false > LogError
geo::OfflineChan OfflineChan() const
Definition: CellHit.h:49
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
float GetExitX() const
Exit point of the particle (position, time and energy)
Definition: FLSHit.h:54
double fMaxEndDistFromEdge
Only keep tracks whose final point is not more than this distance from the edge.
lchan encodeLChan(int detId, plane_t plane, cell_t cell) const
double fBadStepSizeLimit
Only keep tracks without steps outside this limit, see GoodSteps comments.
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void SetTrueW(float aTrueW)
Set True W value.
Definition: PCHit.h:93
"Pre-calibration hit". Common input to calibration procedures
Definition: PCHit.h:16
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
art::ServiceHandle< geo::Geometry > fGeom
void SetView(geo::View_t view)
Definition: CellHit.h:54
void SetPlane(unsigned short plane)
Definition: CellHit.h:53
Provides efficient lookup of CellHits by plane and cell number.
Definition: HitMap.h:22
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
double fCompletenessCut
Only keep tracks with more than this fraction of slice hits included in each view.
void SetCell(unsigned short cell)
Definition: CellHit.h:52
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
Far Detector at Ash River, MN.
CDPStorage service.
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
void SetView(geo::View_t aView)
Set view value.
Definition: PCHit.h:77
double dy[NP][NC]
void SetReadDist(float aDist)
Set distance to cell readout.
Definition: PCHit.h:101
double dx[NP][NC]
void SetPath(float aPath)
Set Path value.
Definition: PCHit.h:83
void SetPE(float pe)
Definition: CellHit.h:55
double DistToEdgeZ(TVector3 vertex)
static DAQChannelMap * getInstance(int detID)
double dz[NP][NC]
double TrueW(const sim::FLSHit &fls) const
std::void_t< T > n
void SetGoodTime(bool aGoodT)
Set quality of timing fit.
Definition: PCHit.h:97
float GetXAverage(const int step) const
Get X-average for the step. This is in local coordinates.
Definition: FLSHit.h:76
T get(std::string const &key) const
Definition: ParameterSet.h:231
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
void reconfigure(fhicl::ParameterSet const &p)
float PE() const
Definition: CellHit.h:42
void GetXYZD(geo::OfflineChan chan, double w, double *xyzd)
Return position in world coordninates and distance to the readout.
Definition: CalibUtil.cxx:294
void SetDCM(int aDCM)
Set dcm value.
Definition: PCHit.h:71
void SetDiblock(int aDiblock)
Set diblock value.
Definition: PCHit.h:69
double DetHalfHeight() const
double median(TH1D *h)
Definition: absCal.cxx:524
art::ServiceHandle< geo::LiveGeometry > fLiveGeom
unsigned short Plane() const
Definition: OfflineChan.h:31
virtual double DistanceFromStart(double z) const
Definition: Track.cxx:229
void FindBelowThresholdCalibCandidates(const std::vector< geo::OfflineChan > &trkchans, const std::map< unsigned int, TVector3 > &directions, std::vector< EPathQuality > &quals, std::vector< double > &paths, std::vector< geo::OfflineChan > &chans)
Find empty cells that would have contributed to the calibration if they&#39;d been hit. We assume they fell below threshold.
Definition: CalibUtil.cxx:232
unsigned int NYCell() const
Number of cells in the y-view.
Definition: Cluster.h:108
Histograms used by attenuation calibration.
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
const std::string path
Definition: plot_BEN.C:43
float GetEntryZ() const
Definition: FLSHit.h:50
std::map< unsigned int, TVector3 > PlaneDirMap() const
map of the direction cosines at each plane (ie z)
Definition: Track.cxx:496
void produce(art::Event &e) override
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
pixel_t getPixel(dchan daqchan) const
Decode the pixel id from a dchan.
const std::vector< sim::PhotonSignal > HitToPhotonSignal(const art::Ptr< rb::CellHit > &hit) const
Returns the PhotonSignals contributing the signal in the specified hit.
const std::vector< sim::PhotonSignal > CellToPhotonSignal(unsigned int const &plane, unsigned int const &cell) const
Returns the PhotonSignals contributing the signal in the specified cell. WARNING: Use with extreme ca...
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
unsigned short Cell() const
Definition: OfflineChan.h:32
bool IsGoodTrack(rb::Track const &track, rb::Cluster const &slice)
void ProcessTrackTrajectory(const rb::Track *track, daqchannelmap::DAQChannelMap *cmap, std::set< double > const &zPlaneBounds, std::vector< caldp::PCHit > *pcHitsTraj) const
double DetHalfWidth() const
void SetPlane(int aPlane)
Set plane value.
Definition: PCHit.h:65
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
art::ServiceHandle< ds::DetectorService > fDS
A (plane, cell) pair.
Definition: OfflineChan.h:17
unsigned int NXCell() const
Number of cells in the x-view.
Definition: Cluster.h:106
bool IsPointLive(TVector3 vertex)
Note the muon catcher is considered bad; use in combination with InMuonCatcher if needed...
double MeanTNS(rb::AveragingScheme scheme=kDefaultScheme) const
Definition: Cluster.cxx:554
bool HaveTruthInfo() const
Is this a file with truth info in? (Is BackTracker going to be any use to you?)
Definition: BackTracker.h:137
#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
double fMaxVtxDistFromEdge
Only keep tracks with an initial point within this dist (cm) from the edge.
double fPathCut
Only keep hits with pathlength less than this.
exit(0)
double fDirZCut
Only keep tracks shallower than this.
cmap::CMap class source code
Definition: CMap.cxx:17
void SetW(float aW)
Set W value.
Definition: PCHit.h:85
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
dcm_id_t getDCM(dchan daqchan) const
Decode the dcm ID from a dchan.
dchan encodeDChan(int detID, diblock_t diblock, dcm_id_t dcm, feb_t feb, pixel_t pixel) const
float GetEdep() const
Get total Energy deposited into the cell for the whole FLSHit.
Definition: FLSHit.h:31
bool GoodTiming() const
Definition: CellHit.h:48
float GetExitY() const
Definition: FLSHit.h:55
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
#define MF_LOG_WARNING(category)
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:291
void FindZBoundaries(std::set< double > &planeZBounds)
Find the boundaries in the z direction of planes in the detector.
Definition: CalibUtil.cxx:332
bool InMuonCatcher(TVector3 vertex)
double fMinPath
Only keep hits with pathlength greater than this.
float GetYAverage(const int step) const
Get Y-average for the step. This is in local coordinates.
Definition: FLSHit.h:78
Float_t e
Definition: plot.C:35
T const * get() const
Definition: Ptr.h:149
uint32_t dchan
< DAQ Channel Map Package
float GetExitZ() const
Definition: FLSHit.h:56
ProductID put(std::unique_ptr< PROD > &&edp, FullSemantic< Level::Run > const semantic)
Definition: DataViewImpl.h:730
void SetTruePath(float aTruePath)
Set True path length value.
Definition: PCHit.h:91
virtual void Clear()
Forget about all owned cell hits.
Definition: Cluster.cxx:279
Float_t w
Definition: plot.C:20
diblock_t getDiBlock(dchan daqchan) const
Decode the diblock ID from a dchan.
Encapsulate the geometry of one entire detector (near, far, ndos)
double fMaxCellsPerPlane
Only keep tracks with fewer <cells/plane> than this limit.
void ProcessTrackForBelowThresholdHits(const rb::Track *track, daqchannelmap::DAQChannelMap *cmap, std::vector< caldp::PCHit > *pcHitsXYThresh) const
art::ServiceHandle< cheat::BackTracker > fBT
void SetCell(int aCell)
Set cell value.
Definition: PCHit.h:67
bool GoodSteps(rb::Track const &track)
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.
feb_t getFEB(dchan daqchan) const
Decode the feb id from a dchan.
float GetEntryY() const
Definition: FLSHit.h:49
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
enum BeamMode string
void SetTrueMeV(float aTrueMeV)
Set True energy (MeV) value.
Definition: PCHit.h:87