Slicer_module.cc
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////
2 // \file Slicer.cxx
3 // \brief TODO
4 // \version $Id: Slicer.cxx,v 1.44 2012-11-09 17:23:01 jpaley Exp $
5 // \authors rbpatter@caltech.edu, bckhouse@caltech.edu
6 ///////////////////////////////////////////////////////////////////////////
7 
13 
15 #include "Geometry/Geometry.h"
16 #include "NovaDAQConventions/DAQConventions.h"
17 #include "RecoBase/CellHit.h"
18 #include "RecoBase/Cluster.h"
19 
20 #include <cmath>
21 #include <string>
22 
23 
24 
25 // Slicer header
26 namespace slicer {
27 
28  class Slicer : public art::EDProducer {
29  public:
30  explicit Slicer(fhicl::ParameterSet const &pset);
31  virtual ~Slicer();
32 
33  void produce(art::Event& evt);
34 
35  protected:
36  // internal methods
37 
38  void NeighborSlice(std::vector<art::Ptr<rb::CellHit> > hitlist,
39  std::unique_ptr< std::vector<rb::Cluster> >& clustercol);
40 
41  void WindowSlice(std::vector<art::Ptr<rb::CellHit> > hitlist,
42  std::unique_ptr< std::vector<rb::Cluster> >& clustercol);
43 
44  // I might want to use RecoHit::T() or CellHit::TNS(). This method
45  // will facilitate that choice. For now, just CellHits...
46  float HitTime(art::Ptr<rb::CellHit> hit) { return hit->TNS(); }
47 
48  // Helper for NeighborSlice
49  bool AreNeighbors(art::Ptr<rb::CellHit> hit_i, art::Ptr<rb::CellHit> hit_j, bool weak = false);
50 
51  // Helper for WindowSlice
52  std::vector<std::vector<art::Ptr<rb::CellHit> > >
53  WindowSliceZ(std::vector<art::Ptr<rb::CellHit> >& hits, bool returnAllSlices = false) const;
54 
55  // Helper for WindowSlice, to make BlockSlices
56  std::vector<std::vector<art::Ptr<rb::CellHit> > >
57  BlockSlice(std::vector<art::Ptr<rb::CellHit> >& hits) const;
58 
59  // configuration settings
60  std::string fAlgorithm; ///< Neighbor/Window
61  bool fSliceByBlock; ///< If true, slices per block are written
62  std::string fCellHitInput; ///< Read CellHits from Cal/[fCellHitFolderIn]
63  bool fDoNotSliceAnything; ///< If true, one big slice is written out
64  unsigned int fMinHits; ///< The minimum number of hits needed to create a slice (this is in any view)
65 
66  // Configuration specific to the Window slicer
67  double fTimeWindow; // ns
68  bool fSeparateZGaps; ///< If true, the neighbor algorithm will separate zGaps
69  double fZGapAllowed; // planes
70  };
71 
72 }
73 
74 
75 
76 
77 
78 namespace slicer
79 {
80  //----------------------------------------------------------------------
82  fAlgorithm (pset.get< std::string >("Algorithm") ),
83  fSliceByBlock (pset.get< bool >("SliceByBlock") ),
84  fCellHitInput (pset.get< std::string >("CellHitInput") ),
85  fDoNotSliceAnything(pset.get< bool >("DoNotSliceAnything")),
86  fMinHits (pset.get< int >("MinHits") ),
87  fTimeWindow (pset.get< double >("TimeWindow") ),
88  fSeparateZGaps (pset.get< bool >("SeparateZGaps") ),
89  fZGapAllowed (pset.get< int >("ZGapAllowed") )
90  {
91  assert(fAlgorithm == "Neighbor" || fAlgorithm == "Window");
92 
93  produces< std::vector<rb::Cluster> >();
94  }
95 
96  //----------------------------------------------------------------------
98  {
99  }
100 
101  //----------------------------------------------------------------------
103  {
104  // Load hit list (CellHits)
106  evt.getByLabel(fCellHitInput, hitcol);
107 
108  std::vector<art::Ptr<rb::CellHit > > hitlist;
109  for(unsigned int i = 0; i < hitcol->size(); ++i){
110  art::Ptr<rb::CellHit> hit(hitcol, i);
111  hitlist.push_back(hit);
112  }
113 
114  std::unique_ptr< std::vector<rb::Cluster> > clustercol(new std::vector<rb::Cluster>);
115 
116  if(fAlgorithm == "Neighbor") NeighborSlice(hitlist, clustercol);
117  if(fAlgorithm == "Window") WindowSlice(hitlist, clustercol);
118 
119  evt.put(std::move(clustercol));
120  }
121 
122  //----------------------------------------------------------------------
123  void Slicer::
124  NeighborSlice(std::vector<art::Ptr<rb::CellHit> > hitlist,
125  std::unique_ptr< std::vector<rb::Cluster> >& clustercol)
126  {
127  // Slicer settings (some detector-specific). Hard-coding these
128  // until I have a better feel for the use-cases and the required
129  // user interface(s).
130 
131  int maxNeighbors = 1;
132  double bigGap = 0.; // ns
133  double lonerGap = 0.; // ns
134  // So that hits with the same exact time make it into a slice, independent
135  // of what order they happen to be sorted in
136  double graceGap = 1; // ns
137  int minNumWeakNeighbors = 0; // Don't require anything
138 
139  //get the geometry
141 
142  // figure out which detector we are using
143  int det = geo->DetId();
144 
145  switch (det) {
147  bigGap = 40;
148  lonerGap = 100;
149  //bigGap = 530; // With current DCS timing quality, we can only trust gaps longer
150  //lonerGap = 530; // than 500 ns (so, using something between 500 and 500+62.5)
151  minNumWeakNeighbors = 0;
152  break;
153  case novadaq::cnv::kNDOS:
154  // bigGap = 80;
155  // lonerGap = 250;
156  bigGap = 530; // With current DCS timing quality, we can only trust gaps longer
157  lonerGap = 530; // than 500 ns (so, using something between 500 and 500+62.5)
158  // Empirically, have seen hits within one tick like this get missed. More
159  // likely in NDOS due to the large number of bad channels. And not too
160  // much danger of dragging the whole detector in, because it's smaller.
161  graceGap = 62.5+1;
162  minNumWeakNeighbors = 0;
163  break;
165  default: // use kFar setting as default
166  bigGap = 100;
167  lonerGap = 350;
168  //bigGap = 530; // With current DCS timing quality, we can only trust gaps longer
169  //lonerGap = 530; // than 500 ns (so, using something between 500 and 500+62.5)
170  minNumWeakNeighbors = 1; // Require at least one neighbor
171  break;
172  }
173 
174  if (fDoNotSliceAnything) {
175  // effectively turn off slicing
176  bigGap = 9E10;
177  }
178 
179 
180  // Time sort the CellHits in the vector
181  rb::SortByTime(hitlist);
182  const int nHits = hitlist.size();
183 
184  // -----------------------------------------------------------------------
185  // These lists are in step with the hit list, holding the chunk number
186  // each hit is associated with and the number of "neighbors" the hit has
187  std::vector<int> hit_chunkID;
188  std::vector<int> hit_numNeighbors;
189 
190 
191  // -----------------------------------------------------------------------
192  // Count neighbors for each hit. Hits with zero or one neighbor can't
193  // bridge the rough-cut time slices
194 
195  for (int i=0;i<nHits;i++) {
196  hit_numNeighbors.push_back(0);
197  }
198  for (int i=0;i<nHits;i++) {
199  art::Ptr<rb::CellHit> hit_i = hitlist[i];
200  for (int j=i+1;j<nHits;j++) {
201  art::Ptr<rb::CellHit> hit_j = hitlist[j];
202  if (AreNeighbors(hit_i, hit_j)) {
203  hit_numNeighbors[i]++;
204  hit_numNeighbors[j]++;
205  }
206  }
207  }
208 
209 
210  // -----------------------------------------------------------------------
211  // Look for significant gaps in time, skipping any loner hits. Assign
212  // a number (chunkID) to each independent set of hits.
213 
214  int currChunkID = -1; // will increment
215  double prevT = -9e9;
216 
217  if (fDoNotSliceAnything) currChunkID = 0;
218 
219  for (int i=0;i<nHits;i++) {
220  art::Ptr<rb::CellHit> hit_i = hitlist[i];
221 
222  if (prevT+bigGap<HitTime(hit_i)) {
223  // -- big gap ==> new chunk
224  currChunkID++;
225  hit_chunkID.push_back(currChunkID);
226  }
227  else {
228  // -- little gap ==> gotta see if it's a loner hit (or loner hits)...
229  // skip any loner hits that are present
230  int j=i;
231  while (j<nHits&&hit_numNeighbors[j]<=maxNeighbors) j++;
232 
233  // if not a loner, automatic bridging:
234  bool bridge_them = (j==i)||fDoNotSliceAnything;
235  if (!bridge_them) {
236  // does the next non-loner hit bridge the gap anyway?
237  if (j<nHits) {
238  art::Ptr<rb::CellHit> hit_j = hitlist[j];
239  if (!(prevT+bigGap<HitTime(hit_j))) {
240  bridge_them = true;
241  }
242  }
243  }
244  for (;i<=j&&i<nHits;i++) {
245  if (!bridge_them) {
246  // new chunk for each loner
247  currChunkID++;
248  }
249  hit_chunkID.push_back(currChunkID);
250  }
251  hit_i = hitlist[--i];
252  }
253 
254  prevT = HitTime(hit_i);
255  }
256  const int nChunks = currChunkID+1;
257 
258 
259  // -----------------------------------------------------------------------
260  // Now look for chunks with too few hits in them. Figure out if these
261  // hits can be associated with a preceeding (preferred) or a following
262  // chunk by seeing if they are neighbors with any of the hits in that
263  // chunk AND if they are close enough in time.
264 
265  // In general, the definition of "neighbor" could be different here.
266  // We'll have to tune this up as we learn more about our events.
267 
268  // Repeat this until nothing new gets saved. (As hits get saved,
269  // other hits may be saved in a second (third?) round. Shouldn't
270  // take more than one or two passes, and speed shouldn't be an issue
271  // here, so I'm coding this up as transparently as I can.
272 
273  bool savedAnyHits = true;
274 
275  // for counting how many hits each chunk has
276  unsigned int numHitsInChunk[nChunks];
277 
278  while (savedAnyHits) {
279 
280  savedAnyHits = false;
281 
282  // count how many hits each chunk has
283  for (int iC=0;iC<nChunks;iC++) {
284  numHitsInChunk[iC] = 0;
285  }
286  for (int i=0;i<nHits;i++) {
287  numHitsInChunk[hit_chunkID[i]]++;
288  }
289 
290  // go through trying to save the hits in the tiny chunks
291  for (int i=0;i<nHits;i++) {
292 
293  art::Ptr<rb::CellHit> hit_i = hitlist[i];
294 
295  if (numHitsInChunk[hit_chunkID[i]]<fMinHits) {
296  // too few hits in this chunk, so let's try to save the hit(s).
297 
298  // Have we done something with this hit yet?
299  bool saved = false;
300 
301  // first, try to attach this hit to a preceeding (good) chunk
302  // if that doesn't work, try to attach to a follow-on (good) chunk
303  for(int step = -1; step <= +1; step += 2){
304  // Either from this hit back to the start, or forwards to the end
305  for (int j=i+step;j>=0 && j<nHits;j+=step) {
306  // j hit isn't in a good chunk, no good to us
307  if (numHitsInChunk[hit_chunkID[j]]<fMinHits) continue;
308 
309  // ok, this j-th hit is in a good chunk
310  art::Ptr<rb::CellHit> hit_j = hitlist[j];
311 
312  const double dt = fabs(HitTime(hit_i)-HitTime(hit_j));
313 
314  // too far in time already; stop looking in this direction
315  if (dt>lonerGap) break;
316 
317  if (dt<graceGap || AreNeighbors(hit_i,hit_j)) {
318  // it's very close in time to the i-th hit or closeish in
319  // time and close in space. Save it!
320  hit_chunkID[i] = hit_chunkID[j];
321  saved = true;
322  savedAnyHits = true;
323  break;
324  }
325  } // j loop
326  if(saved) break;
327  } // step loop
328  } // if small chunk
329  } // i loop
330  } // while (savedAnyHits)
331 
332  // Note: as long as the loop above terminates when (!savedAnyHits),
333  // then numHitsInChunk[] should be up to date (for use below).
334 
335  // // RBP - diagnostics:
336  // for (int i=0;i<nHits;i++) {
337  // const rb::CellHit *hit_i = hitlist[i];
338  // cout << i << " " << hit_chunkID[i] << " "
339  // << hit_numNeighbors[i] << " " << hit_i->View() << " "
340  // << hit_i->Plane() << " " << hit_i->Cell() << " " << HitTime(hit_i) << endl;
341  // }
342  // cout << "----------------" << endl;
343 
344  // Count how many neighbours each hit has from its own chunk, but with
345  // weak conditions.
346  std::vector<int> numWeakNeighbors;
347  numWeakNeighbors.resize(nHits);
348 
349  for (int i=0;i<nHits-1;i++) {
350  currChunkID = hit_chunkID[i];
351  for (int j=i+1;j<nHits;j++) {
352  if(i == j) continue; // A hit doesn't count as its own neighbour
353  if(hit_chunkID[j] != currChunkID) continue; // Must be in same chunk
354 
355  if(AreNeighbors(hitlist[i], hitlist[j], true)){
356  numWeakNeighbors[i]++;
357  numWeakNeighbors[j]++;
358  }
359  }
360  }
361 
362  // Now we remove hits from chunks if they are really far away from all other hits
363  // i.e. no weak neighbors
364  // These hits are relegated to the noise slice
365  for (int i=0;i<nHits;i++)
366  {
367  if(numWeakNeighbors[i] < minNumWeakNeighbors)
368  {
369  currChunkID = hit_chunkID[i];
370  // Decrement the number of hits in that chunk because we are removing it
371  numHitsInChunk[currChunkID] -= 1;
372  // Move that hit to the noise slice
373  hit_chunkID[i] = nChunks;
374  }
375 
376  }
377 
378  // We're going to create vectors of cell hits for our time slices
379  // and try to divide ones with large z-gaps, if desired
380  std::vector<std::vector<art::Ptr<rb::CellHit> > >* tSlices
381  = new std::vector<std::vector<art::Ptr<rb::CellHit> > >(nChunks+1);
382  std::vector<std::vector<art::Ptr<rb::CellHit> > >* slices;
383 
384 
385  for (int i=0;i<nHits;i++) {
386  currChunkID = hit_chunkID[i];
387  art::Ptr<rb::CellHit> hit = hitlist[i];
388  if (currChunkID>=0&&numHitsInChunk[currChunkID]>=fMinHits) {
389  // push to a good time slice
390  (*tSlices)[currChunkID].push_back(hit);
391  }
392  else {
393  // stray hit; push to noise cluster
394  (*tSlices)[nChunks].push_back(hit);
395  }
396  }
397 
398 
399  // Split z-gaps if desired
400  if(fSeparateZGaps){
401  slices = new std::vector<std::vector<art::Ptr<rb::CellHit> > >;
402  // Loop over time slices to separate z-gaps, excluding the noise slice
403  for(unsigned int ts = 0; ts < tSlices->size() - 1; ++ts)
404  {
405  int nHitsTime = (*tSlices)[ts].size();
406  int nHitsZ = 0;
407  // call the z-slicing function with returnAllSlices set to true
408  std::vector<std::vector<art::Ptr<rb::CellHit> > > zSlices =
409  WindowSliceZ((*tSlices)[ts], true);
410 
411  for(unsigned int zs = 0; zs < zSlices.size(); ++zs)
412  {
413  // Order hits by time
414  rb::SortByTime(zSlices[zs]);
415  // Add to slices vector
416  slices->push_back(zSlices[zs]);
417  nHitsZ += zSlices[zs].size();
418 
419  }
420  assert(nHitsTime == nHitsZ);
421  }
422 
423  slices->push_back((*tSlices)[nChunks]);
424  delete tSlices;
425  }else{
426  // Slices are just the time slices
427  slices = tSlices;
428 
429  }
430 
431  // -----------------------------------------------------------------------
432  // For now, call it quits here. Any chunks with fewer than fMinHits hits
433  // will be banished to a "noise" cluster, to be written out last.
434  // Similarly for hits really far from all others in their chunk.
435 
436 
437  // cluster storage:
438  std::vector<rb::Cluster> clusters(slices->size());
439 
440 
441 
442 
443  // index of the noise slice
444  size_t ns = slices->size() - 1;
445  // Add all of the physics slices, the last in the vector is the noise slice
446  for(size_t s = 0; s < ns; ++s)
447  {
448 
449  if ((*slices)[s].size() >= fMinHits)
450  {
451  // add hits from this slice to the cluster
452  for(size_t h = 0; h < (*slices)[s].size(); ++h)
453  {
454  clusters[s].Add((*slices)[s][h]);
455  }
456  }
457  else
458  {
459  // slice is too small; move hits to noise slice
460  for(size_t h = 0; h < (*slices)[s].size(); ++h)
461  {
462 
463  (*slices)[ns].push_back((*slices)[s][h]);
464  }
465  }
466  }
467 
468  // order the noise slice by time
469  rb::SortByTime((*slices)[ns]);
470  // Add the hits in the noise slice to its cluster
471  for(size_t h = 0; h < (*slices)[ns].size(); ++h)
472  {
473  clusters[ns].Add((*slices)[ns][h]);
474  }
475 
476  // tag the noise cluster as such:
477  clusters[ns].SetNoise(true);
478 
479  if(slices) delete slices;
480  // We're going to do a sanity check to make sure slices contain all the hits
481  int nHitsCalHit = hitlist.size();
482  int nHitsSlicer = 0;
483 
484 
485  // Write the clusters.
486  // Always write the noise cluster (iC==ns), even if empty.
487  //iterate over the vector and remove any that fail to have more than
488  //the minimum hits
489  for (size_t iC=0;iC < clusters.size(); ++iC) {
490  if (iC==ns || clusters[iC].NCell()>=fMinHits) {
491  nHitsSlicer += clusters[iC].NCell();
492  clustercol->push_back(clusters[iC]);
493  }
494  else{
495  assert(clusters[iC].NCell() == 0 && " skipping a cluster with a non-zero size" );
496  }
497  }
498  assert(nHitsCalHit == nHitsSlicer);
499 
500 
501  }
502 
503  //----------------------------------------------------------------------
505 
506  // same view & small plane gap & small cell gap ==> neighbors
507 
508  const int allowedPlaneGap = weak ? 20 : 5;
509  const int allowedCellGap = weak ? 24 : 6;
510 
511  return
512  ( hit_i->View()==hit_j->View() &&
513  TMath::Abs(hit_i->Plane()-hit_j->Plane()) <= allowedPlaneGap &&
514  TMath::Abs(hit_i->Cell() -hit_j->Cell() ) <= allowedCellGap );
515 
516  }
517 
518  //----------------------------------------------------------------------
519  void Slicer::
521  std::unique_ptr<std::vector<rb::Cluster> >& slices)
522  {
523  // A timeslice is at most fTimeWindow long, and must contain at least
524  // fMinHits hits. Any fTimeWindow-sized subwindow drawn within it
525  // will also contain at least fMinHits hits.
526 
527  slices->clear();
528  if(hits.size() == 0) return;
529 
530  rb::Cluster noiseClust;
531  noiseClust.SetNoise(true);
532 
533  if(hits.empty()){
534  slices->push_back(noiseClust);
535  return;
536  }
537 
538  std::vector<std::vector<art::Ptr<rb::CellHit> > > tSlices;
539 
541 
542  typedef std::vector<art::Ptr<rb::CellHit> >::iterator it_t;
543  // The inclusive->exclusive range we're considering for the current slice
544  it_t itStart = hits.begin();
545  it_t itEnd = hits.begin();
546 
547  // We use this point to test if we can expand the slice
548  it_t itTest = hits.begin();
549 
550  while(true){
551  // Grow the slice by moving the end until it's fTimeWindow in size
552  while(HitTime(*itEnd) - HitTime(*itStart) < fTimeWindow){
553  ++itEnd;
554  // We ran out of hits altogether
555  if(itEnd == hits.end()) break;
556  }
557  // Break all the way out of the loop
558  if(itEnd == hits.end()) break;
559 
560  // The window doesn't have enough charge, drop first hit and try again
561  if(itEnd - itStart < fMinHits){
562  noiseClust.Add(*itStart);
563  ++itStart;
564  assert(itStart <= itEnd);
565  continue;
566  }
567 
568  // The window is long enough and has enough charge
569  // Catch itTest up until it defines the start of a fTimeWindow-sized
570  // window ending on itEnd
571  while(HitTime(*itEnd) - HitTime(*itTest) > fTimeWindow) ++itTest;
572 
573  // If growing the window more would cause the itTest-itEnd window
574  // not to have enough hits, then write out what we have.
575  if(itEnd+1 - itTest < fMinHits){
576  std::vector<art::Ptr<rb::CellHit> > slice(itStart, itEnd);
577  tSlices.push_back(slice);
578  itStart = itEnd;
579  continue;
580  }
581 
582  // We can still legally grow the window. Do that and try again
583  ++itEnd;
584  if(itEnd == hits.end()) break;
585  }
586 
587  // If there are enough hits in the final slice, make it one too
588  if(itEnd - itStart > fMinHits){
589  std::vector<art::Ptr<rb::CellHit> > slice(itStart, itEnd);
590  tSlices.push_back(slice);
591  }
592  else{
593  for(it_t it = itStart; it != itEnd; ++it) noiseClust.Add(*it);
594  }
595 
596  // Split every time slice up based on z if it has gaps in
597  for(unsigned int s = 0; s < tSlices.size(); ++s){
598  std::vector<std::vector<art::Ptr<rb::CellHit> > > zSlices =
599  WindowSliceZ(tSlices[s]);
600  for(unsigned int n = 0; n < zSlices.size(); ++n){
601  rb::Cluster clust;
602  for(unsigned int h = 0; h < zSlices[n].size(); ++h){
603  clust.Add(zSlices[n][h]);
604  }
605  slices->push_back(clust);
606  }
607  }
608 
609  // Noise cluster goes at the end
610  slices->push_back(noiseClust);
611  }
612 
613  ////////////////////////////////////////////////////////////////////////
614  std::vector<std::vector<art::Ptr<rb::CellHit> > > Slicer::
615  WindowSliceZ(std::vector<art::Ptr<rb::CellHit> >& hits, bool returnAllSlices) const
616  {
617  std::vector<std::vector<art::Ptr<rb::CellHit> > > zSlices;
618 
619  // Shouldn't happen, but still...
620  if(hits.size() == 0) return zSlices;
621 
623 
624  // The extremes, inclusive and exclusive of our candidate slice
625  typedef std::vector<art::Ptr<rb::CellHit> >::iterator it_t;
626  it_t itBegin = hits.begin();
627  it_t itEnd = hits.begin();
628 
629  while(true){
630  const int endPlane = (*itEnd)->Plane();
631  ++itEnd;
632 
633  if(itEnd == hits.end()) break;
634 
635  const int nextPlane = (*itEnd)->Plane();
636 
637  if(nextPlane - endPlane > fZGapAllowed){
638  std::vector<art::Ptr<rb::CellHit> > slice(itBegin, itEnd);
639  if(slice.size() > fMinHits || returnAllSlices) zSlices.push_back(slice);
640  // Start a new slice
641  itBegin = itEnd;
642  }
643  } // end while
644 
645  // Whatever hits are left, give them their own slice
646  if(itBegin != hits.end()){
647  std::vector<art::Ptr<rb::CellHit> > slice(itBegin, itEnd);
648  if(slice.size() > fMinHits|| returnAllSlices) zSlices.push_back(slice);
649  }
650 
651  std::vector<std::vector<art::Ptr<rb::CellHit> > > bSlices;
652  if(fSliceByBlock) {
653  // Split every z slice up based on block position
654  for(unsigned int s = 0; s < zSlices.size(); ++s){
655  bSlices = BlockSlice(zSlices[s]);
656  }
657  }
658  if(fSliceByBlock) return bSlices;
659  else return zSlices;
660  }
661 
662  ////////////////////////////////////////////////////////////////////////
663  std::vector<std::vector<art::Ptr<rb::CellHit> > > Slicer::
665  {
666  // Get the geometry
669 
670  // Figure out which detector we are using
671  // We need to know in order to grab the correct block id
672  const int det = geo->DetId();
673 
674  std::vector<std::vector<art::Ptr<rb::CellHit> > > bSlices;
675 
676  // Shouldn't happen, but still...
677  if(hits.size() == 0) return bSlices;
678 
679  // The extremes, inclusive and exclusive of our candidate slice
680  typedef std::vector<art::Ptr<rb::CellHit> >::iterator it_t;
681  it_t itBegin = hits.begin();
682  it_t itEnd = hits.begin();
683 
684  while(true){
685  daqchannelmap::lchan lchan = cmap->Map()->encodeLChan(det, (*itEnd)->Plane(),1);
686  const int block = cmap->Map()->computeBlock(lchan);
687  ++itEnd;
688 
689  // Check that we have reached end of hit list.
690  // If we have, break out of loop
691  if(itEnd == hits.end()) break;
692 
693  daqchannelmap::lchan lchan2 = cmap->Map()->encodeLChan(det, (*itEnd)->Plane(),1);
694  const int nextBlock = cmap->Map()->computeBlock(lchan2);
695 
696  if(block != nextBlock){
697  std::vector<art::Ptr<rb::CellHit> > slice(itBegin, itEnd);
698  bSlices.push_back(slice);
699  // Start a new slice
700  itBegin = itEnd;
701  }
702  } // end while
703 
704  // Save all of the hits that remain into a block slice
705  if(itBegin != hits.end()){
706  std::vector<art::Ptr<rb::CellHit> > slice(itBegin,itEnd);
707  bSlices.push_back(slice);
708  }
709  return bSlices;
710  }
711 
712 } // end namespace slicer
713 /////////////////////////////////////////////////////////////////////////
714 
716 namespace slicer
717 {
719 }
float TNS() const
Definition: CellHit.h:46
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void SetNoise(bool noise)
Declare the cluster to consist of noise hits or not.
Definition: Cluster.h:161
unsigned short Plane() const
Definition: CellHit.h:39
geo::View_t View() const
Definition: CellHit.h:41
A collection of associated CellHits.
Definition: Cluster.h:47
const daqchannelmap::DAQChannelMap * Map() const
Definition: CMap.h:57
void SortByTime(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in time order (earliest to latest).
Definition: CellHit.cxx:134
DEFINE_ART_MODULE(TestTMapFile)
void NeighborSlice(std::vector< art::Ptr< rb::CellHit > > hitlist, std::unique_ptr< std::vector< rb::Cluster > > &clustercol)
lchan encodeLChan(int detId, plane_t plane, cell_t cell) const
std::vector< std::vector< art::Ptr< rb::CellHit > > > WindowSliceZ(std::vector< art::Ptr< rb::CellHit > > &hits, bool returnAllSlices=false) const
bool fSliceByBlock
If true, slices per block are written.
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
bool fSeparateZGaps
If true, the neighbor algorithm will separate zGaps.
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
unsigned short Cell() const
Definition: CellHit.h:40
const XML_Char * s
Definition: expat.h:262
Definition: Cand.cxx:23
Far Detector at Ash River, MN.
block
print "ROW IS " print row
Definition: elec2geo.py:31
void hits()
Definition: readHits.C:15
bool AreNeighbors(art::Ptr< rb::CellHit > hit_i, art::Ptr< rb::CellHit > hit_j, bool weak=false)
Prototype Near Detector on the surface at FNAL.
unsigned int fMinHits
The minimum number of hits needed to create a slice (this is in any view)
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
bool fDoNotSliceAnything
If true, one big slice is written out.
int evt
Near Detector in the NuMI cavern.
const double j
Definition: BetheBloch.cxx:29
static const double ns
Module that plots metrics from reconstructed cosmic ray data.
float HitTime(art::Ptr< rb::CellHit > hit)
virtual ~Slicer()
void produce(art::Event &evt)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
std::string fAlgorithm
Neighbor/Window.
Definition: structs.h:12
assert(nhit_max >=nhit_nbins)
Class to help Slicer4D manage the collection of hits.
std::string fCellHitInput
Read CellHits from Cal/[fCellHitFolderIn].
std::vector< std::vector< art::Ptr< rb::CellHit > > > BlockSlice(std::vector< art::Ptr< rb::CellHit > > &hits) const
virtual block_t computeBlock(lchan logicalchan) const =0
Which block is this lchan in?
void WindowSlice(std::vector< art::Ptr< rb::CellHit > > hitlist, std::unique_ptr< std::vector< rb::Cluster > > &clustercol)
Helper for AttenCurve.
Definition: Path.h:10
Slicer(fhicl::ParameterSet const &pset)
Encapsulate the geometry of one entire detector (near, far, ndos)
void SortByPlane(std::vector< art::Ptr< rb::CellHit > > &c)
Sort c in plane order (low to high).
Definition: CellHit.cxx:124
enum BeamMode string