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