Track.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file Track.cxx
3 // \brief A track is a \ref Prong with full reconstructed trajectory.
4 // \version $Id: Track.cxx,v 1.22 2012-11-16 19:47:05 bckhouse Exp $
5 // \author Christopher Backhouse - bckhouse@caltech.edu
6 ////////////////////////////////////////////////////////////////////////
7 #include <algorithm>
8 
9 #include "RecoBase/Track.h"
10 #include "Geometry/Geometry.h"
11 #include "GeometryObjects/Geo.h"
13 #include "RecoBase/CellHit.h"
14 
17 #include "cetlib_except/exception.h"
18 
19 namespace rb
20 {
21  //......................................................................
23  TVector3 start,
24  TVector3 dir,
25  int id)
26  : rb::Prong(cells, start, dir, id)
27  {
28  AppendTrajectoryPoint(start);
29  }
30 
31  //..................................................................
34  double v0,
35  double z0,
36  double dv,
37  double dz,
38  int id)
39  : rb::Prong(cells, view, v0, z0, dv, dz, id)
40  {
41  AppendTrajectoryPoint(v0, z0);
42  }
43 
44  //..................................................................
45  Track::Track(const std::vector< art::Ptr<rb::CellHit> >& cells,
47  double v0,
48  double z0,
49  double dv,
50  double dz,
51  int id)
52  : rb::Prong(cells, view, v0, z0, dv, dz, id)
53  {
54  AppendTrajectoryPoint(v0, z0);
55  }
56 
57  //......................................................................
58  Track::Track(const rb::Cluster& clust,
59  TVector3 start,
60  TVector3 dir,
61  int id)
62  : rb::Prong(clust, start, dir, id)
63  {
64  AppendTrajectoryPoint(start);
65  }
66 
67  //......................................................................
68  Track::Track(const rb::Cluster& clust,
69  double v0,
70  double z0,
71  double dv,
72  double dz,
73  int id)
74  : rb::Prong(clust, v0, z0, dv, dz, id)
75  {
76  AppendTrajectoryPoint(v0, z0);
77  }
78 
79  //......................................................................
80  Track::Track(const rb::Prong& prong,
81  int id)
82  : rb::Prong(prong)
83  {
84  SetID(id);
86  }
87 
88  //......................................................................
90  std::vector<double> const& weights,
91  std::vector<TVector3> const& trajectory,
92  TVector3 const& start,
93  TVector3 const& dir,
94  int const& id,
96  : Prong(cells, start, dir, id)
97  , fTraj(trajectory)
98  {
99  fView = view;
100  // set the cell weights
101  for(size_t w = 0; w < weights.size(); ++w) SetWeight(w, weights[w]);
102 
103  return;
104  }
105 
106  //......................................................................
108  {
109  }
110 
111  //......................................................................
113  {
114  fPrecalcTotalGeV = -1; // Invalidate any cached value
115 
116  fTraj.push_back(pt);
117  }
118 
119  //......................................................................
120  void Track::AppendTrajectoryPoint(double pv, double pz)
121  {
122  fPrecalcTotalGeV = -1; // Invalidate any cached value
123 
124  if(fView == geo::kXorY)
125  throw cet::exception("Track") << "view is geo::kXorY, cannot add a 2D trajectroy point";
126 
127  if(fView == geo::kX)
128  fTraj.push_back(TVector3(pv, 0, pz));
129  else
130  fTraj.push_back(TVector3(0, pv, pz));
131  }
132 
133  //......................................................................
135  {
136  fPrecalcTotalGeV = -1; // Invalidate any cached value
137 
138  fTraj.insert(fTraj.begin(), pt);
139 
140  fStart = pt;
141  }
142 
143  //......................................................................
144  void Track::PrependTrajectoryPoint(double pv, double pz)
145  {
146  fPrecalcTotalGeV = -1; // Invalidate any cached value
147 
148  if(fView == geo::kXorY)
149  throw cet::exception("Track") << "view is geo::kXorY, cannot add a 2D trajectroy point";
150 
151  if(fView == geo::kX)
152  PrependTrajectoryPoint(TVector3(pv, 0, pz));
153  else
154  PrependTrajectoryPoint(TVector3(0, pv, pz));
155  }
156 
157  //......................................................................
158  TVector3 Track::TrajectoryPoint(unsigned int i) const
159  {
160  if(i > fTraj.size())
161  throw cet::exception("Track") << "attempt to get trajectory point "
162  << i << " is out of range: " << fTraj.size();
163 
164  return fTraj[i];
165  }
166 
167  //......................................................................
168  void Track::SetStart(TVector3 start)
169  {
170  Prong::SetStart(start);
171 
172  if(fTraj.empty()) fTraj.resize(1);
173  fTraj[0] = fStart;
174  }
175 
176  //......................................................................
177  void Track::SetStart(double v0, double z0)
178  {
179  Prong::SetStart(v0, z0);
180 
181  if(fTraj.empty()) fTraj.resize(1);
182  fTraj[0] = fStart;
183  }
184 
185  //......................................................................
186  TVector3 Track::Stop() const
187  {
188  return fTraj.back();
189  }
190 
191  //......................................................................
192  // Get direction at the stop point
193  TVector3 Track::StopDir() const
194  {
195  const int ntrajectory_points = NTrajectoryPoints();
196 
197  /// If number of trajectory points is less than 2,
198  /// then can't figure out directions
199  if(ntrajectory_points < 2){
200  const TVector3 nullvector3;
201  return nullvector3;
202  }
203 
204  const TVector3 stop_point = TrajectoryPoint(ntrajectory_points - 1);
205  const TVector3 stop_point_m1 = TrajectoryPoint(ntrajectory_points - 2);
206 
207  const TVector3 diff = stop_point - stop_point_m1;
208 
209  return diff.Unit();
210  }
211 
212  //......................................................................
213  double Track::TotalLength() const
214  {
215  if(fTraj.size() == 1) return 0;
216 
217  // We take the sum of the straight lines between the control points
218  double dist = 0;
219 
220  const int N = fTraj.size();
221  for(int n = 0; n < N-1; ++n){
222  dist += (fTraj[n+1]-fTraj[n]).Mag();
223  }
224 
225  return dist;
226  }
227 
228  //......................................................................
229  double Track::DistanceFromStart(double z) const
230  {
231  if(fTraj.size() == 1) return std::abs(z-fStart.Z());
232 
233  int i0, i1;
234  FindNeighbouringPointIndices(z, i0, i1);
235 
236  // Distance along the track to the first of the neighbouring points
237  double dist = 0;
238  for(int i = 0; i < i0; ++i)
239  dist += (fTraj[i+1]-fTraj[i]).Mag();
240 
241  // The extra accumulated between the points. Also correct for the cases
242  // where we're off the start or end of the track.
243 
244  // Skip this step if the two trajectory points have exactly the same
245  // Z. This would happen if the track were exactly vertical, or more likely
246  // if there's a repeated trajectory point somehow.
247  const double denom = (fTraj[i1]-fTraj[i0]).Unit().Z();
248  if(denom != 0)
249  dist += (z-fTraj[i0].Z()) / denom;
250 
251  return dist;
252  }
253 
254  //......................................................................
255  double Track::DistanceFromEnd(double z) const
256  {
257  if(fTraj.size() == 1) return std::abs(z-fStart.Z());
258 
259  int i0, i1;
260  FindNeighbouringPointIndices(z, i0, i1);
261 
262  // Distance from the end the track to the first of the neighbouring points
263  double dist = 0;
264  for(int i = fTraj.size()-1; i > i1; --i)
265  dist += (fTraj[i]-fTraj[i-1]).Mag();
266 
267  // The extra accumulated between the points. Also correct for the cases
268  // where we're off the start or end of the track.
269 
270  // Skip this step if the two trajectory points have exactly the same
271  // Z. This would happen if the track were exactly vertical, or more likely
272  // if there's a repeated trajectory point somehow.
273  const double denom = (fTraj[i1]-fTraj[i0]).Unit().Z();
274  if(denom != 0)
275  dist += (fTraj[i1].Z()-z) / denom;
276 
277  return dist;
278  }
279 
280  //......................................................................
281  double Track::W(const rb::CellHit* chit) const
282  {
283  if(fTraj.empty()) return 0;
284  /// \todo CJB: do we want to allow the case with only one hit?
285  if(fTraj.size() == 1) return Prong::W(chit);
286 
288  double xyz[3];
289  geom->Plane(chit->Plane())->Cell(chit->Cell())->GetCenter(xyz);
290 
291  InterpolateXY(xyz[2], xyz[0], xyz[1]);
292 
293  if(chit->View() == geo::kX) return xyz[1];
294  if(geom->DetId() == novadaq::cnv::kTESTBEAM) return -xyz[0]; // TB horizontal modules (Y-view) have readout on opposite end as ND/FD
295  return xyz[0];
296  /*
297  // Find the points before and after. This part does still assume that
298  // searching along Z is the best way to find such points.
299  TVector3 prev, next;
300  FindNeighbouringPoints(xyz[2], prev, next);
301 
302  // Solve the problem of the closest approach of the trajectory to the line
303  // made by the cell by projecting down to a 2D trajectory and 2D point for
304  // the cell.
305  TVector3 start2D = prev;
306  TVector3 dir2D = next-prev;
307  TVector3 pt2D(xyz);
308  const int otherView = 1-chit->View();
309  start2D[otherView] = dir2D[otherView] = pt2D[otherView] = 0;
310 
311  // Actually find the closest point
312  TVector3 closest;
313  geo::ClosestApproach(pt2D, start2D, dir2D, closest);
314 
315  // Need to solve for the W position of the line, the exact part that we
316  // projected away above. So just figure out Z and evaluate the original
317  // curve there.
318  const double dt = (closest.Z()-prev.Z())/(next-prev).Z();
319 
320  return ((1-dt)*prev+dt*next)[otherView];
321  */
322  }
323 
324  //......................................................................
325  void Track::InterpolateXY(double z, double& x, double& y) const
326  {
327  // The algorithm here is basically valid for z-directed tracks
328 
329  TVector3 a, b;
330  FindNeighbouringPoints(z, a, b);
331 
332  InterpolatePts(z, a, b, x, y);
333  }
334 
335  //......................................................................
337  const TVector3& a, const TVector3& b,
338  double& x, double& y) const
339  {
340  const double denom = b.Z()-a.Z();
341  const double dt = denom ? (z-a.Z())/denom : 0.5;
342 
343  x = (1-dt)*a.X()+dt*b.X();
344  y = (1-dt)*a.Y()+dt*b.Y();
345  }
346 
347  //......................................................................
348  TVector3 Track::InterpolateDir(double z) const
349  {
350  // The algorithm here is basically valid for z-directed tracks
351 
352  TVector3 a, b;
353  FindNeighbouringPoints(z, a, b);
354 
355  return (b-a).Unit();
356  }
357 
358  //......................................................................
360  {
361  fPrecalcTotalGeV = -1; // Invalidate any cached value
362 
363  fTraj.clear();
364  fTraj.resize(0);
365  }
366  //.....................................................................
367  void Track::RemoveTrajectoryPoint(unsigned int i)
368  {
369  fPrecalcTotalGeV = -1; // Invalidate any cached value
370 
371  if(i > fTraj.size())
372  throw cet::exception("Track") << "attepmt to remove trajectory point "
373  << i << " while there are only " << fTraj.size()
374  << " trajectory points";
375 
376  fTraj.erase(fTraj.begin()+i);
377  }
378 
379  //......................................................................
380  void Track::FindNeighbouringPointIndices(double z, int& i0, int& i1) const
381  {
382  if(fTraj.size() < 2)
383  throw cet::exception("Track") << "attempting to find neighboring points with "
384  << "fewer than 2 trajectory points in track";
385 
386  // Track might be going upstream
387  const int zDir = (fTraj.front().Z() > fTraj.back().Z()) ? -1 : 1;
388 
389  // If hit is off the end of the track extrapolate from the two closest points
390  if( zDir*fTraj.front().Z() >= zDir*z ){
391  i0 = 0;
392  i1 = 1;
393  return;
394  }
395  if( zDir*fTraj.back().Z() <= zDir*z ){
396  i0 = fTraj.size()-2;
397  i1 = fTraj.size()-1;
398  return;
399  }
400 
401  bool found = false;
402 
403  // We could do this more efficiently with std::[upper|lower]_bound, but the
404  // presence of reversed tracks complicates that slightly.
405  for(size_t n = 0; n < fTraj.size()-1; ++n){
406  if( zDir*fTraj[n] .Z() <= zDir*z &&
407  zDir*fTraj[n+1].Z() >= zDir*z ){
408  // If we've already found a pair of neighbours then this new pair is
409  // ambiguous which is bad. Probably it means the track reverses
410  // direction wrt Z.
411  if(found){
412  // But, it might also mean the value passed in was right on a
413  // trajectory point. In which case there isn't really a problem.
414  // Either choice will end up just giving the value at this point.
415  if(std::abs(fTraj[n].Z()-z) > 1e-3 && std::abs(fTraj[n+1].Z()-z) > 1e-3){
416  LOG_DEBUG("Track")
417  << "Choice of neighbouring trajectory points is ambiguous";
418  }
419  }
420  i0 = n;
421  i1 = n+1;
422  found = true;
423  break;
424  }
425  }
426 
427  if(!found)
428  throw cet::exception("Track") << "unable to find neighboring point index";
429  }
430 
431  //......................................................................
433  TVector3& pt0, TVector3& pt1) const
434  {
435  int i0, i1;
436  FindNeighbouringPointIndices(z, i0, i1);
437  pt0 = fTraj[i0];
438  pt1 = fTraj[i1];
439  }
440 
441  //......................................................................
443  {
444  if(!Is2D() || !trk.Is2D())
445  throw cet::exception("Track") << "attempting to zip a tracks with unknown dimensions; this track: "
446  << this->View() << " zipping track: " << trk.View();
447 
448  if(this->NTrajectoryPoints() < 2 || trk.NTrajectoryPoints() < 2)
449  throw cet::exception("Track") << "too few points to attempt ZipWith: "
450  << this->NTrajectoryPoints() << " "
451  << trk.NTrajectoryPoints();
452 
453  const rb::Track* xtrk = this;
454  const rb::Track* ytrk = &trk;
455 
456  if(this->View() != geo::kX) std::swap(xtrk, ytrk);
457 
458  if(xtrk->View() != geo::kX || ytrk->View() != geo::kY)
459  throw cet::exception("Track") << "views of 2D tracks are wrong";
460 
461 
462  std::vector<TVector3> zipTraj;
463  zipTraj.reserve(xtrk->NTrajectoryPoints() + ytrk->NTrajectoryPoints());
464  // Give each point a 3rd coordinate interpolated from the other track
465  for(const TVector3& pt: xtrk->fTraj){
466  double junk, y;
467  ytrk->InterpolateXY(pt.Z(), junk, y);
468  zipTraj.emplace_back(pt.X(), y, pt.Z());
469  }
470  for(const TVector3& pt: ytrk->fTraj){
471  double junk, x;
472  xtrk->InterpolateXY(pt.Z(), x, junk);
473  zipTraj.emplace_back(x, pt.Y(), pt.Z());
474  }
475 
476  // Sort output points upstream to downstream. User could reverse result if
477  // they wanted.
478  std::sort(zipTraj.begin(), zipTraj.end(),
479  [](const TVector3& a, const TVector3& b)
480  {
481  return a.Z() < b.Z();
482  });
483 
484  Track ret(xtrk->AllCells(), zipTraj[0], zipTraj[1]-zipTraj[0]);
485  ret.Add(ytrk->AllCells());
486 
487  // First point is already stored as the vertex
488  for(size_t p = 1; p < zipTraj.size(); ++p){
489  ret.AppendTrajectoryPoint(zipTraj[p]);
490  }
491 
492  return ret;
493  }
494 
495  //-----------------------------------------------------------------------------
496  std::map<unsigned int, TVector3> Track::PlaneDirMap() const
497  {
498  std::map<unsigned int, TVector3> planeDirMap;
499 
500  // loop over the planes in the track and get the points on either side
501  // of the plane z
503 
504  // looping over the cells this way does the X and Y
505  // view separately, but the map orders by plane
506  // number using the < operator, so the entries
507  // in the map are linearly increasing, with X and Y
508  // views mixed together.
509  for(size_t c = 0; c < NCell(); ++c){
510  unsigned int plane = Cell(c)->Plane();
511 
512  // check if we already have the direction at this plane
513  // no need to overwrite it if we do.
514  if(planeDirMap.count(plane) > 0) continue;
515 
516  // get the z position for this cell
517  double xyz[3] = {0.};
518  geom->Plane(plane)->Cell(Cell(c)->Cell())->GetCenter(xyz);
519 
520  // find the direction based on the points on either
521  // side of this z location.
522  // set the direction of travel as the direction
523  // at this plane.
524  planeDirMap[plane] = InterpolateDir(xyz[2]);
525 
526  }// end loop over cells
527 
528  return planeDirMap;
529  }
530 
531  //......................................................................
533  const art::Ptr<rb::Track>& b)
534  {
535  return a->MinTNS() < b->MinTNS();
536  }
537 
538  //......................................................................
540  {
541  std::stable_sort(t.begin(), t.end(), CompareTracksByTime);
542  }
543 
544  //......................................................................
546  {
547  std::stable_sort(t.begin(), t.end(), CompareTracksByTime);
548  }
549 
550 } // end namespace rb
551 //////////////////////////////////////////////////////////////////////////////
Track()
Definition: Track.h:27
void FindNeighbouringPointIndices(double z, int &i0, int &i1) const
Definition: Track.cxx:380
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
virtual void SetStart(TVector3 start)
Definition: Track.cxx:168
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
double fPrecalcTotalGeV
-1 = uninitialized
Definition: Cluster.h:324
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Prong.cxx:135
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Track.cxx:281
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
iterator begin()
Definition: PtrVector.h:223
geo::View_t View() const
Definition: CellHit.h:41
const CellGeo * Cell(int icell) const
Definition: PlaneGeo.h:48
const char * p
Definition: xmltok.h:285
Vertical planes which measure X.
Definition: PlaneGeo.h:28
A collection of associated CellHits.
Definition: Cluster.h:47
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
A rb::Prong with full reconstructed trajectory.
Definition: Track.h:20
const PlaneGeo * Plane(unsigned int i) const
void SetWeight(unsigned int globalIdx, double weight)
Set weight of the cell at this index.
Definition: Cluster.cxx:327
virtual ~Track()
Definition: Track.cxx:107
virtual TVector3 Start() const
Definition: Prong.h:73
virtual void SetStart(TVector3 start)
Definition: Prong.cxx:75
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
double dist
Definition: runWimpSim.h:113
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
Float_t Z
Definition: plot.C:38
virtual double TotalLength() const
Length (cm) of all the track segments.
Definition: Track.cxx:213
unsigned short Cell() const
Definition: CellHit.h:40
Track finder for cosmic rays.
Definition: Cand.cxx:23
double MinTNS() const
Definition: Cluster.cxx:482
void SortTracksByTime(std::vector< art::Ptr< rb::Track >> &t)
Sort t in time order (earliest to latest).
Definition: Track.cxx:539
void InterpolatePts(double z, const TVector3 &a, const TVector3 &b, double &x, double &y) const
Helper function for InterpolateXY and ZipWith.
Definition: Track.cxx:336
double dz[NP][NC]
unsigned int NCell() const
Number of cells in either view.
Definition: Cluster.h:110
const double a
geo::View_t fView
view this cluster is in
Definition: Cluster.h:316
void AppendTrajectoryPoint(TVector3 pt)
Definition: Track.cxx:112
void FindNeighbouringPoints(double z, TVector3 &pt0, TVector3 &pt1) const
Helper function for W.
Definition: Track.cxx:432
novadaq::cnv::DetId DetId() const
Prefer ds::DetectorService::DetId() instead.
Definition: GeometryBase.h:243
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
iterator end()
Definition: PtrVector.h:237
void ClearTrajectoryPoints()
Forget about all trajectory points.
Definition: Track.cxx:359
Collect Geo headers and supply basic geometry functions.
Var weights
void SetID(int id)
Definition: Cluster.h:74
void PrependTrajectoryPoint(TVector3 pt)
Support constructing tracks backwards.
Definition: Track.cxx:134
TVector3 Unit() const
Perform a "2 point" Hough transform on a collection of hits.
z
Definition: test.py:28
TVector3 fStart
Start location (xyz, cm)
Definition: Prong.h:111
virtual double DistanceFromStart(double z) const
Definition: Track.cxx:229
std::map< unsigned int, TVector3 > PlaneDirMap() const
map of the direction cosines at each plane (ie z)
Definition: Track.cxx:496
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
A rawdata::RawDigit with channel information decoded.
Definition: CellHit.h:27
A Cluster with defined start position and direction.
Definition: Prong.h:19
virtual TVector3 StopDir() const
Get direction at the stop point.
Definition: Track.cxx:193
TDirectory * dir
Definition: macro.C:5
bool CompareTracksByTime(const art::Ptr< rb::Track > &a, const art::Ptr< rb::Track > &b)
Helper for SortTracksByTime. Is a earlier than b?
Definition: Track.cxx:532
void geom(int which=0)
Definition: geom.C:163
float Mag() const
const hit & b
Definition: hits.cxx:21
virtual TVector3 InterpolateDir(double z) const
Definition: Track.cxx:348
void RemoveTrajectoryPoint(unsigned int i)
Remove the ith trajectory point from the track.
Definition: Track.cxx:367
std::vector< TVector3 > fTraj
Definition: Track.h:160
virtual void InterpolateXY(double z, double &x, double &y) const
Definition: Track.cxx:325
TVector3 Stop() const
Position of the final trajectory point.
Definition: Track.cxx:186
Float_t e
Definition: plot.C:35
virtual double DistanceFromEnd(double z) const
Definition: Track.cxx:255
Float_t w
Definition: plot.C:20
Encapsulate the geometry of one entire detector (near, far, ndos)