WindowTrackingAlg.cxx
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////
2 /// \file WindowTrackingAlg.cxx
3 /// \brief Track finder for cosmic rays
4 /// \author brebel@fnal.gov
5 /// \date
6 ///////////////////////////////////////////////////////////////////////
7 #include <cmath>
8 
9 // NOvASoft includes
10 #include "RecoBase/Track.h"
12 #include "GeometryObjects/Geo.h"
13 #include "Utilities/func/MathUtil.h"
14 #include "RecoBase/CellHit.h"
15 #include "RecoBase/Cluster.h"
17 
18 // ART includes
19 #include "fhiclcpp/ParameterSet.h"
21 #include "cetlib_except/exception.h"
23 
24 #include "TStopwatch.h"
25 
26 bool gt(unsigned short int a, unsigned short int b)
27 {
28  return a > b;
29 }
30 
31 bool lt(unsigned short int a, unsigned short int b)
32 {
33  return a < b;
34 }
35 
36 bool gtvec(std::pair<geo::View_t, TVector3> a, std::pair<geo::View_t, TVector3> b)
37 {
38  return a.second.Z() > b.second.Z();
39 }
40 bool ltvec(std::pair<geo::View_t, TVector3> a, std::pair<geo::View_t, TVector3> b)
41 {
42  return a.second.Z() < b.second.Z();
43 }
44 
46 {
47  return a->Cell() < b->Cell();
48 }
49 
50 
51 namespace trk {
52 
53  //----------------------------------------------------------------------------
54  // Constructor
56  : TrackAlg(pset)
57  , fCellHalfDepth(-1.)
58  {
59  fDHitGood = pset.get<double >("DHitGood" );
60  fMinViewPlanes = pset.get<unsigned int>("MinViewPlanes" );
61  fWindowSize = pset.get<unsigned int>("WindowSize" );
62  fMinHitsTryAgain = pset.get<unsigned int>("MinHitsToTryAgain" );
63  fMinHitFraction = pset.get<double >("MinHitFraction" );
64  fNTrajectory = pset.get<unsigned int>("NTrajectoryForDir" );
65  fMaxPlaneSeparation = pset.get<int >("MaxPlaneSeparation");
66  }
67 
68  //----------------------------------------------------------------------------
69  // Destructor
71 
72  //----------------------------------------------------------------------------
73  // MakeTrack
74  std::vector<rb::Track> WindowTrackingAlg::MakeTrack( const rb::Cluster* slice )
75  {
76  if(fCellHalfDepth < 0.) fCellHalfDepth = fGeo->Plane(0)->Cell(0)->HalfD();
77 
78 
79  //LOG_DEBUG("WindowTrackingAlg") << "fitting new slice";
80 
81  std::vector<rb::Track> tracks;
82 
83  // get the PtrVector of all cells in this slice
84  art::PtrVector<rb::CellHit> allCells = slice->AllCells();
85  allCells.sort();
86  std::vector<art::Ptr<rb::CellHit> > allHitsFor;
87  std::vector<art::Ptr<rb::CellHit> > allHitsRev;
88 
89  for(auto cell : allCells){
90  allHitsFor.push_back(cell);
91  allHitsRev.push_back(cell);
92  }
93 
94  std::vector<art::Ptr<rb::CellHit> > removedHits;
95  std::vector<art::Ptr<rb::CellHit> > removedHitsFor;
96  std::vector<art::Ptr<rb::CellHit> > removedHitsRev;
97  bool tryAnother = true;
98 
99  double chiSqrFor = 0.;
100  double chiSqrRev = 0.;
101 
103 
104  while(tryAnother){
105  tryAnother = false;
106 
107  rb::Track trackForward = this->Make3DTrack(allHitsFor, removedHitsFor, chiSqrFor, lt);
108  rb::Track trackReverse = this->Make3DTrack(allHitsRev, removedHitsRev, chiSqrRev, gt);
109 
110  //LOG_DEBUG("WindowTrackingAlg") << "forward chi^2: " << chiSqrFor << "/" << trackForward.NCell()
111  // << "\nreverse chi^2: " << chiSqrRev << "/" << trackReverse.NCell();
112 
113  if(trackForward.NCell() > 0) chiSqrFor /= trackForward.NCell();
114  if(trackReverse.NCell() > 0) chiSqrRev /= trackReverse.NCell();
115 
116  if(chiSqrFor < chiSqrRev){
117  track = trackForward;
118  removedHits.swap(removedHitsFor);
119  }
120  else{
121  track = trackReverse;
122  removedHits.swap(removedHitsRev);
123  }
124 
125  removedHitsRev.clear();
126  removedHitsFor.clear();
127 
128  // check to see if we made a legit track, it will have at least 1 cell
129  if(track.NCell() > 1){
130  this->SetTrackEndPoints(track);
131 
132  tracks.push_back(track);
133 
134  /* Some debugging code commented out below, leaving it in place for future convenience
135  if(track.Dir().Z() > 0){
136  // get the cells for this track and sort them
137  art::PtrVector<rb::CellHit> trCells = track.AllCells();
138  trCells.sort();
139 
140  double xyzStart[3] = {0.};
141  double xyzEnd[3] = {0.};
142  fGeo->Plane(trCells[0]->Plane())->Cell(trCells[0]->Cell())->GetCenter(xyzStart);
143  fGeo->Plane(trCells[track.NCell()-1]->Plane())->Cell(trCells[track.NCell()-1]->Cell())->GetCenter(xyzEnd);
144 
145  LOG_VERBATIM("WindowTrackingAlg")
146  << "first cell " << xyzStart[2] //<< "," << *(trCells[0])
147  << "\nlast cell " << xyzEnd[2] //<< "," << *(trCells[track.NCell()-1])
148  << "\nstart: " << track.Start().X() << " " << track.Start().Y() << " " << track.Start().Z()
149  << "\nend: " << track.Stop().X() << " " << track.Stop().Y() << " " << track.Stop().Z();
150 
151  for(size_t h = 0; h < track.NCell(); ++h)
152  LOG_VERBATIM("WindowTrackingAlg") << *(track.Cell(h));
153  }*/
154 
155  if(removedHits.size() > fMinHitsTryAgain){
156  tryAnother = true;
157  //LOG_DEBUG("WindowTrackingAlg")
158  //<< "trying to make another track from the "
159  //<< removedHits.size()
160  //<< " removed hits";
161  allHitsFor.clear();
162  allHitsRev.clear();
163  for(auto removed : removedHits){
164  //LOG_DEBUG("WindowTrackingAlg")
165  //<< "removed "
166  //<< *removed;
167  allHitsFor.push_back(removed);
168  allHitsRev.push_back(removed);
169  }
170  removedHits.clear();
171  }// end if there are enough left over hits to look for another track
172  }// end if the returned track was valid
173 
174  }// end loop to look for tracks in the slice
175 
176  // check to see if any of the tracks should be merged
177  if(tracks.size() > 1){
178  //LOG_VERBATIM("WindowTrackingAlg")
179  //<< "there are "
180  //<< tracks.size()
181  //<< " tracks in the slice";
182  TVector3 t1Start, t2Start, t1End, t2End, diff12, diff21;
183  float t1ZDir = 0., t2ZDir = 0.;
184  size_t trks = tracks.size();
185  std::vector<bool> removed(trks, false);
186  for(size_t tr1 = 0; tr1 < trks; ++tr1){
187  // int ctr = 0;
188  if(removed[tr1]) continue;
189  for(size_t tr2 = 0; tr2 < trks; ++tr2){
190  if(tr1 == tr2 || removed[tr2]) continue;
191 
192  // LOG_VERBATIM("WindowTrackingAlg") << "inner loop " << ctr;
193  // ++ctr;
194 
195  // get the start and end points of each track
196  t1Start = tracks[tr1].Start();
197  t2Start = tracks[tr2].Start();
198  t1End = tracks[tr1].Stop();
199  t2End = tracks[tr2].Stop();
200 
201  t1ZDir = tracks[tr1].Dir().Z()/std::abs(tracks[tr1].Dir().Z());
202  t2ZDir = tracks[tr2].Dir().Z()/std::abs(tracks[tr2].Dir().Z());
203 
204  // LOG_VERBATIM("WindowTrackingAlg") << "t1 start: "
205  // << t1Start.X() << " " << t1Start.Y() << " " << t1Start.Z() << " end: "
206  // << t1End.X() << " " << t1End.Y() << " " << t1End.Z() << "\n"
207  // << "t2 start: "
208  // << t2Start.X() << " " << t2Start.Y() << " " << t2Start.Z() << " end: "
209  // << t2End.X() << " " << t2End.Y() << " " << t2End.Z();
210 
211  // now see if any of the end points match up
212  // that is they are within a cell of each other
213  // The total difference in the end points has to be
214  // less than 2*fDHitGood and the z component of that
215  // should be less than fDHitGood
216  // also make sure they are going in the same z direction
217  if(t1ZDir*t2ZDir > 0.){
218  diff12 = t1Start - t2End;
219  diff21 = t2Start - t1End;
220 
221  if( (diff12.Mag() < 2.*fDHitGood && std::abs(diff12.Z()) < fDHitGood) ||
222  (diff21.Mag() < 2.*fDHitGood && std::abs(diff21.Z()) < fDHitGood) ){
223 
224  // get the iterators in the correct order for the z direction of
225  // the track
226  auto firstTr = tracks[tr1];
227  auto secondTr = tracks[tr2];
228  if(t1Start.Z()*t1ZDir > t2Start.Z()*t2ZDir){
229  firstTr = tracks[tr2];
230  secondTr = tracks[tr1];
231  }
232 
233  // merge the tracks
234  std::vector<TVector3> mergedTrajectory;
235  // get a merged trajectory point for the intersection of the two tracks
236  // and set that to the be the first trajectory point in track 1, don't
237  // prepend last point from the trajectory of track 2
238  // the first two points and last two points in a track represent the edge
239  // of and the midpoint of the final plane, so remove the last two points
240  // in firstTr and the first two points in secondTr from the merged track
241  for(size_t tp = 0; tp < firstTr.Trajectory().size()-2; ++tp){
242  mergedTrajectory.push_back(firstTr.TrajectoryPoint(tp));
243  //LOG_VERBATIM("WindowTrackingAlg")
244  //<< mergedTrajectory.back().X() << " "
245  //<< mergedTrajectory.back().Y() << " "
246  //<< mergedTrajectory.back().Z();
247  }
248  mergedTrajectory.push_back(0.5*(secondTr.Trajectory().front() + firstTr.Trajectory().back()));
249  //LOG_VERBATIM("WindowTrackingAlg") << "connection point "
250  //<< mergedTrajectory.back().X() << " "
251  //<< mergedTrajectory.back().Y() << " "
252  //<< mergedTrajectory.back().Z();
253  for(size_t tp = 2; tp < secondTr.Trajectory().size(); ++tp){
254  mergedTrajectory.push_back(secondTr.TrajectoryPoint(tp));
255  //LOG_VERBATIM("WindowTrackingAlg")
256  //<< mergedTrajectory.back().X() << " "
257  //<< mergedTrajectory.back().Y() << " "
258  //<< mergedTrajectory.back().Z();
259  }
260 
261  // now create the merged track with cells from both tracks
262  rb::Track mergedTrack(firstTr.AllCells(), mergedTrajectory.front(), firstTr.Dir());
263  for(auto c : secondTr.AllCells()) mergedTrack.Add(c);
264  for(size_t tp = 1; tp < mergedTrajectory.size(); ++tp) mergedTrack.AppendTrajectoryPoint(mergedTrajectory[tp]);
265  tracks.push_back(mergedTrack);
266 
267  //LOG_VERBATIM("WindowTrackingAlg") << firstTr.Trajectory().size() << " "
268  //<< secondTr.Trajectory().size() << " "
269  //<< mergedTrajectory.size() << " trajectory points now: \n";
270  //for(auto tp : tracks.back().Trajectory())
271  // LOG_VERBATIM("WindowTrackingAlg") << tp.X() << " " << tp.Y() << " " << tp.Z();
272 
273  // remove both track 1 and 2
274  removed[tr1] = true;
275  removed[tr2] = true;
276  //LOG_VERBATIM("WindowTrackingAlg") << "removing tracks " << tr1 << " " << tr2;
277  }// end if the tracks end points are close enough to merge
278  }// end if tracks are going in the same direction in Z
279  }// end inner loop over tracks
280  }// end outer loop over tracks
281 
282  // loop over the removed vector and remove those tracks from the vector
283  size_t rmctr = 0;
284  for(size_t r = 0; r < removed.size(); ++r){
285  if(removed[r]){
286  tracks.erase(tracks.begin() + (r-rmctr));
287  ++rmctr;
288  }
289  }
290  //LOG_VERBATIM("WindowTrackingAlg") << "there are " << tracks.size() << " tracks in the slice";
291  }// end if more than one track in the slice
292 
293  return tracks;
294 
295  } // End of MakeTrack
296 
297  //......................................................................
299  std::vector<art::Ptr<rb::CellHit> > & removedHits,
300  double & chiSqr,
301  bool(*comp)(unsigned short int, unsigned short int))
302  {
303  // loop over all the hits for this slice and get the
304  // x or y cells
305  std::map<unsigned short int, std::vector< art::Ptr<rb::CellHit> >,
306  bool(*)(unsigned short int, unsigned short int) > xViewHits(comp);
307  std::map<unsigned short int, std::vector< art::Ptr<rb::CellHit> >,
308  bool(*)(unsigned short int, unsigned short int) > yViewHits(comp);
309  std::map<unsigned short int, std::pair<float, float>,
310  bool(*)(unsigned short int, unsigned short int) > planeCoG(comp);
311 
312  double xyz[3] = {0.};
313  double v = 0.;
314 
315  // the allCells vector is already sorted properly.
316  std::set<int> planes;
317 
318  //LOG_VERBATIM("WindowTrackingAlg") << "starting with these hits:";
319  for(auto const chit : allCells){
320 
321  planes.insert(chit->Plane());
322 
323  fGeo->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
324 
325  //LOG_VERBATIM("WindowTrackingAlg") << "\t" << chit->Plane() << " " << chit->Cell();
326 
327  if(chit->View() == geo::kX){
328  xViewHits[chit->Plane()].push_back(chit);
329  v = xyz[0];
330  }
331  if(chit->View() == geo::kY){
332  yViewHits[chit->Plane()].push_back(chit);
333  v = xyz[1];
334  }
335 
336  if(planeCoG.find(chit->Plane()) != planeCoG.end()){
337  planeCoG[chit->Plane()].first += chit->ADC()*v;
338  planeCoG[chit->Plane()].second += chit->ADC();
339  }
340  else
341  planeCoG[chit->Plane()] = std::pair<float,float>(chit->ADC()*v,
342  chit->ADC());
343 
344  }// end loop over hits in this slice
345 
346  // loop over planeCoG and see what we have
347  //for(auto pcog : planeCoG)
348  // LOG_VERBATIM("WindowTrackingAlg")
349  // << pcog.first
350  // << " "
351  // << pcog.second.first/pcog.second.second;
352 
353  // remove planes that are too far separated from others in the slice
354  auto prev = planes.begin();
355  auto next = prev;
356  bool removeThisPlane = false;
357 
358  for(auto itr = planes.begin(); itr != planes.end(); ++itr){
359  next = itr; ++next;
360  removeThisPlane = false;
361 
362  // is the first plane separated from the second by too much?
363  if(itr == planes.begin() && next != planes.end()){
364  if(std::abs(*itr - *next) > fMaxPlaneSeparation) removeThisPlane = true;
365  }// end if on the first plane in the set
366  // is the current plane separated by both the previous and next one
367  // too much
368  else if(itr != planes.begin() && next != planes.end()){
369  prev = itr; --prev;
370  if(std::abs(*itr - *prev) > fMaxPlaneSeparation &&
371  std::abs(*itr - *next) > fMaxPlaneSeparation) removeThisPlane = true;
372  }
373  else if(next == planes.end()){
374  prev = itr; --prev;
375  if(std::abs(*itr - *prev) > fMaxPlaneSeparation) removeThisPlane = true;
376  }
377 
378  if(removeThisPlane){
379  xViewHits.erase(*itr);
380  yViewHits.erase(*itr);
381  planeCoG .erase(*itr);
382  // LOG_VERBATIM("WindowTrackingAlg") << "removing plane " << *itr << " from the fit; prev: "
383  // << *prev << " next: " << *next;
384  }
385 
386  }//end loop to look for big gaps between planes.
387 
388  rb::Track retTrack;
389 
390  if(this->ValidTrack(xViewHits) &&
391  this->ValidTrack(yViewHits) ){
392 
393  double chiSqrX = 0.;
394  double chiSqrY = 0.;
395 
396  rb::Track xTrack = FitView(xViewHits, planeCoG, removedHits, geo::kX, chiSqrX);
397  rb::Track yTrack = FitView(yViewHits, planeCoG, removedHits, geo::kY, chiSqrY);
398 
399  chiSqr = chiSqrX + chiSqrY;
400 
401  if(xTrack.View() != geo::kX || yTrack.View() != geo::kY){
402  LOG_DEBUG("WindowTrackingAlg") << "no track found in at least one view, return empty track";
403  return retTrack;
404  }
405 
406  //LOG_DEBUG("WindowTrackingAlg") << " x track view: " << xTrack.View()
407  // << " y track view: " << yTrack.View();
408 
409  // try making a 3D track from these two 2D tracks
410  try{
411  retTrack = xTrack.ZipWith(yTrack);
412  }
413  catch(cet::exception &e){
414  LOG_WARNING("WindowTrackingAlg") << "attempt to zip tracks failed with exception\n"
415  << e
416  << "\n return empty track";
417  return retTrack;
418  }
419 
420  TVector3 start(retTrack.Start());
421  TVector3 end (retTrack.Stop());
422 
423  this->CheckTrackDirectionInY(retTrack, start, end);
424 
425  // determine the direction of the track by fitting the first
426  // several trajectory points - just using the direction from the
427  // first point to the next can be misleading - need a bigger lever arm
428  this->DetermineInitialDirection(retTrack);
429 
430  // restrict track to exist inside detector; if this isn't possible,
431  // the track is no good, so don't write
433  //std::cout << "-------------------------------" << std::endl;
434  //for(auto const& traj : retTrack.Trajectory()) traj.Print();
435  return retTrack;
436  }
437 
438  } // end if we can make a 3D track
439  // see if we can make a short track - only worry about tracks that have at most fMinViewPlanes
440  // in each view - if one view has 12 planes and the other 2, it is probably garbage anyway
441  else if(xViewHits.size() < fMinViewPlanes &&
442  yViewHits.size() < fMinViewPlanes &&
443  xViewHits.size() > 0 &&
444  yViewHits.size() > 0 ){
445  //LOG_VERBATIM("WindowTrackingAlg")
446  //<< "attempting to fit a short track: "
447  //<< xViewHits.size() << " " << yViewHits.size();
448 
449  return this->ShortTrack(xViewHits, yViewHits, planeCoG, chiSqr);
450  }
451  else
452  // can't make a track, so set the chiSqr value to something stupidly big
453  chiSqr = 1.e9;
454 
455  return retTrack;
456  }
457 
458  //......................................................................
459  rb::Track WindowTrackingAlg::FitView(std::map<unsigned short int, std::vector<art::Ptr<rb::CellHit> >,
460  bool(*)(unsigned short int, unsigned short int) > & hits,
461  std::map<unsigned short int, std::pair<float, float>,
462  bool(*)(unsigned short int, unsigned short int) > const& planeCoG,
463  std::vector<art::Ptr<rb::CellHit> > & removedHits,
465  double & chiSqr)
466  {
467  // start off with an N plane window, then step out one plane at a time to add to the window
468  size_t numPlanes = hits.size();
469 
470  // another vector that has all the hits for a given window
471  std::vector< art::Ptr<rb::CellHit> > currWindow;
472  std::vector< art::Ptr<rb::CellHit> > prevWindow;
473 
474  // get the iterator to the hits of the first plane in the window
475  auto startItr = hits.begin();
476 
477  // get the iterator the hits of the last plane in the window
478  auto endItr = hits.begin();
479 
480  if(numPlanes > fWindowSize) std::advance(endItr, fWindowSize);
481  else endItr = hits.end();
482 
483  bool dir = (startItr->first > endItr->first);
484 
485  // fit the initial window
486  unsigned short int startPlane = startItr->second.front()->Plane();
487  for(auto itr = startItr; itr != endItr; ++itr){
488  for(size_t h = 0; h < itr->second.size(); ++h){
489  currWindow.push_back(itr->second[h]);
490  // LOG_VERBATIM("WindowTrackAlg") << "adding to initial window: " << itr->second[h]->Plane() << " " << itr->second[h]->Cell();
491  }
492  }
493 
494  // this is the fit of the initial window, there is nothing in prevWindow
495  // if all the hits are removed from currWindow, return an empty track
496  this->FitWindow(prevWindow, currWindow, removedHits, planeCoG, view);
497  if(currWindow.size() < 1){
498  //LOG_DEBUG("WindowTrack") << "currWindow size is " << currWindow.size()
499  // << " returning empty track";
500  return rb::Track();
501  }
502 
503  // put the current window's hits into the vector we'll use
504  // to make the cluster
505  std::vector< art::Ptr<rb::CellHit> > clHits(currWindow);
506 
507  // if there are no more planes to look at, just make the track
508  // otherwise, endItr is already on the element past the initial
509  // window
510  if(endItr == hits.end() ){
511  // make the track from the window
512  return this->MakeViewTrack(clHits, view, chiSqr, dir, removedHits);
513  }
514 
515  prevWindow.swap(currWindow);
516 
517  //LOG_DEBUG("WindowTrackingAlg") << "start sliding window";
518 
519  // slide the window until you run out of planes
520  while(endItr != hits.end()){
521 
522  //LOG_DEBUG("WindowTrackingAlg") << "prev window size now " << prevWindow.size()
523  // << " curr window has size " << endItr->second.size()
524  // << " starts with plane " << endItr->second.front()->Plane();
525 
526  // remove the first plane of hits from the previousWindow, but only if
527  // we added hits from the current window the last time
528  auto itr = prevWindow.begin();
529  while( itr != prevWindow.end() && currWindow.size() > 0){
530  //LOG_DEBUG("WindowTrackingAlg") << (*itr)->Plane() << " " << (*itr)->Cell();
531  if((*itr)->Plane() == startPlane) itr = prevWindow.erase(itr);
532  else ++itr;
533  }
534  startPlane = prevWindow.front()->Plane();
535 
536  // put the current plane's hits to the the currWindow vector
537  currWindow.swap(endItr->second);
538 
539  this->FitWindow(prevWindow, currWindow, removedHits, planeCoG, view);
540 
541  // add the surviving currWindow hits to the prevWindow and to the clHits
542  for(size_t h = 0; h < currWindow.size(); ++h){
543  prevWindow.push_back(currWindow[h]);
544  clHits.push_back(currWindow[h]);
545  }
546 
547  ++endItr;
548  }// end loop over the windows
549 
550  // make the track using the hits that were found along the
551  // trajectory
552  return this->MakeViewTrack(clHits, view, chiSqr, dir, removedHits);
553  }
554 
555  //......................................................................
557  geo::View_t const& view,
558  double & chiSqr,
559  bool const& dir,
560  std::vector<art::Ptr<rb::CellHit> > & removedHits)
561  {
562  //LOG_DEBUG("WindowTrackingAlg")
563  //<< "view: " << view
564  //<< " has " << trkHits.size()
565  //<< " cell hits";
566  rb::Track retTrack;
567 
568  // make some trajectory points for this track - use the signal weighted
569  // average point in each plane
570 
571  // first make a map of the planes to hits and then find the weighted average for
572  // each plane, set the direction in plane number based on the input bool
573  bool(*comp)(unsigned short int, unsigned short int) = (dir) ? gt : lt;
574  std::map<unsigned short int, std::vector<art::Ptr<rb::CellHit> >,
575  bool(*)(unsigned short int, unsigned short int) > planesToHitLoc(comp);
576 
577  for(auto const hititr : trkHits)
578  planesToHitLoc[hititr->Plane()].push_back(hititr);
579 
580  // check that you ended up with a reasonable track in terms
581  // of the number of hit planes compared to the number of planes crossed
582  // set the chiSqr to something stupidly huge if that is not the case and
583  // return an empty track
584  if( !this->ValidTrack(planesToHitLoc) ){
585  chiSqr = 1.e6;
586  return retTrack;
587  }
588 
589  std::vector<std::pair<double, double> > traj;
590  double xyz[3] = {0.};
591  double &w = (view == geo::kX) ? xyz[0] : xyz[1];
592  double avg = 0.;
593  chiSqr = 0.;
594  unsigned short int plane = 0;
595 
596  for(auto const mapitr : planesToHitLoc){
597  double weightedCells = 0.;
598  double weights = 0.;
599  plane = mapitr.first;
600  for(auto const cellitr : mapitr.second){
601  //LOG_VERBATIM("WindowTrackingAlg")
602  //<< "makeviewtrack: "
603  //<< cellitr->Plane()
604  //<< " " << cellitr->Cell()
605  //<< " " << cellitr->ADC();
606  fGeo->Plane(plane)->Cell(cellitr->Cell())->GetCenter(xyz);
607  weightedCells += w*cellitr->ADC();
608  weights += 1.*cellitr->ADC();
609 
610  } //end loop over cell hits in the current plane
611 
612  // the z position is the same for all cells in a plane, and we are looping
613  // over all cells in each plane
614  if(weights > 0.){
615  avg = weightedCells/weights;
616 
617  for(auto const cellitr : mapitr.second){
618  fGeo->Plane(plane)->Cell(cellitr->Cell())->GetCenter(xyz);
619  chiSqr += std::pow(avg - w, 2.)/3.8; // cell height is 3.8 cm
620  }
621 
622  // LOG_DEBUG("WindowTrackingAlg") << "adding trajectory point in "
623  // << view << " view: " << xyz[2] << " " << avg;
624 
625  traj.push_back(std::pair<double, double>(avg, xyz[2]));
626  }
627 
628  } // end loop over planesToHitLoc
629 
630 
631  double dz = traj[1].second - traj[0].second;
632  double dv = traj[1].first - traj[0].first;
633 
634  //LOG_DEBUG("WindowTrackingAlg") << "there are " << traj.size() << " trajectory points in this view";
635 
636  retTrack = rb::Track(trkHits, view, traj[0].first, traj[0].second, dv, dz);
637  for(size_t t = 1; t < traj.size(); ++t)
638  retTrack.AppendTrajectoryPoint(traj[t].first, traj[t].second);
639 
640  // look for bremstrahlung hits to add back in
641  this->LookForBremsstrahlungHits(retTrack, removedHits);
642 
643  return retTrack;
644  }
645 
646  //......................................................................
647  // this method fits a straight line to the hits it is given and then drops
648  // those that are too far off that line.
649  unsigned int WindowTrackingAlg::FitWindow(std::vector<art::Ptr<rb::CellHit> > const& prevWindow,
650  std::vector<art::Ptr<rb::CellHit> > & currWindow,
651  std::vector<art::Ptr<rb::CellHit> > & removedHits,
652  std::map<unsigned short int, std::pair<float, float>,
653  bool(*)(unsigned short int, unsigned short int) > const& planeCoG,
654  geo::View_t const& view)
655  {
656  size_t prevWinSize = prevWindow.size();
657 
658  double start[2] = {0.};
659  double end[2] = {0.};
660 
661  std::vector<double> x(prevWinSize + currWindow.size(), 0.);
662  std::vector<double> y(prevWinSize + currWindow.size(), 0.);
663  std::vector<double> w(prevWinSize + currWindow.size(), 0.);
664  std::vector<size_t> planeBoundaries;
665  double chiSqr = 0.;
666 
667  // fill the vectors, first with the previous window, then the current one
668  double xyz[3] = {0.};
669  double &v = (view == geo::kX) ? xyz[0] : xyz[1];
670  for(size_t c = 0; c < prevWinSize; ++c){
671  fGeo->Plane(prevWindow[c]->Plane())->Cell(prevWindow[c]->Cell())->GetCenter(xyz);
672  x[c] = xyz[2];
673  y[c] = v;
674  w[c] = 2.;//prevWindow[c]->ADC();
675  //LOG_VERBATIM("WindowTrackingAlg")
676  //<< "previous window points: "
677  //<< prevWindow[c]->Plane() << " "
678  //<< prevWindow[c]->Cell() << " "
679  //<< x[c] << " " << y[c] << " "
680  //<< w[c];
681 
682  }
683 
684  // fit the previous window to be able to figure out which cells
685  // in the current window should go in the initial fit
686  if(prevWinSize > 0)
687  chiSqr = geo::LinFitMinDperp(x, y, w,
688  &start[0], &start[1],
689  &end[0], &end[1]);
690 
691  unsigned int ndrop = 0;
692 
693  // variables to help identify outliers early on
694  unsigned short int currPlane = currWindow.front()->Plane();
695  unsigned short int prevPlane = currWindow.front()->Plane();
696  float cog = 0;
697  for(size_t c = 0; c < currWindow.size(); ++c){
698  currPlane = currWindow[c]->Plane();
699  fGeo->Plane(currPlane)->Cell(currWindow[c]->Cell())->GetCenter(xyz);
700  cog = planeCoG.find(currPlane)->second.first;
701  cog /= planeCoG.find(currPlane)->second.second;
702  x[c + prevWinSize] = xyz[2];
703  y[c + prevWinSize] = v;
704  w[c + prevWinSize] = 1.;//currWindow[c]->ADC();
705 
706  if(currPlane != prevPlane) planeBoundaries.push_back(c);
707 
708  prevPlane = currPlane;
709 
710  // First pass, drop hits that are clear outliers (more than 2x tolerance
711  // from anything else) in the current window. If it turns out they really
712  // are good hits (eg surrounded by bad channels) they will be added back
713  // in by the last step. But this helps prevent distant noise hits from
714  // throwing off the original fit.
715  if(prevWinSize > 0){
716  const double dSq = geo::DsqrToLine(x[c + prevWinSize],
717  y[c + prevWinSize],
718  start[0], start[1],
719  end[0], end[1]);
720 
721  if(dSq > util::sqr(fDHitGood*2)){
722  //LOG_VERBATIM("WindowTrackingAlg")
723  //<< "dropping hit at (z,v) = ("
724  //<< x[c + prevWinSize] << ","
725  //<< y[c + prevWinSize]
726  //<< ") for being too separated from others: "
727  //<< dSq
728  //<< " line start: " << start[0] << " " << start[1]
729  //<< " line end: " << end[0] << " " << end[1]
730  //<< " expected pos "
731  //<< x[c + prevWinSize]*(end[1]-start[1])/(end[0]-start[0]) + (start[1] - start[0]*(end[1]-start[1])/(end[0]-start[0]));
732 
733  w[c + prevWinSize] = 0;
734  ++ndrop;
735  }
736  }// end if we have a previous window to check this hit against
737  else{
738  // set the weight to 0 of cells from the same plane
739  if(std::abs(v - cog) > fDHitGood*4){
740  w[c + prevWinSize] = 0.;
741  ++ndrop;
742  }
743  }
744 
745  //LOG_VERBATIM("WindowTrackingAlg")
746  //<< "current window points: "
747  //<< currWindow[c]->Plane() << " "
748  //<< currWindow[c]->Cell() << " "
749  //<< x[c + prevWinSize] << " "
750  //<< y[c + prevWinSize] << " "
751  //<< w[c + prevWinSize] << " "
752  //<< cog << " "
753  //<< currWindow[c]->ADC();
754 
755  }
756 
757  // check if we have dropped all the hits from the current window
758  // and have none in the previous window. In that case, we are
759  // on the first window of the given slice but have no way to do
760  // a fit. Remove all the hits from the current window
761  if(ndrop == currWindow.size() && prevWinSize < 1){
762  currWindow.clear();
763  return ndrop;
764  }
765 
766  // sometimes you get a plane in the initial window that has
767  // its hits well away from the others and that throws off the
768  // fit. Run through the initial window and drop each plane's
769  // cells in turn. If a plane causes the chi^2 to spike drop it
770  if(prevWinSize < 1){
771 
772  // add the final boundary to the list
773  planeBoundaries.push_back(w.size());
774 
775  size_t prevPB = 0;
776  double bestChi2 = 1.e6;
777  double worstChi2 = 0.;
778  std::pair<size_t, size_t> bestBounds(0, 0);
779  for(auto const p : planeBoundaries){
780 
781  // set the weights for the plane to 0
782  std::vector< std::pair<size_t, double> > prevWeights;
783  for(auto pb = prevPB; pb < p; ++pb){
784  prevWeights.push_back(std::pair<size_t, double>(pb, w[pb]));
785  w[pb] = 0.;
786  }
787 
788  // for(size_t c = prevWinSize; c < x.size(); ++c)
789  // LOG_DEBUG("WindowTrackingAlg") << "look for outliers: " << p << " "
790  // << x[c + prevWinSize] << " "
791  // << y[c + prevWinSize] << " "
792  // << w[c + prevWinSize];
793 
794  try{
795  chiSqr = geo::LinFitMinDperp(x, y, w,
796  &start[0], &start[1],
797  &end[0], &end[1]);
798  }
799  catch(cet::exception &e){
800 
801  //LOG_DEBUG("WindowTrackingAlg") << "caught exception " << e
802  // << " and moving on to next plane";
803 
804  // there were likely no non-zero weights, so just move on to the next plane
805  continue;
806  }
807 
808  //LOG_DEBUG("WindowTrackingAlg") << "chisq = " << chiSqr;
809 
810  if(chiSqr < bestChi2){
811  bestChi2 = chiSqr;
812  bestBounds.first = prevPB;
813  bestBounds.second = p;
814  //LOG_DEBUG("WindowTrackingAlg") << "best chisqr, bounds: " << bestBounds.first << " " << bestBounds.second;
815  }
816  if(chiSqr > worstChi2) worstChi2 = chiSqr;
817 
818  // set the weights back to 1 to evaluate the remaining planes
819  for(auto pw : prevWeights) w[pw.first] = pw.second;
820 
821  // set the previous boundary to the current value
822  prevPB = p;
823 
824  }// end loop over plane boundaries
825 
826  // set the worst plane contributions to 0 weight for now
827  // but only do that if it adds 2 units of chi^2 or more
828  // for each cell in the window
829  if(worstChi2 - bestChi2 > 2.*w.size())
830  for(auto b = bestBounds.first; b < bestBounds.second; ++b) w[b] = 0.;
831 
832  }// end if on the first window for a slice in the chosen view
833 
834  bool droppedAny = true;
835  while(droppedAny){
836 
837  try{
838  chiSqr = geo::LinFitMinDperp(x, y, w,
839  &start[0], &start[1],
840  &end[0], &end[1]);
841  }
842  catch(cet::exception &e){
843 
844  // there were likely no non-zero weights, so stop trying to do a fit
845  break;
846  }
847 
848  //LOG_DEBUG("WindowTracking") << "end points from fit "
849  // << start[0] << " "
850  // << start[1] << " "
851  // << end[0] << " "
852  // << end[1] << " chi^2 = "
853  // << chiSqr;
854 
855  // loop over the hits and for the ones that are more than fDHitGood cm
856  // off the line set the weight to be 0 so that they do not contribute to
857  // future fits. Only drop the worst hit, and then iterate to get the next
858  // worst and so on. Otherwise we can end up dropping lots of hits we
859  // shouldn't due to one outlier.
860  double dSqMax = 0;
861  int dSqMaxIdx = -1;
862  for(size_t n = prevWinSize; n < x.size(); ++n){
863  // Already dropped
864  if(w[n] == 0) continue;
865 
866  const double dSq = geo::DsqrToLine(x[n], y[n],
867  start[0], start[1],
868  end[0], end[1]);
869  if(dSq > dSqMax){
870  dSqMax = dSq;
871  dSqMaxIdx = n;
872  //LOG_DEBUG("WindowTracking") << "highest dSq is at: "
873  // << x[n] << " " << y[n]
874  // << " " << dSq;
875  }
876  } // end for n
877 
878  droppedAny = false;
879  if(dSqMax > util::sqr(fDHitGood)){
880  //LOG_VERBATIM("WindowTracking")
881  //<< "dropping hit at ("
882  //<< x[dSqMaxIdx] << ", "
883  //<< y[dSqMaxIdx] << ")"
884  //<< " dSq = " << dSqMax << " expected "
885  //<< x[dSqMaxIdx]*(end[1]-start[1])/(end[0]-start[0]) + (start[1] - start[0]*(end[1]-start[1])/(end[0]-start[0]));
886 
887  w[dSqMaxIdx] = 0.;
888  droppedAny = true;
889  ++ndrop;
890  }
891  }// end while
892 
893  for(size_t n = prevWinSize; n < x.size(); ++n){
894  // Already part of the fit
895  if(w[n] != 0) continue;
896 
897  const double dSq = geo::DsqrToLine(x[n], y[n],
898  start[0], start[1],
899  end[0], end[1]);
900  if(dSq < util::sqr(fDHitGood)){
901  //LOG_DEBUG("WindowTracking")
902  //<< "adding back hit at ("
903  //<< x[n] << ", " << y[n]
904  //<< ")" << " dSq = " << dSq;
905 
906  w[n] = 1.;
907  --ndrop;
908  }
909  } // end for n
910 
911  // remove any hits from the current window that have a weight of 0
912  // erase all entries from w that are from the previous window
913  auto witr = w.begin();
914  if(prevWinSize > 0)
915  witr = w.erase(w.begin(), w.begin()+prevWinSize);
916  auto citr = currWindow.begin();
917  while( witr != w.end() && citr != currWindow.end()){
918  if( (*witr) < 1. ){
919  //LOG_VERBATIM("WindowTrackingAlg")
920  //<< "dropping hit at "
921  //<< (*citr)->Plane() << ","
922  //<< (*citr)->Cell();
923  // add the hit to the removed hits vector
924  removedHits.push_back(*citr);
925 
926  citr = currWindow.erase(citr);
927  witr = w.erase(witr);
928  }
929  else{
930  ++witr;
931  ++citr;
932  }
933  }
934 
935  //LOG_DEBUG("WindowTracking") << "end points from fit "
936  // << start[0] << " " << start[1] << " "
937  // << end[0] << " " << end[1];
938 
939  return ndrop;
940  }
941 
942  //......................................................................
944  std::vector<art::Ptr<rb::CellHit> > & removedHits)
945  {
946  // map the z position of each trajectory point to the position in the view
947  std::map<double, double> zToViewPos;
948 
949  for(auto tp : track.Trajectory()){
950  if (track.View() == geo::kX) zToViewPos[tp.Z()] = tp.X();
951  else if(track.View() == geo::kY) zToViewPos[tp.Z()] = tp.Y();
952  }// end loop to map z to trajectory in view
953 
954  // loop over the rejected cells to see if any of them are close enough
955  // to the other cells in the plane to be added back to the track
956  double xyz[3] = {0.};
957  double &v = (track.View() == geo::kX) ? xyz[0] : xyz[1];
958  auto rhitr = removedHits.begin();
959  bool useHit = false;
960  while( rhitr != removedHits.end() ){
961  useHit = false;
962  auto c = *rhitr;
963  fGeo->Plane(c->Plane())->Cell(c->Cell())->GetCenter(xyz);
964 
965  if(zToViewPos.count(xyz[2]) > 0){
966 
967  // if the cell is close enough to the trajectory at this location
968  // add it to the track and remove it from the
969  if(std::abs(zToViewPos[xyz[2]] - v) < 8.*fDHitGood){
970  track.Add(c);
971  rhitr = removedHits.erase(rhitr);
972  useHit = true;
973  }
974  }// end if the z position is in the map
975 
976  // the erase method returns an iterator to the position
977  // after the removed item, ie effectively incrementing the iterator
978  // if we didn't remove the hit, we need to increment the iterator
979  if(!useHit) ++rhitr;
980 
981  }// end loop over previously rejected hits
982 
983  return;
984  }
985 
986 
987  //......................................................................
989  {
990 
991  // use at most the first N (default = 20) trajectory points on a track to get
992  // the incoming direction
993  size_t numTrajPts = (track.NTrajectoryPoints() > fNTrajectory) ? fNTrajectory : track.NTrajectoryPoints();
994 
995  // put the first n trajectory points into vectors
996  std::vector<double> x(numTrajPts, 0.);
997  std::vector<double> y(numTrajPts, 0.);
998  std::vector<double> z(numTrajPts, 0.);
999  std::vector<double> w(numTrajPts, 1.);
1000 
1001  TVector3 pos;
1002  for(size_t t = 0; t < numTrajPts; ++t){
1003  pos = track.TrajectoryPoint(t);
1004  x[t] = pos.X();
1005  y[t] = pos.Y();
1006  z[t] = pos.Z();
1007  }
1008 
1009  double xzStart[2] = {0.};
1010  double xzEnd[2] = {0.};
1011  double yzStart[2] = {0.};
1012  double yzEnd[2] = {0.};
1013 
1014  geo::LinFitMinDperp(z, x, w, &xzStart[0], &xzStart[1], &xzEnd[0], &xzEnd[1]);
1015  geo::LinFitMinDperp(z, y, w, &yzStart[0], &yzStart[1], &yzEnd[0], &yzEnd[1]);
1016 
1017  if(xzStart[0] < z.front()) std::swap(xzStart, xzEnd);
1018  if(yzStart[0] < z.front()) std::swap(yzStart, yzEnd);
1019 
1020  TVector3 start(xzStart[1], yzStart[1], z.front());
1021  TVector3 end (xzEnd[1], yzEnd[1], z.back());
1022 
1023  track.SetDir(end - start);
1024 
1025  return;
1026  }
1027 
1028  //......................................................................
1029  // method to make 3D tracks from views with only 1 or 2 hit planes in the slice
1030  rb::Track WindowTrackingAlg::ShortTrack(std::map<unsigned short int, std::vector<art::Ptr<rb::CellHit> >,
1031  bool(*)(unsigned short int, unsigned short int) > & xViewHits,
1032  std::map<unsigned short int, std::vector<art::Ptr<rb::CellHit> >,
1033  bool(*)(unsigned short int, unsigned short int) > & yViewHits,
1034  std::map<unsigned short int, std::pair<float, float>,
1035  bool(*)(unsigned short int, unsigned short int) > const& planeCoG,
1036  double & chiSqr)
1037  {
1038 
1039  // the assumption for these tracks is that they are very steep and leave lots of
1040  // hits in only at most 3 planes in a view. Just use the two planes with the
1041  // most hits to set the track direction and then add in the hits from the third
1042  // plane if they fall along that line
1043  // make a copy of the input maps and then drop the planes with
1044  // the fewest hits in each view. At this point there should
1045  // be only two planes in each view
1046  auto cogCopy(planeCoG);
1047  auto xHitsCopy(xViewHits);
1048  auto yHitsCopy(yViewHits);
1049  unsigned short int removedXPlane = 0;
1050  unsigned short int removedYPlane = 0;
1051  if(xViewHits.size() > 2) this->ShortTrackRemoveExtraPlane(xHitsCopy, cogCopy, removedXPlane);
1052  if(yViewHits.size() > 2) this->ShortTrackRemoveExtraPlane(yHitsCopy, cogCopy, removedYPlane);
1053 
1054  // make sure the two plane hit view has the planes separated by only 1 plane in the other view, and that
1055  // the one hit plane view is between them
1056  int firstXPlane = xHitsCopy.cbegin()->first;
1057  int secondXPlane = 0;
1058  int firstYPlane = yHitsCopy.cbegin()->first;
1059  int secondYPlane = 0;
1060  float firstXCoG = cogCopy.find(firstXPlane)->second.first/cogCopy.find(firstXPlane)->second.second;
1061  float firstYCoG = cogCopy.find(firstYPlane)->second.first/cogCopy.find(firstYPlane)->second.second;
1062  float secondXCoG = 0.;
1063  float secondYCoG = 0.;
1064  double xyz[3] = {0.};
1065  fGeo->Plane(xHitsCopy.cbegin()->first)->Cell((xHitsCopy.cbegin()->second.front())->Cell())->GetCenter(xyz);
1066  float firstXZ = xyz[2];
1067  fGeo->Plane(yHitsCopy.cbegin()->first)->Cell((yHitsCopy.cbegin()->second.front())->Cell())->GetCenter(xyz);
1068  float firstYZ = xyz[2];
1069  float secondXZ = xyz[2];
1070  float secondYZ = xyz[2];
1071  bool twoXPlanes = (xHitsCopy.size() == 2) ? true : false;
1072  bool twoYPlanes = (yHitsCopy.size() == 2) ? true : false;
1073 
1074  rb::Track retTrack;
1075 
1076 
1077  if(twoXPlanes){
1078  secondXPlane = xHitsCopy.crbegin()->first;
1079  secondXCoG = cogCopy.find(secondXPlane)->second.first/cogCopy.find(secondXPlane)->second.second;
1080  fGeo->Plane(xHitsCopy.crbegin()->first)->Cell((xHitsCopy.crbegin()->second.front())->Cell())->GetCenter(xyz);
1081  secondXZ = xyz[2];
1082  }
1083  if(twoYPlanes){
1084  secondYPlane = yHitsCopy.crbegin()->first;
1085  secondYCoG = cogCopy.find(secondYPlane)->second.first/cogCopy.find(secondYPlane)->second.second;
1086  fGeo->Plane(yHitsCopy.crbegin()->first)->Cell((yHitsCopy.crbegin()->second.front())->Cell())->GetCenter(xyz);
1087  secondYZ = xyz[2];
1088  }
1089 // LOG_VERBATIM("WindowTrackingAlg") << firstXZ << " " << firstXCoG << " " << secondXZ << " " << secondXCoG;
1090 // LOG_VERBATIM("WindowTrackingAlg") << firstYZ << " " << firstYCoG << " " << secondYZ << " " << secondYCoG;
1091 
1092 
1093  bool makeTrack = false;
1094 
1095  // check that the planes in each view are contiguous
1096  if(twoXPlanes && !twoYPlanes){
1097  if(std::abs(firstXPlane - secondXPlane) < 3 &&
1098  std::abs(firstXPlane - firstYPlane) < 2 &&
1099  std::abs(secondXPlane - firstYPlane) < 2)
1100  makeTrack = true;
1101  }
1102  else if(!twoXPlanes && twoYPlanes){
1103  if(std::abs(firstYPlane - secondYPlane) < 3 &&
1104  std::abs(firstYPlane - firstXPlane) < 2 &&
1105  std::abs(secondYPlane - firstXPlane) < 2)
1106  makeTrack = true;
1107  }
1108  else if(twoXPlanes && twoYPlanes){
1109  if(std::abs(firstYPlane - secondYPlane) < 3 &&
1110  std::abs(firstXPlane - secondXPlane) < 3 &&
1111  std::abs(firstYPlane - firstXPlane ) < 2 &&
1112  std::abs(secondYPlane - secondXPlane) < 2)
1113  makeTrack = true;
1114  }
1115 
1116  if(makeTrack){
1117 
1118 // LOG_VERBATIM("WindowTrackingAlg") << "making the short tracks in each view";
1119 
1120  rb::Track xTrack = this->ShortViewTrack(xHitsCopy, cogCopy, twoXPlanes, firstXPlane, secondXPlane, geo::kX);
1121  rb::Track yTrack = this->ShortViewTrack(yHitsCopy, cogCopy, twoYPlanes, firstYPlane, secondYPlane, geo::kY);
1122 
1123  // for a short track all the points are used, make the chiSqr 1
1124  chiSqr = 1.;
1125 
1126  // don't use ZipWith to make the 3D track - it is simple enough that
1127  // we can do it here
1128  // get the direction cosines in dx/dz and dy/dz
1129  TVector3 dxdz = xTrack.Dir();
1130  TVector3 dydz = yTrack.Dir();
1131 
1132  // figure out the trajectory points for each plane
1133  std::vector<std::pair<geo::View_t, TVector3> > trajPoints(4);
1134  trajPoints[0] = std::pair<geo::View_t, TVector3>(geo::kX, TVector3(firstXCoG, 0., firstXZ ));
1135  trajPoints[1] = std::pair<geo::View_t, TVector3>(geo::kY, TVector3(0., firstYCoG, firstYZ ));
1136  trajPoints[2] = std::pair<geo::View_t, TVector3>(geo::kX, TVector3(secondXCoG, 0., secondXZ));
1137  trajPoints[3] = std::pair<geo::View_t, TVector3>(geo::kY, TVector3(0., secondYCoG, secondYZ));
1138 
1139  // remove the second point in a view if it doesn't exist
1140  if(!twoXPlanes) trajPoints.erase(trajPoints.begin()+2);
1141  if(!twoYPlanes) trajPoints.erase(trajPoints.begin()+trajPoints.size()-1);
1142 
1143  // check if we need to add in an extra plane in either view
1144  if(xViewHits.size() > 2) this->ShortTrackExtraPlane(xViewHits, planeCoG, trajPoints, xTrack, removedXPlane);
1145  if(yViewHits.size() > 2) this->ShortTrackExtraPlane(yViewHits, planeCoG, trajPoints, yTrack, removedYPlane);
1146 
1147  // now sort the list according to the z value of the planes
1148  if(cogCopy.cbegin()->first < cogCopy.crbegin()->first)
1149  std::sort(trajPoints.begin(), trajPoints.end(), ltvec);
1150  else
1151  std::sort(trajPoints.begin(), trajPoints.end(), gtvec);
1152 
1153  double deltaZ = 0.;
1154  TVector3 neighbor;
1156  // loop over the trajectory points and set the x & y positions
1157  // be sure to get the closest neighboring trajectory point in the
1158  // opposite view to be able to set the 3D positions correctly
1159  // and use the direction cosines from the 2D tracks
1160  for(size_t tp = 0; tp < trajPoints.size(); ++tp){
1161  view = trajPoints[tp].first;
1162  if(tp < 1){
1163  deltaZ = trajPoints[tp].second.Z() - trajPoints[tp+1].second.Z();
1164  neighbor = trajPoints[tp+1].second;
1165  }
1166  else{
1167  deltaZ = trajPoints[tp].second.Z() - trajPoints[tp-1].second.Z();
1168  neighbor = trajPoints[tp-1].second;
1169  }
1170  if (view == geo::kX) trajPoints[tp].second.SetY(neighbor.Y() + deltaZ*dydz.Y()/dydz.Z());
1171  else if(view == geo::kY) trajPoints[tp].second.SetX(neighbor.X() + deltaZ*dxdz.X()/dxdz.Z());
1172 
1173  //trajPoints[tp].second.Print();
1174  } // end loop over trajectory points
1175 
1176  // make the 3D track
1177  rb::Track retTrack(xTrack.AllCells(),
1178  trajPoints.front().second,
1179  TVector3(trajPoints.back().second.X()-trajPoints.front().second.X(),
1180  trajPoints.back().second.Y()-trajPoints.front().second.Y(),
1181  trajPoints.back().second.Z()-trajPoints.front().second.Z())
1182  );
1183  retTrack.Add(yTrack.AllCells());
1184  for(size_t tp = 1; tp < trajPoints.size(); ++tp) retTrack.AppendTrajectoryPoint(trajPoints[tp].second);
1185 
1186 
1187  this->CheckTrackDirectionInY(retTrack, trajPoints.front().second, trajPoints.back().second);
1188 
1189  // restrict track to exist inside detector; if this isn't possible,
1190  // the track is no good, so don't write
1191  if(geo::ClampRayToDetectorVolume(&trajPoints.front().second,
1192  &trajPoints.back().second,
1193  &(*fGeo))) return retTrack;
1194  //else LOG_DEBUG("WindowTrackAlg") << "Cannot clamp short track to detector";
1195 
1196  }// end if this works
1197 
1198  return retTrack;
1199 
1200  }
1201 
1202  //----------------------------------------------------------------------
1203  // method to determine which plane in a set of 3 (for a short track) has the least signal
1204  // and should be removed to determine the track
1205  void WindowTrackingAlg::ShortTrackRemoveExtraPlane(std::map<unsigned short int, std::vector<art::Ptr<rb::CellHit> >,
1206  bool(*)(unsigned short int, unsigned short int) > & viewHits,
1207  std::map<unsigned short int, std::pair<float, float>,
1208  bool(*)(unsigned short int, unsigned short int) > & planeCoG,
1209  unsigned short int & removedPlane)
1210  {
1211  unsigned short int minPlane = 0;
1212 
1213  size_t minHits = 1e6;
1214  for(auto vh : viewHits){
1215  if(vh.second.size() < minHits){
1216  minPlane = vh.first;
1217  minHits = vh.second.size();
1218  }
1219  }
1220 
1221  for(auto cogItr : planeCoG){
1222  if(cogItr.first == minPlane){
1223  removedPlane = cogItr.first;
1224  viewHits.erase(removedPlane);
1225  planeCoG.erase(removedPlane);
1226  }
1227  }
1228 
1229  return;
1230  }
1231 
1232 
1233  //----------------------------------------------------------------------
1234  rb::Track WindowTrackingAlg::ShortViewTrack(std::map<unsigned short int, std::vector<art::Ptr<rb::CellHit> >,
1235  bool(*)(unsigned short int, unsigned short int) > & viewHits,
1236  std::map<unsigned short int, std::pair<float, float>,
1237  bool(*)(unsigned short int, unsigned short int) > const& planeCoG,
1238  bool const& twoPlanes,
1239  int const& firstPlane,
1240  int const& secondPlane,
1241  geo::View_t const& view)
1242 
1243  {
1244  // make 2d tracks using the center of gravity for each plane
1245  std::vector<art::Ptr<rb::CellHit> > hits;
1246  for(auto itr : viewHits){
1247  for(auto ch : itr.second) hits.push_back(ch);
1248  }
1249 
1250  double xyz1[3] = {0.};
1251  double xyz2[3] = {0.};
1252  auto chVecItr = viewHits.cbegin()->second.cbegin();
1253  fGeo->Plane(firstPlane) ->Cell((*chVecItr)->Cell())->GetCenter(xyz1);
1254  float cog1 = planeCoG.find(firstPlane) ->second.first/planeCoG.find(firstPlane) ->second.second;
1255  float deltaW = 0.;
1256  float deltaZ = 2.*fCellHalfDepth;
1257  if(twoPlanes){
1258  chVecItr = viewHits.crbegin()->second.cbegin();
1259  fGeo->Plane(secondPlane)->Cell((*chVecItr)->Cell())->GetCenter(xyz2);
1260  deltaW = (planeCoG.find(secondPlane)->second.first/planeCoG.find(secondPlane)->second.second);
1261  deltaZ = xyz2[2]-xyz1[2];
1262  }
1263  else{
1264  // hits only on one plane, so use the center of the first and last cells in the plane to find the
1265  // spread in the current view
1266  chVecItr = viewHits.cbegin()->second.cbegin();
1267  fGeo->Plane(firstPlane)->Cell((*chVecItr)->Cell())->GetCenter(xyz1);
1268  cog1 = (view == geo::kX) ? xyz1[0] : xyz1[1];
1269  fGeo->Plane(firstPlane)->Cell((viewHits.cbegin()->second.back())->Cell())->GetCenter(xyz2);
1270  deltaW = (view == geo::kX) ? xyz2[0] : xyz2[1];
1271  }
1272  deltaW -= cog1;
1273 
1274  // be sure we don't just have one cell in the view
1275  if(deltaW == 0.) deltaW = fCellHalfDepth;
1276 
1277 // LOG_VERBATIM("WindowTrackingAlg") << view << " " << cog1 << " " << xyz1[2] << " " << deltaW << " " << deltaZ;
1278 
1279  rb::Track track(hits, view, cog1, xyz1[2], deltaW, deltaZ);
1280  track.AppendTrajectoryPoint(deltaW+cog1, deltaZ);
1281 
1282  return track;
1283  }
1284 
1285  //......................................................................
1286  // method to add back in an extra plane to a short track - short tracks
1287  // are initially found using only the two planes with the most signal, it is
1288  // possible there is another plane to add back in
1289  void WindowTrackingAlg::ShortTrackExtraPlane(std::map<unsigned short int, std::vector<art::Ptr<rb::CellHit> >,
1290  bool(*)(unsigned short int, unsigned short int) > const& viewHits,
1291  std::map<unsigned short int, std::pair<float, float>,
1292  bool(*)(unsigned short int, unsigned short int) > const& planeCoG,
1293  std::vector<std::pair<geo::View_t, TVector3> > & trajPoints,
1294  rb::Track & viewTrack,
1295  unsigned short int const& removedPlane)
1296  {
1297 
1298  // get the expected trajectory point for that view
1299  double xyz[3] = {0.};
1300  fGeo->Plane(removedPlane)->Cell(0)->GetCenter(xyz);
1301  TVector3 trajPoint;
1302  float dirZ = viewTrack.Dir().Z();
1303  float dirV = 0.;
1304  float cog = planeCoG.find(removedPlane)->second.first/planeCoG.find(removedPlane)->second.second;
1305  float expectedV = 0.;
1306  float dV = 0.;
1307  if(viewTrack.View() == geo::kX){
1308  trajPoint = trajPoints[2].second;
1309  dirV = viewTrack.Dir().X();
1310  expectedV = trajPoint.X() + (dirV/dirZ)*std::abs(xyz[2] - trajPoint.Z());
1311  dV = cog - trajPoint.X();
1312 
1313  // set the trajectory point for this plane in case we want to keep it
1314  trajPoint.SetXYZ(cog, 0., xyz[2]);
1315  }
1316  else{
1317  trajPoint = trajPoints[3].second;
1318  dirV = viewTrack.Dir().Y();
1319  expectedV = trajPoint.Y() + (dirV/dirZ)*std::abs(xyz[2] - trajPoint.Z());
1320  dV = cog - trajPoint.Y();
1321 
1322  // set the trajectory point for this plane in case we want to keep it
1323  trajPoint.SetXYZ(0., cog, xyz[2]);
1324  }
1325 
1326  // compare the expected trajectory location to the center of gravity for that plane, and if they are
1327  // within 8*fDHitGood, ie a very steep track, add it back in
1328  if(std::abs(expectedV - cog) < 8*fDHitGood && dV*dirV > 0){
1329  trajPoints.push_back(std::pair<geo::View_t, TVector3>(viewTrack.View(), trajPoint));
1330  for(auto h : viewHits.at(removedPlane)) viewTrack.Add(h);
1331  }
1332 
1333  return;
1334  }
1335 
1336 
1337  //......................................................................
1338  // check that the track is going the correct direction for a
1339  // cosmic ray, ie starts at high y and ends at low y. The end
1340  // points have to vary by at least fDHitGood to trigger a switch
1342  TVector3 & start,
1343  TVector3 & end)
1344  {
1345 
1346  if(start.Y() < end.Y()){
1347 
1348  TVector3 pt2 = track.TrajectoryPoint(track.NTrajectoryPoints()-2);
1349  rb::Track ctrack(track.AllCells(), end, pt2-end);
1350  for(int p = track.NTrajectoryPoints()-2; p >= 0; --p)
1351  ctrack.AppendTrajectoryPoint(track.TrajectoryPoint(p));
1352 
1353  track = ctrack;
1354 
1355  start = ctrack.Start();
1356  end = ctrack.Stop();
1357  }// end if we have to flip the direction on this track
1358 
1359  return;
1360  }
1361 
1362  //......................................................................
1363  // this method is for a single view, XZ or YZ
1364  bool WindowTrackingAlg::ValidTrack(std::map<unsigned short int, std::vector<art::Ptr<rb::CellHit> >,
1365  bool(*)(unsigned short int, unsigned short int) > const& planeToHit)
1366  {
1367  if(planeToHit.size() < fMinViewPlanes) return false;
1368 
1369  unsigned int const& startPlane = planeToHit.cbegin()->first;
1370  unsigned int const& endPlane = planeToHit.crbegin()->first;
1371 
1372  // the 0.5 in the MinHitFraction accounts for the fact that
1373  // every other plane is X or Y as we are using the full extent
1374  // of the slice
1375  if(planeToHit.size() < fMinHitFraction*0.5*std::abs(1.*startPlane - 1.*endPlane) )
1376  return false;
1377 
1378  return true;
1379  }
1380 
1381  //----------------------------------------------------------------------------
1383  {
1384  // look for any trajectory points that have a Nan
1385 // std::cout << " checking out the track before endpoints set" << std::endl;
1386 // for(auto tp : track.Trajectory())
1387 // if(std::isnan(tp.X()) ||
1388 // std::isnan(tp.Y()) ||
1389 // std::isnan(tp.Z()) ) tp.Print();
1390 
1391  // extrapolate the track along the direction cosines
1392  // first get the cells from the first and last planes
1393  std::vector<art::Ptr<rb::CellHit> > startPlaneCells;
1394  std::vector<art::Ptr<rb::CellHit> > endPlaneCells;
1395  unsigned int minPlane = track.MinPlane();
1396  unsigned int maxPlane = track.MaxPlane();
1397  for(auto ch : track.AllCells() ){
1398  if(ch->Plane() == maxPlane) endPlaneCells.push_back(ch);
1399  if(ch->Plane() == minPlane) startPlaneCells.push_back(ch);
1400  }
1401 
1402  std::sort(startPlaneCells.begin(), startPlaneCells.end(), sortCellPtr);
1403  std::sort(endPlaneCells.begin(), endPlaneCells.end(), sortCellPtr);
1404 
1405  // figure out the direction from the next to last point to the last point
1406  TVector3 endPt = track.Trajectory().back();
1407  TVector3 endDir = endPt - track.Trajectory()[track.NTrajectoryPoints()-2];
1408 
1409  // figure out the direction from the second point to the first point
1410  TVector3 startPt = track.Trajectory().front();
1411  TVector3 startDir = track.Trajectory()[1] - startPt;
1412 
1413  geo::View_t endPlaneView = fGeo->Plane(maxPlane)->View();
1414  geo::View_t startPlaneView = fGeo->Plane(minPlane)->View();
1415 
1416  // determining which is the startPlane of the track and the endPlane
1417  // depends on the z direction. Swap the min and max vectors if the
1418  // direction cosine in z is negative
1419  if(endDir.Z() < 0.){
1420  startPlaneCells.swap(endPlaneCells);
1421  std::swap(startPlaneView, endPlaneView);
1422  }
1423 
1424  // get the highest cell in the start plane and the lowest in the
1425  // end plane
1426  unsigned short int startCell = startPlaneCells.back()->Cell();
1427  unsigned short int endCell = endPlaneCells.front()->Cell();
1428  // if the start or end plane is in the X view, then we need to
1429  // determine the extrema for the cells based on the track
1430  // trajectory
1431  if(startPlaneView == geo::kX && startDir.X() > 0.) startCell = startPlaneCells.front()->Cell();
1432  if(endPlaneView == geo::kX && endDir.X() > 0.) endCell = endPlaneCells.back()->Cell();
1433 
1434  // point the track back from the starting trajectory point to
1435  // a better location accounting for the potential difference between
1436  // the center of gravity in a plane and its final cell in the plane
1437  double xyz[3] = {0.};
1438  fGeo->Plane(startPlaneCells.front()->Plane())->Cell(startCell)->GetCenter(xyz);
1439  this->FindEndPoint(startPlaneView, xyz, track.Start(), track.Dir(), true, track);
1440 
1441  // carry the track further in its direction of travel
1442  fGeo->Plane(endPlaneCells.front()->Plane())->Cell(endCell)->GetCenter(xyz);
1443  this->FindEndPoint(endPlaneView, xyz, track.Stop(), track.StopDir(), false, track);
1444 
1445  return;
1446  }
1447 
1448  //----------------------------------------------------------------------------
1450  double * cellCenter,
1451  TVector3 const& endPoint,
1452  TVector3 const& epDir,
1453  bool const& start,
1454  rb::Track & track)
1455  {
1456  double deltaZ = fGeo->Plane(0)->Cell(0)->HalfD();
1457 
1458  double &v = cellCenter[0];
1459  double epV = endPoint.X();
1460  double epVDir = epDir.X();
1461  int signDirV = 1.;
1462 
1463 // std::string startEnd("start");
1464 // if(!start) startEnd = "end";
1465 
1466  if(view == geo::kY){
1467  v = cellCenter[1];
1468  epV = endPoint.Y();
1469  epVDir = epDir.Y();
1470  }
1471 
1472 // LOG_VERBATIM("WindowTrackingAlg") << startEnd << " point: "
1473 // << endPoint.X() << " " << endPoint.Y() << " " << endPoint.Z()
1474 // << " " << epDir.X() << " " << epDir.Y() << " " << epDir.Z()
1475 // << " cell at "
1476 // << v << " " << cellCenter[2] << " " << deltaZ;
1477 
1478  // bail if the direction in the view is 0
1479  // the current end point will have to be good enough
1480  if(epVDir == 0.) return;
1481 
1482  double deltaV = std::abs(v - epV);
1483  TVector3 dV(epDir.X()/epVDir, epDir.Y()/epVDir, epDir.Z()/epVDir);
1484  // point the track backwards
1485  if (epVDir > 0. && start) signDirV = -1;
1486  else if(epVDir < 0. && !start) signDirV = -1;
1487  TVector3 pt(endPoint.X() + deltaV*signDirV*dV.X(),
1488  endPoint.Y() + deltaV*signDirV*dV.Y(),
1489  endPoint.Z() + deltaV*signDirV*dV.Z());
1490 
1491  // only prepend if the point is different from the original starting point
1492  if(pt != endPoint){
1493 
1494  // if the new start point is more than half a cell width
1495  // from the z position of the first cell in the start plane
1496  // just carry the track to the boundary of the plane
1497  if(std::abs(pt.Z() - cellCenter[2]) > deltaZ){
1498  signDirV = 1.;
1499  if (epDir.Z() > 0. && start) signDirV = -1.;
1500  else if(epDir.Z() < 0. && !start) signDirV = -1.;
1501  pt.SetXYZ(endPoint.X() + deltaZ*signDirV*epDir.X()/epDir.Z(),
1502  endPoint.Y() + deltaZ*signDirV*epDir.Y()/epDir.Z(),
1503  endPoint.Z() + deltaZ*signDirV);
1504  }
1505 
1506  // check that we don't have nans for the new coordinates
1507  if( !std::isnan(pt.X()) &&
1508  !std::isnan(pt.Y()) &&
1509  !std::isnan(pt.Z())){
1510  if(start) track.PrependTrajectoryPoint(pt);
1511  else track.AppendTrajectoryPoint(pt);
1512 
1513  // LOG_VERBATIM("WindowTrackingAlg") << "new " << startEnd << " point is: "
1514  // << pt.X() << " "
1515  // << pt.Y() << " "
1516  // << pt.Z();
1517 
1518  }
1519 
1520  } // end if the new point and the previous endPoint are different
1521 
1522  return;
1523  }
1524 
1525 
1526 }
1527 
1528 
1529 
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
An algorithm to perform cosmic ray track fitting.
Definition: TrackAlg.h:25
double HalfD() const
Definition: CellGeo.cxx:205
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
TVector3 TrajectoryPoint(unsigned int i) const
The ith point on the trajectory, a 3-vector in cm.
Definition: Track.cxx:158
ctrack
Definition: test3.py:11
rb::Track ShortViewTrack(std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &viewHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, bool const &twoPlanes, int const &firstPlane, int const &secondPlane, geo::View_t const &view)
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
constexpr T pow(T x)
Definition: pow.h:75
Vertical planes which measure X.
Definition: PlaneGeo.h:28
rb::Track MakeViewTrack(std::vector< art::Ptr< rb::CellHit > > const &trkHits, geo::View_t const &view, double &chiSqr, bool const &dir, std::vector< art::Ptr< rb::CellHit > > &removedHits)
Definition: event.h:19
A collection of associated CellHits.
Definition: Cluster.h:47
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
WindowTrackingAlg(fhicl::ParameterSet const &pset)
unsigned int fMinHitsTryAgain
Minimum number of hits left to look for another track.
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
bool ClampRayToDetectorVolume(TVector3 *p0, TVector3 *p1, const GeometryBase *geom)
If either endpoint is outside the detector move it to the edge.
Definition: Geo.cxx:640
const PlaneGeo * Plane(unsigned int i) const
virtual TVector3 Start() const
Definition: Prong.h:73
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
unsigned int FitWindow(std::vector< art::Ptr< rb::CellHit > > const &prevWindow, std::vector< art::Ptr< rb::CellHit > > &currWindow, std::vector< art::Ptr< rb::CellHit > > &removedHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, geo::View_t const &view)
rb::Track ShortTrack(std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &xViewHits, std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &yViewHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, double &chiSqr)
float abs(float number)
Definition: d0nt_math.hpp:39
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
art::PtrVector< rb::CellHit > AllCells() const
Get all cells from both views.
Definition: Cluster.cxx:180
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
virtual void Add(const art::Ptr< rb::CellHit > &cell, double weight=1)
Definition: Cluster.cxx:84
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
Float_t Z
Definition: plot.C:38
rb::Track FitView(std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &hits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, std::vector< art::Ptr< rb::CellHit > > &removedHits, geo::View_t view, double &chiSqr)
void CheckTrackDirectionInY(rb::Track &track, TVector3 &start, TVector3 &end)
art::ServiceHandle< geo::Geometry > fGeo
handle to the geometry
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
std::vector< rb::Track > MakeTrack(const rb::Cluster *slice)
rb::Track Make3DTrack(std::vector< art::Ptr< rb::CellHit > > &allHits, std::vector< art::Ptr< rb::CellHit > > &removedHits, double &chiSqr, bool(*comp)(unsigned short int, unsigned short int))
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:53
static const double pb
Definition: Units.h:90
void hits()
Definition: readHits.C:15
void prev()
Definition: show_event.C:91
double fCellHalfDepth
half of cell size in z
double dz[NP][NC]
void LookForBremsstrahlungHits(rb::Track &track, std::vector< art::Ptr< rb::CellHit > > &removedHits)
double LinFitMinDperp(const std::vector< double > &x, const std::vector< double > &y, const std::vector< double > &w, double *x1, double *y1, double *x2, double *y2)
Find the best-fit line to a collection of points in 2-D by minimizing the squared perpendicular dista...
Definition: Geo.cxx:449
bool sortCellPtr(art::Ptr< rb::CellHit > const a, art::Ptr< rb::CellHit > const b)
const double a
void SetTrackEndPoints(rb::Track &track)
double DsqrToLine(double x0, double y0, double x1, double y1, double x2, double y2)
In two dimensions, return the perpendicular distance from a point (x0,y0) to a line defined by end po...
Definition: Geo.cxx:338
virtual void SetDir(TVector3 dir)
Definition: Prong.cxx:104
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
T get(std::string const &key) const
Definition: ParameterSet.h:231
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
unsigned int fMinViewPlanes
Minimum number of planes in a view to fit.
Collect Geo headers and supply basic geometry functions.
Var weights
void FindEndPoint(geo::View_t const &view, double *cellCenter, TVector3 const &endPoint, TVector3 const &endDir, bool const &start, rb::Track &track)
void PrependTrajectoryPoint(TVector3 pt)
Support constructing tracks backwards.
Definition: Track.cxx:134
z
Definition: test.py:28
double fMinHitFraction
Minimum value of hits/plane to look for a track.
#define LOG_WARNING(category)
void ShortTrackRemoveExtraPlane(std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > &viewHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > &planeCoG, unsigned short int &removedPlane)
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
void DetermineInitialDirection(rb::Track &track)
reference front()
Definition: PtrVector.h:379
bool ltvec(std::pair< geo::View_t, TVector3 > a, std::pair< geo::View_t, TVector3 > b)
virtual TVector3 StopDir() const
Get direction at the stop point.
Definition: Track.cxx:193
TDirectory * dir
Definition: macro.C:5
bool gtvec(std::pair< geo::View_t, TVector3 > a, std::pair< geo::View_t, TVector3 > b)
unsigned int fWindowSize
Number of planes in the sliding window.
void ShortTrackExtraPlane(std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > const &viewHits, std::map< unsigned short int, std::pair< float, float >, bool(*)(unsigned short int, unsigned short int) > const &planeCoG, std::vector< std::pair< geo::View_t, TVector3 > > &trajPoints, rb::Track &viewTrack, unsigned short int const &removedPlane)
std::vector< TVector3 > const & Trajectory() const
return a constant reference to the track trajectory points
Definition: Track.h:86
const hit & b
Definition: hits.cxx:21
TRandom3 r(0)
unsigned int MaxPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:508
bool ValidTrack(std::map< unsigned short int, std::vector< art::Ptr< rb::CellHit > >, bool(*)(unsigned short int, unsigned short int) > const &planeToHit)
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
bool gt(unsigned short int a, unsigned short int b)
unsigned int fNTrajectory
number of trajectory points to use in determining incoming direction
Float_t e
Definition: plot.C:35
Track finder for cosmic rays.
Float_t track
Definition: plot.C:35
bool lt(unsigned short int a, unsigned short int b)
int fMaxPlaneSeparation
Maximum number of planes allowed to separate planes with hits in slice.
Float_t w
Definition: plot.C:20
void next()
Definition: show_event.C:84