ViewMerger_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file ViewMerger.cxx
3 // \author Christopher Backhouse - bckhouse@caltech.edu
4 ////////////////////////////////////////////////////////////////////////
5 
6 #include <algorithm>
7 
8 // NOvASoft includes
10 #include "Geometry/Geometry.h"
11 #include "GeometryObjects/Geo.h"
12 
13 #include "RecoBase/Track.h"
14 #include "RecoBase/Vertex.h"
15 #include "RecoBase/Cluster.h"
16 #include "RecoBase/HitMap.h"
17 #include "RecoBase/CellHit.h"
18 #include "RecoBase/RecoHit.h"
19 #include "RecoBase/FilterList.h"
20 
21 #include "Utilities/AssociationUtil.h"
22 #include "Utilities/func/MathUtil.h"
23 
24 // Framework includes
28 #include "fhiclcpp/ParameterSet.h"
34 
35 namespace dt
36 {
37  /// Match 2D tracks into 3D tracks
38  class ViewMerger : public art::EDProducer
39  {
40  public:
41  explicit ViewMerger(const fhicl::ParameterSet& pset);
42  ~ViewMerger();
43 
44  virtual void reconfigure(const fhicl::ParameterSet& pset);
45 
46  virtual void beginJob();
47  virtual void endJob();
48 
49  virtual void produce(art::Event& evt);
50 
51  protected:
52  /// \brief Do most of the work. Called from \ref produce
53  ///
54  /// \param sliceMap All the hits in this slice
55  /// \param xtracks The 2D tracks from the x-view to merge. Those used will
56  /// be removed from the list. Any left unmerged will remain in the list.
57  /// \param ytracks Ditto
58  /// \return The 3D tracks produced
59  std::vector<rb::Track>
60  MatchTracks(const rb::HitMap& sliceMap,
61  std::list<rb::Track>& xtracks,
62  std::list<rb::Track>& ytracks) const;
63 
64  /// Find vertical tracks seperated by one plane. Return results of merge
65  std::vector<rb::Track>
66  MatchVerticals(std::list<rb::Track>& xtracks,
67  std::list<rb::Track>& ytracks) const;
68 
69  /// Tilt tracks arbitrarily and merge them into a 3D track
71  const rb::Track& ytrk) const;
72 
73  /// Corrected for direction cosines, so if you divide this by extent in Z
74  /// you'd get dE/dx
75  void TotalChargePerView(const rb::Track& zip,
76  double& totx, double& toty) const;
77 
78  /// If we extend tracks to match in length, how many cells do they pass
79  /// over but don't hit? Use the mean energy deposit in these tracks to
80  /// estimate how much charge has gone missing, in GeV.
81  double MissingChargeByExtension(const rb::HitMap& sliceMap,
82  const rb::Track& aTrk,
83  const rb::Track& bTrk) const;
84 
85  double ScoreForExtension(const rb::HitMap& sliceMap,
86  const rb::Track& xTrk,
87  const rb::Track& yTrk) const;
88 
89  double ScoreForJoinPlusExtension(const rb::HitMap& sliceMap,
90  const rb::Track& single,
91  const rb::Track& parta,
92  const rb::Track& partb) const;
93 
94  /// Append two tracks in the same view (crudely).
95  rb::Track JoinTracks(const rb::Track& parta, const rb::Track& partb) const;
96 
97  /// \a hmap should contain cells that are hit and so OK to pass through
98  int CountUnexplainedOnLine(const rb::HitMap& hmap,
99  TVector3 start, TVector3 end,
100  geo::View_t view) const;
101 
102  double ChargeDifferenceBetweenViews(const rb::Track& zip) const;
103 
104  /// Extend each track to the vertex point, only if it's completely
105  /// uncontroversial.
106  void MaybeExtendToVertex(const rb::Vertex* vtx,
107  std::vector<rb::Track>& tracks) const;
108 
109  /// Find number of overlapping planes, return unexp/overlap
110  double MatchScore(double unexp,
111  const rb::Track& aTrk,
112  const rb::Track& bTrk) const;
113 
114  /// Make a 3D track based on two 2D tracks
115  rb::Track ZipTracks(const rb::Track& aTrk, const rb::Track& bTrk) const;
116 
117  /// Helper for \ref Adjacent. Fills \a mins and \a maxs with minimum and
118  /// maximum cell number in each plane of \a trk
119  void Extremes(const rb::Track& trk,
120  std::map<int, int>& mins, std::map<int, int>& maxs) const;
121 
122  /// Helper for \ref MergeAdjacentTracks. Are all hits on one track adjacent
123  /// to the other one?
124  bool Adjacent(const rb::Track& a, const rb::Track& b) const;
125 
126  /// Any tracks that are adjacent along the whole length of the shorter one
127  /// will be merged
128  void MergeAdjacentTracks(std::list<rb::Track>& trks) const;
129 
130  /// \brief Toy vertex finder
131  ///
132  /// Uses the intersection of the 3D tracks whose intersection is
133  /// proportionally nearest their ends.
134  rb::Vertex* FindVertex(std::vector<rb::Track> trks) const;
135 
136  /// \brief Turn a 2D track into a 3D track
137  ///
138  /// Uses \a merged to guess a W position
140  const std::vector<rb::Track>& merged) const;
141 
142  std::string fSliceInput; ///< Module name to look for slices under
143  std::string fTrackInput; ///< Module name to look for 2D tracks under
144 
145  bool fObeyPreselection; ///< Check rb::IsFiltered?
146 
147  FILE* fDebugFile;
148  };
149 
150 
152  : fDebugFile(0)
153  {
154  produces<std::vector<rb::Track> >();
155  produces<std::vector<rb::Vertex> >();
156  // Tracks are associated with the slice they came from
157  produces<art::Assns<rb::Track, rb::Cluster> >();
158 
159  reconfigure(pset);
160  }
161 
162  //......................................................................
164  {
165  }
166 
167  //......................................................................
169  {
170  fSliceInput = pset.get<std::string>("SliceInput");
171  fTrackInput = pset.get<std::string>("TrackInput");
172 
174  fDebugFile = fopen(pset.get<std::string>("DebugFile").c_str(), "w");
175 
176  fObeyPreselection = pset.get<bool>("ObeyPreselection");
177  }
178 
179  //......................................................................
181  {
182  }
183 
184  //......................................................................
186  {
188  }
189 
190  //......................................................................
192  {
193  std::unique_ptr<std::vector<rb::Track> > trackcol(new std::vector<rb::Track>);
194  std::unique_ptr<std::vector<rb::Vertex> > vtxcol(new std::vector<rb::Vertex>);
195  // Associate tracks with the slice they came from
196  std::unique_ptr<art::Assns<rb::Track, rb::Cluster> > assns(new art::Assns<rb::Track, rb::Cluster>);
197 
199  evt.getByLabel(fSliceInput, hSlices);
200 
202  for(unsigned int i = 0; i < hSlices->size(); ++i)
203  slices.push_back(art::Ptr<rb::Cluster>(hSlices, i));
204 
205  const unsigned int sliceMax = slices.size();
207  for(unsigned int sliceIdx = 0; sliceIdx < sliceMax; ++sliceIdx){
208  if(slices[sliceIdx]->IsNoise()) continue;
209 
210  if(fObeyPreselection && rb::IsFiltered(evt, slices[sliceIdx])) continue;
211 
212  // Get all the tracks associated with this slice
213  std::vector< art::Ptr<rb::Track> > pTracks = fmt.at(sliceIdx);
214 
215  if(pTracks.empty()) continue;
216 
217  const rb::HitMap sliceMap(slices[sliceIdx]->AllCells());
218 
219  // In a nicer form
220  std::list<rb::Track> xtracks, ytracks;
221  for(unsigned int trkIdx = 0; trkIdx < pTracks.size(); ++trkIdx){
222  // No sensible way to merge these, and they make a mess
223  if(pTracks[trkIdx]->NCell() == 1) continue;
224  if(pTracks[trkIdx]->View() == geo::kX)
225  xtracks.push_back(*pTracks[trkIdx]);
226  else
227  ytracks.push_back(*pTracks[trkIdx]);
228  } // end for trkIdx
229 
230  MergeAdjacentTracks(xtracks);
231  MergeAdjacentTracks(ytracks);
232 
233  // Actually do all the work
234  std::vector<rb::Track> matchedTracks = MatchTracks(sliceMap,
235  xtracks, ytracks);
236  std::vector<rb::Track> verts = MatchVerticals(xtracks, ytracks);
237  matchedTracks.insert(matchedTracks.end(), verts.begin(), verts.end());
238 
239  // Write summary to debug file
240  fprintf(fDebugFile, "%d: ", evt.event());
241 
242  if(!matchedTracks.empty()){
243  fprintf(fDebugFile, "3d: ");
244  for(unsigned int n = 0; n < matchedTracks.size(); ++n){
245  fprintf(fDebugFile, "%u ", matchedTracks[n].NCell());
246  }
247  }
248 
249  if(!xtracks.empty() || !ytracks.empty()){
250  fprintf(fDebugFile, "2d: ");
251  typedef std::list<rb::Track>::iterator it_t;
252  for(it_t it = xtracks.begin(); it != xtracks.end(); ++it)
253  fprintf(fDebugFile, "%u(x) ", it->NCell());
254  for(it_t it = ytracks.begin(); it != ytracks.end(); ++it)
255  fprintf(fDebugFile, "%u(y) ", it->NCell());
256  }
257 
258  fprintf(fDebugFile, "\n");
259 
260  const rb::Vertex* vtx = FindVertex(matchedTracks);
261  if(vtx){
262  MaybeExtendToVertex(vtx, matchedTracks);
263  vtxcol->push_back(*vtx);
264  }
265  delete vtx;
266 
267  // This can look pretty silly and I'm not sure it adds anything.
268  /*
269  // Invent a third dimension for the unmatched 2D tracks and include them
270  // in the list.
271  std::vector<rb::Track> fakes;
272 
273  for(std::list<rb::Track>::iterator it = xtracks.begin(); it != xtracks.end(); ++it){
274  fakes.push_back(FakeThirdDimension(*it, matchedTracks));
275  }
276  for(std::list<rb::Track>::iterator it = ytracks.begin(); it != ytracks.end(); ++it){
277  fakes.push_back(FakeThirdDimension(*it, matchedTracks));
278  }
279 
280  matchedTracks.insert(matchedTracks.end(), fakes.begin(), fakes.end());
281  */
282 
283  for(unsigned int n = 0; n < matchedTracks.size(); ++n){
284  trackcol->push_back(matchedTracks[n]);
285  util::CreateAssn(*this, evt, *trackcol, slices[sliceIdx], *assns);
286  }
287  } // end for sliceIdx
288 
289  evt.put(std::move(trackcol));
290  evt.put(std::move(vtxcol));
291  evt.put(std::move(assns));
292  }
293 
294  //......................................................................
295  std::vector<rb::Track> ViewMerger::
296  MatchTracks(const rb::HitMap& sliceMap,
297  std::list<rb::Track>& xtracks,
298  std::list<rb::Track>& ytracks) const
299  {
300  // This is what we return
301  std::vector<rb::Track> matchedTracks;
302 
303  typedef std::list<rb::Track>::iterator it_t;
304 
305  // Keep going until we stop being able to match anything
306  while(true){
307  double bestScore = 1e10;
308  it_t bestX = xtracks.end();
309  it_t bestY = ytracks.end();
310  it_t bestJoin = xtracks.end(); // Picking this end() is just convention
311 
312  for(it_t xTrk = xtracks.begin(); xTrk != xtracks.end(); ++xTrk){
313  assert(xTrk->View() == geo::kX);
314  // No way to sensibly extend the trajectory for these right now
315  if(xTrk->ExtentPlane() == 1) continue;
316 
317  for(it_t yTrk = ytracks.begin(); yTrk != ytracks.end(); ++yTrk){
318  assert(yTrk->View() == geo::kY);
319  if(yTrk->ExtentPlane() == 1) continue;
320 
321  // Simple merge, join in x, join in y
322  for(int mode = 0; mode < 3; ++mode){
323  // Don't try to join tracks within a view. DiscreteTracker got
324  // better at this, and right now we're just making a mess when we
325  // try.
326  if(mode == 1 || mode == 2) continue;
327 
328  it_t joinTrk;
329  if(mode == 0) joinTrk = xtracks.end();
330  if(mode == 1){joinTrk = xTrk; ++joinTrk;}
331  if(mode == 2){joinTrk = yTrk; ++joinTrk;}
332 
333  // For simple mode, go once with joinTrk invalid. For other modes,
334  // loop through all joinable tracks in the right view.
335  for(; mode == 0 ||
336  (mode == 1 && joinTrk != xtracks.end()) ||
337  (mode == 2 && joinTrk != ytracks.end()); ++joinTrk){
338  // This could be done smarter if we interfaced with
339  // DiscreteTracker properly. It has a lot more information about
340  // these tracks that we could use.
341  double score;
342  if(mode == 0) score = ScoreForExtension(sliceMap, *xTrk, *yTrk);
343  if(mode == 1) score = ScoreForJoinPlusExtension(sliceMap, *yTrk,
344  *xTrk, *joinTrk);
345  if(mode == 2) score = ScoreForJoinPlusExtension(sliceMap, *xTrk,
346  *yTrk, *joinTrk);
347 
348  if(score < bestScore){
349  bestScore = score;
350  bestX = xTrk;
351  bestY = yTrk;
352  bestJoin = joinTrk;
353  }
354 
355  if(mode == 0) break;
356  } // end for joinTrk
357  } // end for mode
358  } // end for yTrk
359  } // end for xTrk
360 
361  // No matches left
362  if(bestX == xtracks.end()) break;
363 
364  LOG_DEBUG("ViewMerger") << "Best score: " << bestScore << std::endl;
365 
366  if(bestJoin != xtracks.end()){
367  if(bestJoin->View() == geo::kX)
368  matchedTracks.push_back(ZipTracks(JoinTracks(*bestX, *bestJoin),
369  *bestY));
370  else
371  matchedTracks.push_back(ZipTracks(*bestX,
372  JoinTracks(*bestY, *bestJoin)));
373  }
374  else{
375  matchedTracks.push_back(ZipTracks(*bestX, *bestY));
376  }
377 
378  // Can't use either 2D track again
379  xtracks.erase(bestX);
380  ytracks.erase(bestY);
381  if(bestJoin != xtracks.end()){
382  if(bestJoin->View() == geo::kX)
383  xtracks.erase(bestJoin);
384  else
385  ytracks.erase(bestJoin);
386  }
387  } // end while
388 
389  return matchedTracks;
390  }
391 
392  //......................................................................
393  std::vector<rb::Track> ViewMerger::
394  MatchVerticals(std::list<rb::Track>& xtracks,
395  std::list<rb::Track>& ytracks) const
396  {
397  std::vector<rb::Track> ret;
398 
399  typedef std::list<rb::Track>::iterator it_t;
400 
401  bool progress = true;
402  while(progress){
403  progress = false;
404 
405  for(it_t xTrk = xtracks.begin(); xTrk != xtracks.end(); ++xTrk){
406  assert(xTrk->View() == geo::kX);
407  if(xTrk->ExtentPlane() != 1) continue;
408 
409  for(it_t yTrk = ytracks.begin(); yTrk != ytracks.end(); ++yTrk){
410  assert(yTrk->View() == geo::kY);
411  if(yTrk->ExtentPlane() != 1) continue;
412 
413  if(std::abs(1.*xTrk->MinPlane() - 1.*yTrk->MinPlane()) != 1.) continue;
414 
415  ret.push_back(ZipVerticalTracks(*xTrk, *yTrk));
416  xtracks.erase(xTrk);
417  ytracks.erase(yTrk);
418  progress = true;
419  break;
420  } // end for yTrk
421  if(progress) break;
422  } // end for xTrk
423  } // end while
424 
425  return ret;
426  }
427 
428  //......................................................................
430  const rb::Track& ytrk) const
431  {
433 
434  double xyz[3], dxyz[3];
435  geom->CellInfo(xtrk.MinPlane(), xtrk.MinCell(geo::kX), 0, xyz, dxyz);
436  const double x0 = xyz[0];
437  const double xz0 = xyz[2]-dxyz[2];
438  geom->CellInfo(xtrk.MinPlane(), xtrk.MaxCell(geo::kX), 0, xyz, dxyz);
439  const double x1 = xyz[0];
440  const double xz1 = xyz[2]+dxyz[2];
441 
442  geom->CellInfo(ytrk.MinPlane(), ytrk.MinCell(geo::kY), 0, xyz, dxyz);
443  const double y0 = xyz[1];
444  const double yz0 = xyz[2]-dxyz[2];
445  geom->CellInfo(ytrk.MinPlane(), ytrk.MaxCell(geo::kY), 0, xyz, dxyz);
446  const double y1 = xyz[1];
447  const double yz1 = xyz[2]+dxyz[2];
448 
449  // If we had any idea where the vertex was we could potentially swap x0, x1
450  // or y0, y1.
451 
452  const double dxdz = (x1-x0)/(xz1-xz0);
453  const double dydz = (y1-y0)/(yz1-yz0);
454 
455  TVector3 start, stop;
456 
457  if(xz0+xz1 < yz0+yz1){
458  // x is upstream of y
459  start.SetX(x0);
460  start.SetY(y0+dydz*(xz0-yz0));
461  start.SetZ(xz0);
462 
463  stop.SetX(x0+dxdz*(yz1-xz0));
464  stop.SetY(y1);
465  stop.SetZ(yz1);
466  }
467  else{
468  // y is upstream of x
469  start.SetX(x0+dxdz*(yz0-xz0));
470  start.SetY(y0);
471  start.SetZ(yz0);
472 
473  stop.SetX(x1);
474  stop.SetY(y1+dydz*(xz1-yz1));
475  stop.SetZ(xz1);
476  }
477 
478  rb::Track ret(xtrk.AllCells(), start, stop-start);
479  ret.Add(ytrk.AllCells());
480  ret.AppendTrajectoryPoint(stop);
481 
482  return ret;
483  }
484 
485  //......................................................................
487  double& totx, double& toty) const
488  {
489  // TODO: accounting for uncalibrated hits and dead channels
490 
491  totx = 0;
492  toty = 0;
493 
494  for(unsigned int n = 0; n < zip.NCell(); ++n){
495  const art::Ptr<rb::CellHit> chit = zip.Cell(n);
496  const rb::RecoHit rhit = zip.RecoHit(chit);
497  if(rhit.IsCalibrated()){
498  // Estimate track direction here
499  const double eps = 1e-3;
500  double x0, y0, x1, y1;
501  zip.InterpolateXY(rhit.Z(), x0, y0);
502  zip.InterpolateXY(rhit.Z()+eps, x1, y1);
503  const double dz = eps/util::pythag(eps, x1-x0, y1-y0);
504 
505  if(chit->View() == geo::kX)
506  totx += rhit.GeV()*dz;
507  else
508  toty += rhit.GeV()*dz;
509  }
510  }
511  }
512 
513  //......................................................................
515  const rb::Track& xTrk,
516  const rb::Track& yTrk) const
517  {
518  // UnexplainedByJoinPlusExtension doesn't know any better. Easier to fix it
519  // here.
520  if(xTrk.View() == geo::kY){
521  assert(yTrk.View() == geo::kX);
522  return MissingChargeByExtension(slice, yTrk, xTrk);
523  }
524 
525  const rb::Track& zip = ZipTracks(xTrk, yTrk);
526 
527  const int unexpX =
528  CountUnexplainedOnLine(slice, zip.Start(), xTrk.Start(), geo::kX)+
529  CountUnexplainedOnLine(slice, zip.Stop(), xTrk.Stop(), geo::kX);
530 
531  const int unexpY =
532  CountUnexplainedOnLine(slice, zip.Start(), yTrk.Start(), geo::kY)+
533  CountUnexplainedOnLine(slice, zip.Stop(), yTrk.Stop(), geo::kY);
534 
535  // End up passing through too many live cells that were unhit
536  if(unexpX+unexpY > 4) return 1e12;
537 
538  double xdens, ydens;
539  TotalChargePerView(zip, xdens, ydens);
540  xdens /= xTrk.NCell();
541  ydens /= yTrk.NCell();
542 
543  return unexpX*xdens+unexpY*ydens;
544  }
545 
546  //......................................................................
548  const rb::Track& xTrk,
549  const rb::Track& yTrk) const
550  {
551  // No overlap: MatchScore requires that anyway
552  if(xTrk.MinPlane() > yTrk.MaxPlane()) return 1e12;
553  if(yTrk.MinPlane() > xTrk.MaxPlane()) return 1e12;
554 
555  const double miss = MissingChargeByExtension(sliceMap, xTrk, yTrk);
556 
557  // End up passing through too many live cells that were unhit
558  if(miss > 5e11) return 1e12;
559 
560  // Disable charge difference component for now.
561  // const double diff = ChargeDifferenceBetweenViews(ZipTracks(xTrk, yTrk));
562 
563  const double score = MatchScore(miss/*+diff*/, xTrk, yTrk);
564 
565  LOG_DEBUG("ViewMerger") << "Missing: " << miss
566  << ", score = " << score;
567  // << ", charge difference " << diff << std::endl;
568 
569  return score;
570  }
571 
572  //......................................................................
574  const rb::Track& single,
575  const rb::Track& parta,
576  const rb::Track& partb) const
577  {
578  // Must be in this order, must not overlap
579  if(partb.MinPlane() < parta.MaxPlane()) return 1e12;
580 
581  // No overlap: MatchScore requires that anyway
582  if(parta.MinPlane() > single.MaxPlane()) return 1e12;
583  if(single.MinPlane() > partb.MaxPlane()) return 1e12;
584 
585  const rb::Track& join = JoinTracks(parta, partb);
586 
587  const rb::Track& zip = ZipTracks(single, join);
588  double totx, toty;
589  TotalChargePerView(zip, totx, toty);
590  const double avgDens = (totx+toty)/zip.NCell();
591 
592  // Don't allow join to use existing hits, only extension
593  const int unexpJoin = CountUnexplainedOnLine(rb::HitMap(),
594  parta.Stop(), partb.Start(),
595  join.View());
596  // Join shouldn't go through any unhit cells or others' cells. This
597  // provides some protection against crazy merges.
598  if(unexpJoin > 0) return 1e12;
599 
600  // Charge missing at the ends
601  const double miss = MissingChargeByExtension(sliceMap, single, join);
602 
603  // End up passing through too many live cells that were unhit
604  if(miss > 4*avgDens) return 1e12;
605 
606  const double diff = ChargeDifferenceBetweenViews(ZipTracks(single, join));
607 
608  // Small penalty for having a join
609  const double score = MatchScore(miss+avgDens+diff, single, join);
610 
611  LOG_DEBUG("ViewMerger") << "Join: missing charge " << miss
612  << ", score = " << score
613  << ", charge difference: " << diff << std::endl;
614 
615  return score;
616  }
617 
618  //......................................................................
620  const rb::Track& partb) const
621  {
622  assert(parta.View() == partb.View());
623 
624  // Must be in this order, must not overlap
625  assert(partb.MinPlane() >= parta.MaxPlane());
626 
627  // This is pretty crude. Could do a neater job with tighter interface into
628  // the 2D tracker.
629 
630  rb::Track join(parta);
631  // Right now, no two tracks can share hits, so this is safe. Might not be
632  // so in future, would have to remove duplicates.
633  join.Add(partb.AllCells());
634  for(unsigned int n = 0; n < partb.NTrajectoryPoints(); ++n){
635  join.AppendTrajectoryPoint(partb.TrajectoryPoint(n));
636  }
637 
638  return join;
639  }
640 
641  //......................................................................
643  TVector3 start, TVector3 end,
644  geo::View_t view) const
645  {
646  // TODO: we need to penalize the case where the track now runs outside the
647  // detector for some of its length. This is actually leading to
648  // mismerges...
651 
652  std::vector<geo::OfflineChan> xcells, ycells;
653  geom->CountCellsOnLineFast(start.X(), start.Y(), start.Z(),
654  end.X(), end.Y(), end.Z(),
655  xcells, ycells);
656 
657  int ret = 0;
658  if(view != geo::kY){ // X or 3D
659  for(unsigned int n = 0; n < xcells.size(); ++n){
660  const int plane = xcells[n].Plane();
661  const int cell = xcells[n].Cell();
662  if(!hmap.CellExists(plane, cell) &&
663  !bcl->IsBad(plane, cell) &&
664  !bcl->IsLowEfficiency(plane, cell)) ++ret;
665  }
666  }
667  if(view != geo::kX){ // Y or 3D
668  for(unsigned int n = 0; n < ycells.size(); ++n){
669  const int plane = ycells[n].Plane();
670  const int cell = ycells[n].Cell();
671  if(!hmap.CellExists(plane, cell) &&
672  !bcl->IsBad(plane, cell) &&
673  !bcl->IsLowEfficiency(plane, cell)) ++ret;
674  }
675  }
676 
677  return ret;
678  }
679 
680  //......................................................................
682  {
683  double totx, toty;
684  TotalChargePerView(zip, totx, toty);
685 
686  return fabs(totx-toty);
687  }
688 
689  //......................................................................
691  std::vector<rb::Track>& tracks) const
692  {
693  rb::HitMap hmap;
694  for(unsigned int n = 0; n < tracks.size(); ++n) hmap.Add(&tracks[n]);
695 
696  for(unsigned int n = 0; n < tracks.size(); ++n){
697  // Only if we don't pass through a single non-dead cell
698  if(CountUnexplainedOnLine(hmap, tracks[n].Start(), vtx->GetXYZ(),
699  tracks[n].View()) == 0){
700 
701  TVector3 xyz = vtx->GetXYZ();
702  if(!tracks[n].Is3D()) xyz[1-tracks[n].View()] = 0;
703 
704  // And only if the dot-product is reasonable
705  if((tracks[n].Start()-xyz).Unit().Dot(tracks[n].Dir()) < .94) continue;
706 
707  tracks[n].PrependTrajectoryPoint(xyz);
708  if(tracks[n].Is3D()){
709  // Overwrites the point we just prepended
710  tracks[n].SetStart(vtx->GetXYZ());
711  }
712  else{
713  // Overwrites the point we just prepended
714  tracks[n].SetStart(xyz[tracks[n].View()], xyz.Z());
715  }
716  }
717  } // end for n
718  }
719 
720  //......................................................................
721  double ViewMerger::MatchScore(double miss,
722  const rb::Track& aTrk,
723  const rb::Track& bTrk) const
724  {
725  // These casts to int are absolutely necessary :-@
726  const int dStart = abs(int(aTrk.MinPlane())-int(bTrk.MinPlane()));
727  assert(dStart > 0);
728  const int dEnd = abs(int(aTrk.MaxPlane())-int(bTrk.MaxPlane()));
729  assert(dEnd > 0);
730  const int totL = aTrk.ExtentPlane()+bTrk.ExtentPlane();
731  const int overlap = totL-dStart-dEnd;
732  if(overlap <= 0) return 1e12;
733 
734  // Add a small amount to miss so we can distinguish different zeroes
735  return (miss+1e-3)/overlap;
736  }
737 
738  //......................................................................
740  const rb::Track& yTrk) const
741  {
742  return xTrk.ZipWith(yTrk);
743 
744  // TODO: restore protection against going up to a plane too far. Can just
745  // truncate that returned track.
746 
747  /*
748  if(xTrk.View() == geo::kY){
749  assert(yTrk.View() == geo::kX);
750  return ZipTracks(yTrk, xTrk);
751  }
752 
753  art::ServiceHandle<geo::Geometry> geom;
754 
755  std::vector<TVector3> trajPts;
756 
757  double minZ = std::min(xTrk.Start().Z(), yTrk.Start().Z());
758  double maxZ = std::max(xTrk.Stop().Z(), yTrk.Stop().Z());
759 
760  const int minPlane = std::min(xTrk.MinPlane(), yTrk.MinPlane());
761  const int maxPlane = std::max(xTrk.MaxPlane(), yTrk.MaxPlane());
762 
763  // 2D tracks often stick halfway into the next plane. Don't allow the final
764  // track to go further than we have any evidence of hits. In that case put
765  // it to within 1mm of the end of the last plane with hits.
766  double xyz[3], dxyz[3];
767  geom->CellInfo(minPlane, 0, 0, xyz, dxyz);
768  minZ = std::max(minZ, xyz[2]-dxyz[2]+.1);
769 
770  geom->CellInfo(maxPlane, 0, 0, xyz, dxyz);
771  maxZ = std::min(maxZ, xyz[2]+dxyz[2]-.1);
772 
773  const double kZStep = 5; // cm
774  const int stepMax = int((maxZ-minZ)/kZStep);
775  for(int step = 0; step <= stepMax; ++step){
776  const double z = minZ+(double(step)/stepMax)*(maxZ-minZ);
777 
778  double xyz[3];
779  xyz[2] = z;
780 
781  double junk;
782  xTrk.InterpolateXY(z, xyz[0], junk);
783  yTrk.InterpolateXY(z, junk, xyz[1]);
784 
785  trajPts.push_back(xyz);
786  } // end for step
787 
788  rb::Track ret(art::PtrVector<rb::CellHit>(),
789  trajPts[0], trajPts[1]-trajPts[0]);
790  ret.Add(xTrk.AllCells());
791  ret.Add(yTrk.AllCells());
792 
793  for(unsigned int n = 1; n < trajPts.size(); ++n)
794  ret.AppendTrajectoryPoint(trajPts[n]);
795 
796  return ret;
797  */
798  }
799 
800  //......................................................................
802  std::map<int, int>& mins,
803  std::map<int, int>& maxs) const
804  {
805  for(unsigned int n = 0; n < trk.NCell(); ++n){
806  const int p = trk.Cell(n)->Plane();
807  const int c = trk.Cell(n)->Cell();
808 
809  if(c > maxs[p]) maxs[p] = c;
810  // Careful with zero default value
811  if(mins.find(p) == mins.end() || c < mins[p]) mins[p] = c;
812  }
813  }
814 
815  //......................................................................
816  bool ViewMerger::Adjacent(const rb::Track& a, const rb::Track& b) const
817  {
818  assert(a.View() == b.View());
819  assert(a.Is2D());
820 
821  if(a.ExtentPlane() < b.ExtentPlane()) return Adjacent(b, a);
822  // Guaranteed b is the shorter one
823 
824  // Don't attach verticals
825  if(b.ExtentPlane() == 1) return false;
826 
827  std::map<int, int> amins, amaxs, bmins, bmaxs;
828  Extremes(a, amins, amaxs);
829  Extremes(b, bmins, bmaxs);
830 
831  // Maybe a is all below b
832  bool all = true;
833  for(std::map<int, int>::iterator it = bmins.begin(); it != bmins.end(); ++it){
834  if(amaxs[it->first] != it->second-1){
835  all = false;
836  break;
837  }
838  }
839  if(all) return true;
840 
841  // Or maybe a is all above b
842  for(std::map<int, int>::iterator it = bmaxs.begin(); it != bmaxs.end(); ++it){
843  if(amins[it->first] != it->second+1) return false;
844  }
845  return true;
846  }
847 
848  //......................................................................
849  void ViewMerger::MergeAdjacentTracks(std::list<rb::Track>& trks) const
850  {
851  for(std::list<rb::Track>::iterator iti = trks.begin(); iti != trks.end(); ++iti){
852  for(std::list<rb::Track>::iterator itj = trks.begin(); itj != trks.end(); ++itj){
853  if(iti == itj) continue;
854 
855  if(iti->ExtentPlane() < itj->ExtentPlane()) continue;
856 
857  if(Adjacent(*iti, *itj)){
858  LOG_DEBUG("ViewMerger") << "Adjacent: " << *iti << " and " << *itj << std::endl;
859 
860  iti->Add(itj->AllCells());
861  trks.erase(itj);
862  MergeAdjacentTracks(trks); // Start again from the top
863  return;
864  }
865  } // end for itj
866  } // end for iti
867  }
868 
869  // Links against the version in DiscreteTracker
870  bool CompareByLengthVM(const rb::Track& a, const rb::Track& b)
871  {
872  return a.TotalLength() < b.TotalLength();
873  }
874 
875  //......................................................................
876  rb::Vertex* ViewMerger::FindVertex(std::vector<rb::Track> trks) const
877  {
878  if(trks.empty()) return 0;
879 
880  // Longest tracks first. Only used for fallback at the end
881  std::sort(trks.rbegin(), trks.rend(), CompareByLengthVM);
882 
883  double bestVx, bestVy, bestVz;
884  // I don't believe there's a path that doesn't set these, but need this to
885  // shut the compiler up in optimized builds.
886  bestVx = bestVy = bestVz = 0;
887 
888  double bestScore = 1e10;
889  for(unsigned int idx1 = 0; idx1 < trks.size()-1; ++idx1){
890  if(trks[idx1].View() != geo::kXorY) continue;
891  if(trks[idx1].ExtentPlane() < 3) continue;
892  for(bool end1: {false, true}){
893  for(unsigned int idx2 = idx1+1; idx2 < trks.size(); ++idx2){
894  if(trks[idx2].View() != geo::kXorY) continue;
895  if(trks[idx2].ExtentPlane() < 3) continue;
896  for(bool end2: {false, true}){
897 
898  // Start and direction defining point of these two lines
899  const TVector3 s1 = end1 ? trks[idx1].Stop() : trks[idx1].Start();
900  const TVector3 s2 = end2 ? trks[idx2].Stop() : trks[idx2].Start();
901  const TVector3 d1 = trks[idx1].TrajectoryPoint(end1 ? trks[idx1].NTrajectoryPoints()-2 : 1);
902  const TVector3 d2 = trks[idx2].TrajectoryPoint(end2 ? trks[idx2].NTrajectoryPoints()-2 : 1);
903 
904  double vx, vy, vz1, vz2;
905  bool ok = geo::LineIntersection(s1.X(), s1.Z(), d1.X(), d1.Z(),
906  s2.X(), s2.Z(), d2.X(), d2.Z(),
907  vx, vz1);
908  if(!ok) continue;
909  ok = geo::LineIntersection(s1.Y(), s1.Z(), d1.Y(), d1.Z(),
910  s2.Y(), s2.Z(), d2.Y(), d2.Z(),
911  vy, vz2);
912  if(!ok) continue;
913  const double vz = (vz1+vz2)/2;
914 
915  const double dist1 = util::pythag(s1.X()-vx, s1.Y()-vy, s1.Z()-vz);
916  const double dist2 = util::pythag(s2.X()-vx, s2.Y()-vy, s2.Z()-vz);
917 
918  // Favour both long tracks and hits close to the ends
919  const double score = dist1/trks[idx1].TotalLength()+dist2/trks[idx2].TotalLength();
920  if(score < bestScore){
921  bestScore = score;
922  bestVx = vx;
923  bestVy = vy;
924  bestVz = vz;
925  } // end if better
926  } // end for end2
927  } // end for idx2
928  } // end for end1
929  } // end for idx1
930  // Too far from the ends, fallback
931  if(bestScore > 6){
932  bestVx = trks[0].Start().X();
933  bestVy = trks[0].Start().Y();
934  bestVz = trks[0].Start().Z();
935  }
936 
937  // Time is an afterthought
938  return new rb::Vertex(bestVx, bestVy, bestVz, trks[0].MinTNS());
939  }
940 
941  //......................................................................
943  const std::vector<rb::Track>& merged) const
944  {
945  assert(!trk.Is3D());
946 
947  const double z = trk.Start().Z();
948 
949  const rb::Track* best = 0;
950 
951  // The best source for W information is the longest overlapping track
952  for(unsigned int n = 0; n < merged.size(); ++n){
953  if(merged[n].Start().Z()-5 < z &&
954  merged[n].Stop().Z()+5 > z){
955  if(!best || merged[n].ExtentPlane() > best->ExtentPlane())
956  best = &merged[n];
957  }
958  }
959 
960  // Failing that, just the longest track
961  if(!best){
962  for(unsigned int n = 0; n < merged.size(); ++n){
963  if(!best || merged[n].ExtentPlane() > best->ExtentPlane())
964  best = &merged[n];
965  }
966  }
967 
968  const geo::View_t view = trk.View();
969 
970  double w = 0;
971  if(best){
972  double x, y;
973  best->InterpolateXY(z, x, y);
974  w = (view == geo::kX) ? y : x;
975  }
976 
977  // Fill in W information everywhere
978 
979  TVector3 start = trk.Start();
980  start[1-view] = w;
981 
982  rb::Track ret(trk.AllCells(), start, trk.Dir());
983 
984  for(unsigned int n = 1; n < trk.NTrajectoryPoints(); ++n){
985  TVector3 pt = trk.TrajectoryPoint(n);
986  pt[1-view] = w;
987  ret.AppendTrajectoryPoint(pt);
988  }
989 
990  return ret;
991  }
992 
994 
995 } //namespace dt
996 ////////////////////////////////////////////////////////////////////////
virtual void endJob()
virtual geo::View_t View() const
kXorY for 3D clusters.
Definition: Cluster.h:99
#define LOG_DEBUG(stream)
Definition: Messenger.h:149
size_t NTrajectoryPoints() const
Definition: Track.h:83
rb::Track ZipVerticalTracks(const rb::Track &xtrk, const rb::Track &ytrk) const
Tilt tracks arbitrarily and merge them into a 3D track.
A 3D position and time representing an interaction vertex.
Definition: Vertex.h:15
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 MergeAdjacentTracks(std::list< rb::Track > &trks) const
Match 2D tracks into 3D tracks.
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
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
set< int >::iterator it
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
virtual void produce(art::Event &evt)
float Dot(const Proxy &v) const
Float_t y1[n_points_granero]
Definition: compare.C:5
rb::Track ZipTracks(const rb::Track &aTrk, const rb::Track &bTrk) const
Make a 3D track based on two 2D tracks.
X or Y views.
Definition: PlaneGeo.h:30
bool Is2D() const
Definition: Cluster.h:96
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
unsigned short Plane() const
Definition: CellHit.h:39
Float_t x1[n_points_granero]
Definition: compare.C:5
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
float Z() const
Definition: RecoHit.h:38
std::string fTrackInput
Module name to look for 2D tracks under.
void Extremes(const rb::Track &trk, std::map< int, int > &mins, std::map< int, int > &maxs) const
unsigned int MaxCell(geo::View_t view) const
Definition: Cluster.cxx:518
void abs(TH1 *hist)
bool progress
Insert implicit nodes #####.
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
DEFINE_ART_MODULE(TestTMapFile)
virtual TVector3 Start() const
Definition: Prong.h:73
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
void CellInfo(unsigned int ip, unsigned int ic, View_t *view=0, double *pos=0, double *dpos=0) const
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
void CountCellsOnLineFast(double x1, double y1, double z1, double x2, double y2, double z2, std::vector< OfflineChan > &Xhitsonline, std::vector< OfflineChan > &Yhitsonline)
Make list of cells in each view traversed by a line segment.
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
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
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
fclose(fg1)
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
Definition: Cand.cxx:23
void TotalChargePerView(const rb::Track &zip, double &totx, double &toty) const
double dz[NP][NC]
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
const double a
TVector3 GetXYZ() const
Definition: Vertex.cxx:45
T get(std::string const &key) const
Definition: ParameterSet.h:231
rb::Track FakeThirdDimension(const rb::Track &trk, const std::vector< rb::Track > &merged) const
Turn a 2D track into a 3D track.
rb::Track ZipWith(const rb::Track &trk) const
Combine with a 2D track from the other view to make a 3D track.
Definition: Track.cxx:442
virtual TVector3 Dir() const
Unit vector describing prong direction.
Definition: Prong.h:77
int evt
rb::RecoHit RecoHit(const art::Ptr< rb::CellHit > &chit) const
Return calibrated hit based on assumed W coordinate.
Definition: Cluster.cxx:259
ViewMerger(const fhicl::ParameterSet &pset)
Collect Geo headers and supply basic geometry functions.
bool IsCalibrated() const
You MUST check here before accessing PECorr, MIP or GeV.
Definition: RecoHit.cxx:35
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
double ScoreForExtension(const rb::HitMap &sliceMap, const rb::Track &xTrk, const rb::Track &yTrk) const
bool Adjacent(const rb::Track &a, const rb::Track &b) const
EventNumber_t event() const
Definition: Event.h:67
std::string fSliceInput
Module name to look for slices under.
TVector3 Unit() const
Vertex location in position and time.
Definition: View.py:1
z
Definition: test.py:28
size_type size() const
Definition: PtrVector.h:308
bool CompareByLengthVM(const rb::Track &a, const rb::Track &b)
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
float GeV() const
Definition: RecoHit.cxx:69
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double ScoreForJoinPlusExtension(const rb::HitMap &sliceMap, const rb::Track &single, const rb::Track &parta, const rb::Track &partb) const
double pythag(double x, double y)
2D Euclidean distance
Definition: MathUtil.h:29
bool fObeyPreselection
Check rb::IsFiltered?
std::vector< rb::Track > MatchVerticals(std::list< rb::Track > &xtracks, std::list< rb::Track > &ytracks) const
Find vertical tracks seperated by one plane. Return results of merge.
double MatchScore(double unexp, const rb::Track &aTrk, const rb::Track &bTrk) const
Find number of overlapping planes, return unexp/overlap.
unsigned int MinCell(geo::View_t view) const
Definition: Cluster.cxx:472
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)
void MaybeExtendToVertex(const rb::Vertex *vtx, std::vector< rb::Track > &tracks) const
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
rb::Track JoinTracks(const rb::Track &parta, const rb::Track &partb) const
Append two tracks in the same view (crudely).
virtual void reconfigure(const fhicl::ParameterSet &pset)
virtual void InterpolateXY(double z, double &x, double &y) const
Definition: Track.cxx:325
int CountUnexplainedOnLine(const rb::HitMap &hmap, TVector3 start, TVector3 end, geo::View_t view) const
hmap should contain cells that are hit and so OK to pass through
bool CellExists(unsigned int planeIdx, unsigned int cellIdx) const
Does the map contain any cell at this position?
Definition: HitMap.cxx:81
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Float_t e
Definition: plot.C:35
virtual void beginJob()
bool IsBad(int plane, int cell)
Float_t w
Definition: plot.C:20
std::vector< rb::Track > MatchTracks(const rb::HitMap &sliceMap, std::list< rb::Track > &xtracks, std::list< rb::Track > &ytracks) const
Do most of the work. Called from produce.
rb::Vertex * FindVertex(std::vector< rb::Track > trks) const
Toy vertex finder.
double ChargeDifferenceBetweenViews(const rb::Track &zip) const
Encapsulate the geometry of one entire detector (near, far, ndos)
virtual bool Is3D() const
Definition: Prong.h:71
bool IsLowEfficiency(int plane, int cell)
void Add(const art::Ptr< rb::CellHit > &chit)
Definition: HitMap.cxx:37
double MissingChargeByExtension(const rb::HitMap &sliceMap, const rb::Track &aTrk, const rb::Track &bTrk) const