DiscreteTracker_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file DiscreteTracker.cxx
3 // \author Christopher Backhouse - bckhouse@caltech.edu
4 ////////////////////////////////////////////////////////////////////////
5 
6 // ROOT includes
7 #include "TH1.h"
8 
9 // NOvASoft includes
10 #include "GeometryObjects/Geo.h"
11 
12 #include "RecoBase/Track.h"
13 #include "RecoBase/Vertex.h"
14 #include "RecoBase/Cluster.h"
15 #include "RecoBase/CellHit.h"
16 #include "RecoBase/HitMap.h"
17 #include "RecoBase/FilterList.h"
18 
19 
21 #include "Utilities/AssociationUtil.h"
22 #include "Utilities/func/MathUtil.h"
23 
24 // Package includes
25 #include "DiscreteTracker/Cand.h"
26 #include "DiscreteTracker/Chain.h"
27 #include "DiscreteTracker/Chunk.h"
28 #include "DiscreteTracker/Types.h"
29 #include "DiscreteTracker/View.h"
30 
31 // Framework includes
37 #include "fhiclcpp/ParameterSet.h"
38 
39 namespace dt
40 {
41  /// Job module performing "discrete" track reconstruction
43  {
44  public:
45  explicit DiscreteTracker(const fhicl::ParameterSet& pset);
47 
48  /// Splits hits according to slice and view for feeding to \ref
49  /// RecoDiscrete. Writes the returned tracks into the event.
50  void produce(art::Event& evt);
51 
52  void reconfigure(const fhicl::ParameterSet& pset);
53 
54  void beginJob();
55  void endJob();
56 
57  protected:
58  /// \brief Cut to prevent choking on large air shower events
59  ///
60  /// Adapted from the version in KalmanTrack
61  bool HighActivity(const art::Ptr<rb::Cluster>& clust) const;
62 
64  const Cand& prev,
65  const View& v,
66  const dt::View::ChunkMap& chunk_map,
67  Direction dir,
68  bool exactHits) const;
69 
70  void ExtendChain(const View& v,
71  const dt::View::ChunkMap& chunk_map,
72  Chain& chain, Direction dir) const;
73 
74  Chain RecoOneTrack(const View& v,
75  const dt::View::ChunkMap& chunk_map,
76  std::vector<rb::Track>* seeds_out,
77  bool& ok) const;
78 
79  void SpliceChains(std::vector<Chain>& chains) const;
80 
81  bool SpliceIsSane(const Chain& a, const Chain& b) const;
82 
83  std::vector<rb::Track>
85  std::vector<rb::Track>* seeds_out,
86  std::vector<rb::Track>* unspliced_out,
87  std::vector<rb::Track>* chains_out) const;
88 
89  std::vector<Segment>
90  PossibleNextSegs(const View& v,
91  const dt::View::ChunkMap& chunk_map,
92  const Cand& cand, Direction dir,
93  std::vector<Chunk>& chunks_out,
94  bool shw, bool exactHits) const;
95 
96  double DotScoreHelper(double m1, double m2) const;
97 
98  double DotScore(const Cand& c1, const Cand& c2,
99  bool allowFlip = false) const;
100 
101  bool TryExtendCandOne(Cand& cand,
102  const View& v,
103  const dt::View::ChunkMap& chunk_map,
104  Direction dir,
105  const Cand& prev_dir,
106  bool shw,
107  bool exactHits) const;
108 
109  bool TryExtendCand(Cand& cand,
110  const Cand& old_cand,
111  const View& v,
112  const dt::View::ChunkMap& chunk_map,
113  Direction dir,
114  bool exactHits) const;
115 
116  Cand BestSeed(const dt::View& v, bool& ok) const;
117 
118  // Modifies \a chain in-place, returns bits clipped off, if any
119  std::vector<Chain> CheckAndFixSparseChain(Chain& chain) const;
120 
121  /// \brief Toy vertex finder
122  ///
123  /// Uses the intersection of the two tracks in each view whose intersection
124  /// is proportionally nearest their ends.
125  rb::Vertex* FindVertex(std::vector<rb::Track> xTrks,
126  std::vector<rb::Track> yTrks) const;
127 
128  void MergeParasiteTracks(std::vector<rb::Track>& tracks) const;
129 
130  std::string fClusterInput; ///< Module name to look for Clusters under
131  bool fConsiderBadChans; ///< Passed to \ref View
132  bool fObeyPreselection; ///< Check rb::IsFiltered?
133 
134  /// \brief Joined Cands that only have one chunk in common must have a
135  /// dot-product better than this.
137  /// \brief The number of adjacent cells in a single plane, otherwise not
138  /// fit into any track, required to declare them a vertical track
140 
141  bool fOnlyOneTrack; ///< Stop after the first track in each view
142 
143  TH1* fAngle; ///< Opening angle of all cands considered for joining
144  TH1* fSinglesAngle; ///< For all cands with only one overlapping chunk
145 
146  FILE* fDebugFile; ///< File that track details are written to
147  };
148 
149 
150  //......................................................................
152  : fDebugFile(0)
153  {
154  produces<std::vector<rb::Track> >();
155  produces<std::vector<rb::Track> >("seeds");
156  produces<std::vector<rb::Track> >("unspliced");
157  produces<std::vector<rb::Track> >("chains");
158  produces<std::vector<rb::Vertex> >();
159  // Tracks are associated with the slice they came from
160  produces<art::Assns<rb::Track, rb::Cluster> >();
161 
162  reconfigure(pset);
163  }
164 
165  //......................................................................
167  {
168  }
169 
170  //......................................................................
172  {
173  fClusterInput = pset.get<std::string>("ClusterInput");
174  fConsiderBadChans = pset.get<bool>("ConsiderBadChans");
175  fMinDotProduct = pset.get<double>("MinDotProduct");
176  fVerticalTrackMinCells = pset.get<int>("VerticalTrackMinCells");
177  fOnlyOneTrack = pset.get<bool>("OnlyOneTrack");
178  fObeyPreselection = pset.get<bool>("ObeyPreselection");
179 
181  fDebugFile = fopen(pset.get<std::string>("DebugFile").c_str(), "w");
182  }
183 
184  //......................................................................
186  {
188  fAngle = tfs->make<TH1F>("angle", ";#Delta#theta", 100, 0, 1);
189  fSinglesAngle = tfs->make<TH1F>("singles_angle", ";#Delta#theta",
190  100, 0, 1);
191  }
192 
193  //......................................................................
195  {
197  }
198 
199  //......................................................................
201  {
202  std::unique_ptr<std::vector<rb::Track> > debug_seeds_col(new std::vector<rb::Track>);
203  std::unique_ptr<std::vector<rb::Track> > debug_unspliced_col(new std::vector<rb::Track>);
204  std::unique_ptr<std::vector<rb::Track> > debug_chains_col(new std::vector<rb::Track>);
205  std::unique_ptr<std::vector<rb::Track> > trackcol(new std::vector<rb::Track>);
206  std::unique_ptr<std::vector<rb::Vertex> > vtxcol(new std::vector<rb::Vertex>);
207  // Associate tracks with the slice they came from
208  std::unique_ptr<art::Assns<rb::Track, rb::Cluster> > assns(new art::Assns<rb::Track, rb::Cluster>);
209 
211  evt.getByLabel(fClusterInput, clusts);
212 
213  const int clustMax = clusts->size();
214  for(int clustIdx = 0; clustIdx < clustMax; ++clustIdx){
215  art::Ptr<rb::Cluster> clust(clusts, clustIdx);
216  if(fObeyPreselection && rb::IsFiltered(evt, clust)) continue;
217 
218  if(clust->IsNoise()) continue;
219 
220  if(HighActivity(clust)){
221  mf::LogWarning("DiscreteTracker") << "Skipping event due to very high activity";
222  continue;
223  }
224 
225  std::vector<rb::Track> tracks[2];
226  for(geo::View_t view: {geo::kX, geo::kY}){
227  const char viewC = (view == geo::kX) ? 'x' : 'y';
228 
230  const int Nx = clust->NCell(view);
231  for(int n = 0; n < Nx; ++n) cells.push_back(clust->Cell(view, n));
232 
233  fprintf(fDebugFile, "%d: %c: ", evt.event(), viewC);
234 
235  tracks[view] = RecoDiscrete(cells,
236  debug_seeds_col.get(),
237  debug_unspliced_col.get(),
238  debug_chains_col.get());
239 
240  for(unsigned int n = 0; n < tracks[view].size(); ++n)
241  fprintf(fDebugFile, "%u ", tracks[view][n].NCell());
242  fprintf(fDebugFile, "\n");
243 
244  for(std::vector<rb::Track>::const_iterator it = tracks[view].begin();
245  it != tracks[view].end(); ++it){
246  trackcol->push_back(*it);
247  // Associate most recently added track with the slice
248  util::CreateAssn(*this, evt, *trackcol, clust, *assns);
249  } // end for it
250  } // end for view
251 
252  const rb::Vertex* vtx = FindVertex(tracks[geo::kX], tracks[geo::kY]);
253  if(vtx) vtxcol->push_back(*vtx);
254  delete vtx;
255  } // end for clustIdx
256 
257  evt.put(std::move(debug_seeds_col), "seeds");
258  evt.put(std::move(debug_unspliced_col), "unspliced");
259  evt.put(std::move(debug_chains_col), "chains");
260  evt.put(std::move(trackcol));
261  evt.put(std::move(vtxcol));
262  evt.put(std::move(assns));
263  }
264 
265  //......................................................................
267  {
270 
271  // get the total number of non masked off channels by looping over all
272  // planes/cells in the detector
273  int numGoodChannels = 0;
274 
275  for(unsigned int iPlane = 0; iPlane < geom->NPlanes(); ++iPlane){
276  const geo::PlaneGeo* p = geom->Plane(iPlane);
277  for(unsigned int iCell = 0; iCell < p->Ncells(); ++iCell){
278  if(!bcl->IsBad(iPlane, iCell)) ++numGoodChannels;
279  }
280  }
281 
282  return clust->NCell() > 0.15*numGoodChannels;
283  }
284 
285  //......................................................................
286  std::vector<Segment>
288  const dt::View::ChunkMap& chunk_map,
289  const Cand& cand, Direction dir,
290  std::vector<Chunk>& chunks_out,
291  bool shw, bool exactHits) const
292  {
293  std::vector<Segment> ret;
294 
295  // If the current chunk is only partially taken then the next thing to do
296  // must be to try and complete it.
297  if(dir == kUpstream && !cand.FirstSeg().left){
298  ret.push_back(cand.FirstChunk().GetSegs().first);
299  return ret;
300  }
301 
302  if(dir == kDownstream && cand.LastSeg().left){
303  ret.push_back(cand.LastChunk().GetSegs().second);
304  return ret;
305  }
306 
307  int nextPlane = v.AdjacentPlane(cand.ExtremalPlane(dir), dir);
308  if(nextPlane == -1){
309  if(!exactHits) return ret;
310  if(dir == kDownstream){
311  nextPlane = chunk_map.begin()->first;
312  if(nextPlane <= int(cand.LastPlane())) return ret;
313  }
314  else{
315  nextPlane = chunk_map.rbegin()->first;
316  if(nextPlane >= int(cand.FirstPlane())) return ret;
317  }
318  }
319 
320  dt::View::ChunkMap::const_iterator it = chunk_map.find(nextPlane);
321  if(exactHits){
322  while(it == chunk_map.end()){
323  nextPlane = v.AdjacentPlane(nextPlane, dir);
324  if(nextPlane == -1) return ret;
325  it = chunk_map.find(nextPlane);
326  }
327  }
328  if(it == chunk_map.end()) return ret;
329 
330  const std::vector<Chunk>& chunks = it->second;
331  const unsigned int chunkMax = chunks.size();
332 
333  for(unsigned int chunkIdx = 0; chunkIdx < chunkMax; ++chunkIdx){
334  Chunk chunk = chunks[chunkIdx];
335  if(exactHits) chunk.SetUp(cand.Up());
336 
337  if(chunk.Up() != cand.Up()) continue;
338  if(shw && !chunk.IsWorthShowering()) continue;
339  if(shw) chunk.SetShowerChunk();
340 
341  chunks_out.push_back(chunk);
342  if(dir == kUpstream)
343  ret.push_back(chunk.GetSegs().second);
344  else
345  ret.push_back(chunk.GetSegs().first);
346  }
347 
348  return ret;
349  }
350 
351  //......................................................................
352  double DiscreteTracker::DotScoreHelper(double m1, double m2) const
353  {
354  return fabs(1+m1*m2)/sqrt((1+m1*m1)*(1+m2*m2));
355  }
356 
357  //......................................................................
358  double DiscreteTracker::DotScore(const Cand& c1, const Cand& c2,
359  bool allowFlip) const
360  {
361  // Can't even determine a direction
362  if(c1.NSegs() < 2 || c2.NSegs() < 2) return -1;
363 
364  double bestDot = -2;
365 
366  for(bool flip1: {false, true}){
367  if(flip1 && !allowFlip) continue;
368 
369  bool ok = true;
370  const Cand c1f = flip1 ? c1.MaybeFlip(ok) : c1;
371  if(!ok) continue;
372 
373  for(bool flip2: {false, true}){
374  if(flip2 && !allowFlip) continue;
375 
376  bool ok = true;
377  const Cand c2f = flip2 ? c2.MaybeFlip(ok) : c2;
378  if(!ok) continue;
379 
380  const double m1lo = c1f.MinGradient();
381  const double m1hi = c1f.MaxGradient();
382  const double m2lo = c2f.MinGradient();
383  const double m2hi = c2f.MaxGradient();
384 
385  assert(m1lo <= m1hi);
386  assert(m2lo <= m2hi);
387 
388  // Overlapping ideas of direction, best is perfect
389  if((m1hi > m2lo && m1lo < m2hi) ||
390  (m1lo < m2hi && m1hi > m2lo)) return 1;
391 
392  bestDot = std::max(bestDot, DotScoreHelper(m1hi, m2lo));
393  bestDot = std::max(bestDot, DotScoreHelper(m1lo, m2hi));
394  } // end for flip2
395  } // end for flip1
396 
397  return bestDot;
398  }
399 
400  //......................................................................
402  const View& v,
403  const dt::View::ChunkMap& chunk_map,
404  Direction dir,
405  const Cand& prev_dir,
406  bool shw, bool exactHits) const
407  {
408  // TODO: hoist this out of while loop below
409  std::vector<Chunk> chunks;
410  const std::vector<Segment> segs = PossibleNextSegs(v, chunk_map,
411  cand, dir, chunks, shw,
412  exactHits);
413  if(segs.empty()) return false;
414 
415  // Backup of cand, so we can restore it in case TryAddSeg succeeds
416  Cand bak = cand;
417 
418  std::vector<Cand> extended;
419 
420  const unsigned int segMax = segs.size();
421  for(unsigned int segIdx = 0; segIdx < segMax; ++segIdx){
422  const Segment& seg = segs[segIdx];
423 
424  if(cand.TryAddSeg(seg, dir)){
425 
426  // May need to add the chunk to match
427  if((dir == kUpstream && !seg.left) ||
428  (dir == kDownstream && seg.left)){
429  cand.AddChunk(chunks[segIdx], dir);
430  }
431 
432  const bool smallOverlap = (bak.NSegs() <= 2);
433  if(true){//smallOverlap){
434  // If there's only a small overlap, ensure dot product is sensible
435  const double dotscore = DotScore(prev_dir, cand);
436  fAngle->Fill(acos(dotscore));
437  if(smallOverlap) fSinglesAngle->Fill(acos(dotscore));
438  if(smallOverlap && dotscore < fMinDotProduct){
439  cand = bak;
440  continue;
441  }
442  }
443 
444  // Otherwise it's good
445  extended.push_back(cand);
446  cand = bak;
447  } // end if added seg
448  } // end for segIdx
449 
450  if(extended.empty()) return false;
451 
452  // TODO: Arbitrary for now
453  cand = extended.front();
454  if(extended.size() > 1){
455  // std::cout << "***WARNING: MORE THAN ONE EXTENSION CANDIDATE (" << extended.size() << "). SHOULD HAVE TRIED TO MAKE A SMARTER DECISION***" << std::endl;
456  }
457 
458  return true;
459  }
460 
461  //......................................................................
463  const Cand& old_cand,
464  const View& v,
465  const dt::View::ChunkMap& chunk_map,
466  Direction dir,
467  bool exactHits) const
468  {
469  const Cand bak = cand;
470 
471  for(bool shwNext: {false, true}){
472  if(shwNext && !cand.IsExtremalChunkComplete(dir)) continue;
473  for(bool shwLast: {false, true}){
474  if(shwLast && cand.ExtremalChunk(dir).IsShowerChunk()) continue;
475  if(shwLast && !cand.ExtremalChunk(dir).IsWorthShowering()) continue;
476 
477  if(shwLast && cand.NSegs() <= 2) continue;
478 
479  if(shwLast) cand.MarkExtremalChunkShower(dir);
480  // Reset back by hand at end of loop. Careful about early exits
481 
482  bool any = false;
483  while(TryExtendCandOne(cand, v, chunk_map, dir, old_cand, shwNext, exactHits)) any = true;
484 
485  if(any) return true;
486 
487  // TODO: can probably hoist this higher with care
488  bool ok;
489  Cand flipped = cand.MaybeFlip(ok);
490  if(!ok){
491  if(shwLast) cand = bak;
492  continue;
493  }
494 
495  while(TryExtendCandOne(flipped, v, chunk_map, dir, old_cand, shwNext, exactHits)) any = true;
496 
497  if(any){
498  cand = flipped;
499  return true;
500  }
501 
502  if(shwLast) cand = bak;
503  } // end for shwLast
504  } // end for shwNext
505 
506  cand = bak;
507  return false;
508  }
509 
510  //......................................................................
512  const Cand& prev,
513  const View& v,
514  const dt::View::ChunkMap& chunk_map,
515  Direction dir,
516  bool exactHits) const
517  {
518  while(true){
519  if(cur.NSegs() < 2) return false;
520  if(TryExtendCand(cur, prev, v, chunk_map, dir, exactHits)) return true;
521 
522  cur.ShortenOne(dir);
523 
524  // Too much freedom
525  if(cur.IsAllShower()) return false;
526  if(cur.NChunks() == 0) return false; // all dead
527  if(cur.AllHits().empty()) return false;
528  }
529  }
530 
531  //......................................................................
533  const dt::View::ChunkMap& chunk_map,
534  Chain& chain, Direction dir) const
535  {
536  Cand cur = chain.ExtremalCand(dir);
537 
538  while(TryExtendChainOne(cur, chain.ExtremalCand(dir), v, chunk_map, dir, false)){
539  chain.Add(cur, dir);
540  }
541  }
542 
543  //......................................................................
545  const dt::View::ChunkMap& chunk_map,
546  std::vector<rb::Track>* seeds_out,
547  bool& ok) const
548  {
549  const Cand& seed = BestSeed(v, ok);
550  if(!ok) return Chain();
551  seeds_out->push_back(seed.ToTrack());
552 
553  Chain chain(seed);
554  ExtendChain(v, chunk_map, chain, kUpstream);
555  ExtendChain(v, chunk_map, chain, kDownstream);
556  ok = true;
557  return chain;
558  }
559 
560  //......................................................................
561  std::vector<rb::Track>
563  std::vector<rb::Track>* seeds_out,
564  std::vector<rb::Track>* unspliced_out,
565  std::vector<rb::Track>* chains_out) const
566  {
567  std::vector<rb::Track> ret;
568  if(cells.empty()) return ret;
569 
570  std::vector<Chain> chains;
571 
572  dt::View v(cells, fConsiderBadChans);
573 
574  while(true){
575  // TODO, this should be shared
576  const dt::View::ChunkMap chunk_map = v.MakeChunks();
577 
578  bool ok;
579  Chain chain = RecoOneTrack(v, chunk_map, seeds_out, ok);
580 
581  if(!ok) break;
582 
583  chain.TrimEnds();
584  // TODO: how can this happen?
585  if(chain.FirstCand().NSegs() < 2) break;
586  if(chain.LastCand().NSegs() < 2) break;
587 
588  chains.push_back(chain);
589  unspliced_out->push_back(chain.ToTrack());
590 
591  // We're going to try and find another track, but first remove all the
592  // hits we've already explained.
594  for(unsigned int n = 0; n < hits.size(); ++n){
595  v.Channel(hits[n]->Plane(), hits[n]->Cell()).MarkHitFound();
596  }
597 
598  if(fOnlyOneTrack) break;
599  } // end while
600 
601  SpliceChains(chains);
602 
603  unsigned int N = chains.size();
604  for(unsigned int n = 0; n < N; ++n){
605  const std::vector<Chain> remnants = CheckAndFixSparseChain(chains[n]);
606  chains.insert(chains.end(), remnants.begin(), remnants.end());
607  }
608 
609  N = chains.size();
610  for(unsigned int n = 0; n < N; ++n){
611  // Comes out of Truncate() at least sometimes
612  if(chains[n].AllHits().empty()){
613  mf::LogWarning("DiscreteTracker") << "Track with no hits. Skipping.";
614  continue;
615  }
616 
617  if(chains[n].FirstSeg().zL > chains[n].LastSeg().zL){
618  mf::LogWarning("DiscreteTracker") << "Reversed track. Skipping.";
619  continue;
620  }
621 
622  ret.push_back(chains[n].ToTrack());
623  const std::vector<rb::Track> debug = chains[n].ToDebugTracks();
624  chains_out->insert(chains_out->end(), debug.begin(), debug.end());
625  }
626 
627  MergeParasiteTracks(ret);
628 
629  return ret;
630  }
631 
632  //......................................................................
633  bool DiscreteTracker::SpliceIsSane(const Chain& a, const Chain& b) const
634  {
635  // Must not overlap. Can't jump more than half the whole total length
636 
637  if(a.LastPlane() < b.FirstPlane()){
638  if(DotScore(a.LastCand(), b.FirstCand(), true) < fMinDotProduct) return false;
639 
640  return (2*(b.FirstPlane()-a.LastPlane()) < b.LastPlane()-a.FirstPlane());
641  }
642  if(b.LastPlane() < a.FirstPlane()){
643  if(DotScore(b.LastCand(), a.FirstCand(), true) < fMinDotProduct) return false;
644 
645  return (2*(a.FirstPlane()-b.LastPlane()) < a.LastPlane()-b.FirstPlane());
646  }
647 
648  return false;
649  }
650 
651  struct SpliceDesc
652  {
653  SpliceDesc(int cIdx, int pIdx, Chain& ch, int g, int t)
654  : childIdx(cIdx), parentIdx(pIdx), spliced(ch), gap(g), total(t) {}
655 
656  bool operator<(const SpliceDesc& sd) const
657  {
658  // Smaller gaps are better. Break ties with total length of spliced
659  // pieces.
660  if(gap == sd.gap) return total > sd.total;
661  return gap < sd.gap;
662  }
663 
664  int childIdx, parentIdx;
666  int gap;
667  int total;
668  };
669 
670  //......................................................................
671  void DiscreteTracker::SpliceChains(std::vector<Chain>& chains) const
672  {
673  while(true){
674  std::vector<int> parentIdxs(chains.size());
675  for(unsigned int n = 0; n < chains.size(); ++n) parentIdxs[n] = n;
676 
677  std::vector<int> childIdxs(chains.size());
678  for(unsigned int n = 0; n < chains.size(); ++n) childIdxs[n] = n;
679 
680  // Keep the best splice first
681  std::multiset<SpliceDesc> splices;
682 
683  for(unsigned int i = 0; i < chains.size(); ++i){
684  const unsigned int parentIdx = parentIdxs[i];
685  const Chain& parent = chains[parentIdx];
686 
688  for(unsigned int j = 0; j < chains.size(); ++j){
689  const unsigned int childIdx = childIdxs[j];
690 
691  if(childIdx == parentIdx) continue;
692 
693  const Chain& child = chains[childIdx];
694 
695  int gap;
696  if(dir == kDownstream)
697  gap = child.FirstPlane() - parent.LastPlane();
698  else
699  gap = parent.FirstPlane() - child.LastPlane();
700 
701  // Have to get to correct side
702  if(gap < 0) continue;
703 
704  // Bail out early if we know we can't beat the best splice
705  if(!splices.empty() && gap > splices.begin()->gap) continue;
706 
707  if(!SpliceIsSane(parent, child)) continue;
708 
709  View v(child.AllHits(), fConsiderBadChans);
710  // Doing it this way preserves the "shower" flag on chunks
711  const dt::View::ChunkMap chunk_map = child.AsChunkMap();
712 
713  Cand join = parent.ExtremalCand(dir);
714  if(!TryExtendChainOne(join, parent.ExtremalCand(dir),
715  v, chunk_map,
716  dir, true)) continue;
717 
718  // Also need to check that the join is willing to be made in the
719  // other direction. Otherwise you can get a one-hit cand on one
720  // side, that's willing to match almost anything, in crazy ways.
721  View pview(parent.AllHits(), fConsiderBadChans);
722  const dt::View::ChunkMap pmap = parent.AsChunkMap();
723  Cand join2 = child.ExtremalCand(OtherDir(dir));
724  if(!TryExtendChainOne(join2, child.ExtremalCand(OtherDir(dir)),
725  pview, pmap, OtherDir(dir), true)) continue;
726 
727  Chain spliced = parent;
728  spliced.Add(join, dir);
729  // The splice can often subsume the whole child chain
730  if(child.ExtremalPlane(dir) != join.ExtremalPlane(dir)){
731  spliced.Add(child, dir);
732  }
733 
734  splices.insert(SpliceDesc(childIdx, parentIdx, spliced, gap,
735  child.NHitPlanes()+parent.NHitPlanes()));
736  } // end for j (childIdx)
737  } // end for dir
738  } // end for i (parentIdx)
739 
740  if(splices.empty()) return;
741 
742  const SpliceDesc sd = *splices.begin();
743  chains[sd.parentIdx] = sd.spliced;
744  chains.erase(chains.begin()+sd.childIdx);
745  }
746  }
747 
748  //......................................................................
749  Cand DiscreteTracker::BestSeed(const dt::View& v, bool& ok) const
750  {
751  ok = true;
752 
753  // TODO, we could iterate over these instead
754  // const std::set<unsigned int> planes = v.GetPlanes();
755  const dt::View::ChunkMap chunks_map = v.MakeChunks();
756 
757  Cand best(v.GeoView(), true);
758 
759  if(chunks_map.empty()){
760  ok = false;
761  return best;
762  }
763 
764  std::list<Cand> live;
765 
766  // Go through planes
767  for(dt::View::ChunkMap::const_iterator it = chunks_map.begin(); it != chunks_map.end(); ++it){
768  const int plane = it->first;
769 
770  const std::vector<Chunk>& chunks = it->second;
771 
772  std::list<Cand> newLive;
773 
774  for(std::list<Cand>::iterator live_it = live.begin(); live_it != live.end(); ++live_it){
775 
776  if(plane != v.NextPlane(live_it->LastChunk().Plane())){
777  if(live_it->NSegs() > best.NSegs()) best = *live_it;
778  continue;
779  }
780 
781  // TODO: make a less wasteful version of TryAddSeg? Try to keep just
782  // the one list.
783  Cand curCand = *live_it; // Take a copy
784 
785  // Were any chunks any use?
786  bool any = false;
787 
788  for(unsigned int chunkIdx = 0; chunkIdx < chunks.size(); ++chunkIdx){
789  const Chunk& chunk = chunks[chunkIdx];
790  // TODO: potentially real data could have seeds that never had two
791  // adjacent hits?
792  if(chunk.AllHits().empty()) continue;
793 
794  // TODO: have chunks given back in two lists?
795  if(chunk.Up() != curCand.Up()) continue;
796 
797  const std::pair<Segment, Segment>& segs = chunk.GetSegs();
798 
799  if(curCand.TryAddSeg(segs.first, kDownstream)){
800  curCand.AddChunk(chunk, kDownstream);
801  if(curCand.TryAddSeg(segs.second, kDownstream)){
802  // Good, keep going
803  newLive.insert(newLive.end(), curCand);
804  curCand = *live_it; // New one to play with
805  }
806  else{
807  // Failed second half, first might be useful
808  if(curCand.NSegs() > best.NSegs()) best = curCand;
809  curCand = *live_it;
810  }
811  // Got at least the first seg in
812  any = true;
813  }
814  } // end for chunkIdx
815 
816  if(!any){
817  // Got stuck, but it might still be the best
818  if(curCand.NSegs() > best.NSegs()) best = curCand;
819  // No need to make a new curCand anymore
820  }
821  } // end for live_it
822 
823  live.swap(newLive);
824  newLive.clear();
825 
826  // Seed new cands from this chunk
827  for(unsigned int chunkIdx = 0; chunkIdx < chunks.size(); ++chunkIdx){
828  const Chunk& chunk = chunks[chunkIdx];
829 
830  const std::pair<Segment, Segment>& segs = chunk.GetSegs();
831 
832  Cand newcand(v.GeoView(), chunk.Up());
833  newcand.AddChunk(chunk, kDownstream);
834  bool startok = newcand.TryAddSeg(segs.first, kDownstream);
835  assert(startok);
836  startok = startok && newcand.TryAddSeg(segs.second, kDownstream);
837  assert(startok);
838  live.insert(live.end(), newcand);
839 
840  // Only the second seg
841  Cand newcand2(v.GeoView(), chunk.Up());
842  newcand2.AddChunk(chunk, kDownstream);
843  startok = startok && newcand2.TryAddSeg(segs.second, kDownstream);
844  live.insert(live.end(), newcand2);
845 
846  assert(startok);
847  } // end for chunkIdx
848  } // end for it
849 
850  // Maybe one that made it all the way to the end is the best?
851  for(std::list<Cand>::iterator live_it = live.begin(); live_it != live.end(); ++live_it){
852  if(live_it->NSegs() > best.NSegs()) best = *live_it;
853  }
854 
855  if(best.NSegs() == 0) ok = false;
856  // TODO: how can this even happen? But if it does we'll get stuck in an
857  // endless loop, so bail out of the rest of the event.
858  if(best.AllHits().empty()) ok = false;
859  return best;
860  }
861 
862  //......................................................................
864  {
865  std::vector<Chain> remnants;
866 
867  const std::vector<Chunk> chunks = chain.AllChunks();
868  assert(!chunks.empty());
869 
870  std::set<int> planes;
871  for(unsigned int n = 0; n < chunks.size(); ++n)
872  planes.insert(chunks[n].Plane());
873 
874  bool progress;
875  do{
876  progress = false;
878  const Direction dir = OtherDir(end);
879 
880  int nTot = 0;
881  int nMiss = 0;
882 
883  int curPlane = chain.ExtremalPlane(end);
884 
885  if(chain.AllHits().empty()) break;
886  View v(chain.AllHits(), fConsiderBadChans);
887 
888  do{
889  ++nTot;
890 
891  const bool hitHere = planes.count(curPlane);
892  if(!hitHere) ++nMiss;
893 
894  // If we've accumulated more than 1/2 gaps, and have now come back
895  // into hits
896  if(hitHere && nMiss > 0){
897  // Want to allow HGH... and ban GH... and HGH
898  if(nTot == nMiss || (nMiss > 1 && nMiss*2 >= nTot-1)){
899  Chain remnant = chain;
900  remnant.Truncate(curPlane+((dir == kUpstream) ? +1 : -1), dir);
901  remnants.push_back(remnant);
902 
903  chain.Truncate(curPlane, end);
904  progress = true;
905  // Shouldn't look any further from this end, but try the other one
906  break;
907  }
908  }
909 
910  curPlane = v.AdjacentPlane(curPlane, dir);
911  } while(curPlane != -1);
912  } // end for end/dir
913  } while(progress);
914 
915  return remnants;
916  }
917 
918  bool CompareByLength(const rb::Track& a, const rb::Track& b)
919  {
920  return a.TotalLength() < b.TotalLength();
921  }
922 
923  //......................................................................
924  rb::Vertex* DiscreteTracker::FindVertex(std::vector<rb::Track> xTrks,
925  std::vector<rb::Track> yTrks) const
926  {
927  if(xTrks.empty() || yTrks.empty()) return 0;
928 
929  // Longest tracks first. Only used for fallback at the end
930  std::sort(xTrks.rbegin(), xTrks.rend(), CompareByLength);
931  std::sort(yTrks.rbegin(), yTrks.rend(), CompareByLength);
932 
933  double bestVx, bestVy, bestVz1, bestVz2;
934  // I don't believe there's a path that doesn't set these, but need this to
935  // shut the compiler up in optimized builds.
936  bestVx = bestVy = bestVz1 = bestVz2 = 0;
937 
938  for(geo::View_t view: {geo::kX, geo::kY}){
939  double bestScore = 1e10;
940  const std::vector<rb::Track>& trks = (view == geo::kX) ? xTrks : yTrks;
941  for(unsigned int idx1 = 0; idx1 < trks.size()-1; ++idx1){
942  if(trks[idx1].ExtentPlane() < 3) continue;
943  for(bool end1: {false, true}){
944  for(unsigned int idx2 = idx1+1; idx2 < trks.size(); ++idx2){
945  if(trks[idx2].ExtentPlane() < 3) continue;
946  for(bool end2: {false, true}){
947 
948  // Start and direction defining point of these two lines
949  const TVector3 s1 = end1 ? trks[idx1].Stop() : trks[idx1].Start();
950  const TVector3 s2 = end2 ? trks[idx2].Stop() : trks[idx2].Start();
951  const TVector3 d1 = trks[idx1].TrajectoryPoint(end1 ? trks[idx1].NTrajectoryPoints()-2 : 1);
952  const TVector3 d2 = trks[idx2].TrajectoryPoint(end2 ? trks[idx2].NTrajectoryPoints()-2 : 1);
953 
954  double vxy, vz;
955  const bool ok = geo::LineIntersection(s1[view], s1.Z(), d1[view], d1.Z(),
956  s2[view], s2.Z(), d2[view], d2.Z(),
957  vxy, vz);
958  if(!ok) continue;
959 
960  const double dist1 = util::pythag(s1[view]-vxy, s1.Z()-vz);
961  const double dist2 = util::pythag(s2[view]-vxy, s2.Z()-vz);
962 
963  // Favour both long tracks and hits close to the ends
964  const double score = dist1/trks[idx1].TotalLength()+dist2/trks[idx2].TotalLength();
965  if(score < bestScore){
966  bestScore = score;
967  if(view == geo::kX){
968  bestVx = vxy;
969  bestVz1 = vz;
970  }
971  else{
972  bestVy = vxy;
973  bestVz2 = vz;
974  }
975  } // end if better
976  } // end for end2
977  } // end for idx2
978  } // end for end1
979  } // end for idx1
980  // Too far from the ends, fallback
981  if(bestScore > 6){
982  if(view == geo::kX){
983  bestVx = trks[0].Start().X();
984  bestVz1 = trks[0].Start().Z();
985  }
986  else{
987  bestVy = trks[0].Start().Y();
988  bestVz2 = trks[0].Start().Z();
989  }
990  }
991  } // end for view
992 
993  // Time is an afterthought
994  return new rb::Vertex(bestVx, bestVy, (bestVz1+bestVz2)/2, xTrks[0].MinTNS());
995  }
996 
997  //......................................................................
998  void DiscreteTracker::MergeParasiteTracks(std::vector<rb::Track>& tracks) const
999  {
1000  // Worried this could be too aggressive around the vertex
1001 
1002  const unsigned int N = tracks.size();
1003  if(N < 2) return;
1004 
1005  // Longest first
1006  std::sort(tracks.begin(), tracks.end(),
1007  [](const rb::Track& a, const rb::Track& b){
1008  return a.ExtentPlane() > b.ExtentPlane();
1009  });
1010 
1011  for(unsigned int hostIdx = 0; hostIdx < N-1; ++hostIdx){
1012  const rb::HitMap hmap(&tracks[hostIdx]);
1013  for(unsigned int paraIdx = hostIdx+1; paraIdx < N; ++paraIdx){
1014  const rb::Track& para = tracks[paraIdx];
1015  bool match = true;
1016  for(unsigned int i = 0; i < para.NCell(); ++i){
1017  const art::Ptr<rb::CellHit> chit = para.Cell(i);
1018  if(!hmap.CellExists(chit->Plane(), chit->Cell()-1) &&
1019  !hmap.CellExists(chit->Plane(), chit->Cell()+1)){
1020  match = false;
1021  break;
1022  }
1023  } // end for i
1024  if(match){
1025  // Make sure there are no duplicates in the final host track
1026  rb::HitMap hmap(tracks[hostIdx].AllCells());
1027  for(unsigned int n = 0; n < para.NCell(); ++n){
1028  const art::Ptr<rb::CellHit> chit = para.Cell(n);
1029  if(!hmap.CellExists(chit->Plane(), chit->Cell()))
1030  tracks[hostIdx].Add(chit);
1031  }
1032  tracks.erase(tracks.begin()+paraIdx);
1033  MergeParasiteTracks(tracks);
1034  return;
1035  }
1036  } // end for paraIdx
1037  } // end for hostIdx
1038  }
1039 
1041 
1042 } //namespace dt
1043 ////////////////////////////////////////////////////////////////////////
T max(const caf::Proxy< T > &a, T b)
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
double DotScore(const Cand &c1, const Cand &c2, bool allowFlip=false) const
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.
void MergeParasiteTracks(std::vector< rb::Track > &tracks) const
bool LineIntersection(double x0, double y0, double x1, double y1, double X0, double Y0, double X1, double Y1, double &x, double &y)
Find the intersection between two line segments.
Definition: Geo.cxx:184
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
TH1 * fAngle
Opening angle of all cands considered for joining.
int NHitPlanes() const
Definition: Chain.cxx:266
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
bool Up() const
Definition: Cand.h:61
FILE * fDebugFile
File that track details are written to.
std::pair< Segment, Segment > GetSegs() const
Definition: Chunk.cxx:102
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
Sequence of contiguous hits and dead cells all on the same plane.
Definition: Chunk.h:17
bool IsShowerChunk() const
Definition: Chunk.h:43
double MinGradient() const
Definition: Cand.h:36
Representation of the channels in one view.
Definition: View.h:22
const char * p
Definition: xmltok.h:285
bool operator<(const SpliceDesc &sd) const
T sqrt(T number)
Definition: d0nt_math.hpp:156
bool CompareByLength(const rb::Track &a, const rb::Track &b)
unsigned int FirstPlane() const
Definition: Chain.cxx:167
Vertical planes which measure X.
Definition: PlaneGeo.h:28
Definition: Chain.py:1
bool fConsiderBadChans
Passed to View.
unsigned int Ncells() const
Number of cells in this plane.
Definition: PlaneGeo.h:43
unsigned int ExtremalPlane(Direction dir) const
Definition: Cand.h:68
dt::Channel Channel(int plane, int cell) const
Definition: View.cxx:128
unsigned int FirstPlane() const
Definition: Cand.h:66
geo::View_t GeoView() const
Definition: View.h:26
bool IsAllShower() const
Definition: Cand.cxx:415
const Chunk & ExtremalChunk(Direction dir) const
Definition: Cand.h:46
Chain RecoOneTrack(const View &v, const dt::View::ChunkMap &chunk_map, std::vector< rb::Track > *seeds_out, bool &ok) const
std::map< int, std::vector< Chunk > > ChunkMap
Definition: View.h:27
bool progress
Insert implicit nodes #####.
bool left
Definition: Segment.h:41
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
const Chunk & FirstChunk() const
Definition: Cand.cxx:353
const Segment & LastSeg() const
Definition: Cand.cxx:371
const PlaneGeo * Plane(unsigned int i) const
bool Up() const
Definition: Chunk.cxx:49
DEFINE_ART_MODULE(TestTMapFile)
double DotScoreHelper(double m1, double m2) const
int AdjacentPlane(unsigned int plane, Direction dir) const
Definition: View.cxx:121
T acos(T number)
Definition: d0nt_math.hpp:54
void ExtendChain(const View &v, const dt::View::ChunkMap &chunk_map, Chain &chain, Direction dir) const
double fMinDotProduct
Joined Cands that only have one chunk in common must have a dot-product better than this...
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
unsigned int NChunks() const
Definition: Cand.h:41
SpliceDesc(int cIdx, int pIdx, Chain &ch, int g, int t)
Provides efficient lookup of CellHits by plane and cell number.
Definition: HitMap.h:22
unsigned int LastPlane() const
Definition: Chain.cxx:173
int NextPlane(unsigned int plane) const
Definition: View.cxx:111
DiscreteTracker(const fhicl::ParameterSet &pset)
c2
Definition: demo5.py:33
ProductID put(std::unique_ptr< PROD > &&product)
Definition: Event.h:102
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
fclose(fg1)
double sd(Eigen::VectorXd x)
unsigned short Cell() const
Definition: CellHit.h:40
std::vector< Chain > CheckAndFixSparseChain(Chain &chain) const
bool SpliceIsSane(const Chain &a, const Chain &b) const
std::vector< Segment > PossibleNextSegs(const View &v, const dt::View::ChunkMap &chunk_map, const Cand &cand, Direction dir, std::vector< Chunk > &chunks_out, bool shw, bool exactHits) const
Definition: Cand.cxx:23
std::vector< Chunk > AllChunks() const
Definition: Chain.cxx:245
ChunkMap MakeChunks() const
Definition: View.cxx:148
rb::Track ToTrack() const
Definition: Chain.cxx:29
unsigned int LastPlane() const
Definition: Cand.h:67
void hits()
Definition: readHits.C:15
art::PtrVector< rb::CellHit > AllHits() const
Definition: Chunk.h:37
Calculation and representation of a straight line passing through several "segment" windows...
Definition: Cand.h:25
void produce(art::Event &evt)
void prev()
Definition: show_event.C:91
bool TryExtendCandOne(Cand &cand, const View &v, const dt::View::ChunkMap &chunk_map, Direction dir, const Cand &prev_dir, bool shw, bool exactHits) const
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
unsigned int seed
Definition: runWimpSim.h:102
bool HighActivity(const art::Ptr< rb::Cluster > &clust) const
Cut to prevent choking on large air shower events.
Job module performing "discrete" track reconstruction.
const double a
void SetShowerChunk()
Definition: Chunk.h:44
void SpliceChains(std::vector< Chain > &chains) const
T get(std::string const &key) const
Definition: ParameterSet.h:231
Geometry information for a single readout plane.
Definition: PlaneGeo.h:36
rb::Track ToTrack() const
Definition: Cand.cxx:453
chain
Check that an output directory exists.
rb::Vertex * FindVertex(std::vector< rb::Track > xTrks, std::vector< rb::Track > yTrks) const
Toy vertex finder.
bool fObeyPreselection
Check rb::IsFiltered?
Collect Geo headers and supply basic geometry functions.
art::PtrVector< rb::CellHit > AllHits() const
Definition: Cand.cxx:424
bool IsWorthShowering() const
Definition: Chunk.cxx:194
bool IsFiltered(const art::Event &evt, art::Ptr< T > x, const std::vector< std::string > &labels)
Is this Ptr marked "filtered out"?
Definition: FilterList.h:96
const Segment & FirstSeg() const
Definition: Cand.cxx:365
const double j
Definition: BetheBloch.cxx:29
EventNumber_t event() const
Definition: Event.h:67
Direction OtherDir(Direction d)
Definition: Types.h:7
int NSegs() const
Definition: Cand.h:85
Cand MaybeFlip(bool &ok) const
Definition: Cand.cxx:547
std::vector< rb::Track > RecoDiscrete(const art::PtrVector< rb::CellHit > &cells, std::vector< rb::Track > *seeds_out, std::vector< rb::Track > *unspliced_out, std::vector< rb::Track > *chains_out) const
Vertex location in position and time.
Direction
Definition: Types.h:5
void Add(const Cand &c)
Definition: Chain.cxx:142
Definition: View.py:1
dt::View::ChunkMap AsChunkMap() const
Definition: Chain.cxx:337
bool empty() const
Definition: PtrVector.h:336
size_type size() const
Definition: PtrVector.h:308
bool IsExtremalChunkComplete(Direction dir) const
Definition: Cand.cxx:377
void MarkExtremalChunkShower(Direction dir)
Definition: Cand.cxx:384
int fVerticalTrackMinCells
The number of adjacent cells in a single plane, otherwise not fit into any track, required to declare...
const Chunk & LastChunk() const
Definition: Cand.cxx:359
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
void reconfigure(const fhicl::ParameterSet &pset)
TH1 * fSinglesAngle
For all cands with only one overlapping chunk.
static constexpr Double_t m2
Definition: Munits.h:145
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void MarkHitFound()
Definition: Channel.cxx:32
bool TryAddSeg(const Segment &seg, Direction dir, bool check=true)
Appends seg to fSegs and updates internal state reflecting limits on the possible straight lines...
Definition: Cand.cxx:150
art::PtrVector< rb::CellHit > AllHits() const
Definition: Chain.cxx:220
T * make(ARGS...args) const
void SetUp(bool up)
Definition: Chunk.cxx:55
unsigned int ExtremalPlane(Direction dir) const
Definition: Chain.cxx:179
void TrimEnds()
Remove Cands that add no actual hits. After building but before writing.
Definition: Chain.cxx:277
TDirectory * dir
Definition: macro.C:5
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
void AddChunk(const Chunk &chunk, Direction dir)
Simply store the chunk in the list. Need to call AddSeg too.
Definition: Cand.cxx:36
Cand BestSeed(const dt::View &v, bool &ok) const
bool TryExtendCand(Cand &cand, const Cand &old_cand, const View &v, const dt::View::ChunkMap &chunk_map, Direction dir, bool exactHits) const
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
void geom(int which=0)
Definition: geom.C:163
const hit & b
Definition: hits.cxx:21
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
assert(nhit_max >=nhit_nbins)
c1
Definition: demo5.py:24
const Cand & FirstCand() const
Definition: Chain.h:45
bool fOnlyOneTrack
Stop after the first track in each view.
bool TryExtendChainOne(Cand &cur, const Cand &prev, const View &v, const dt::View::ChunkMap &chunk_map, Direction dir, bool exactHits) const
bool CellExists(unsigned int planeIdx, unsigned int cellIdx) const
Does the map contain any cell at this position?
Definition: HitMap.cxx:81
void Truncate(unsigned int plane, Direction end)
Definition: Chain.cxx:306
Window the line must pass through from (z,y0)-(z,y1)
Definition: Segment.h:24
unsigned int NPlanes() const
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
void ShortenOne(Direction dir)
Definition: Cand.cxx:533
bool IsBad(int plane, int cell)
const Cand & ExtremalCand(Direction dir) const
Definition: Chain.h:47
const Cand & LastCand() const
Definition: Chain.h:46
double MaxGradient() const
Definition: Cand.h:37
std::string fClusterInput
Module name to look for Clusters under.