LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
Public Member Functions | Private Member Functions | Private Attributes | List of all members
trkf::TrackTrajectoryAlg Class Reference

#include "TrackTrajectoryAlg.h"

Public Member Functions

 TrackTrajectoryAlg ()
 
virtual ~TrackTrajectoryAlg ()
 
void TrackTrajectory (std::array< std::vector< geo::WireID >, 3 > trkWID, std::array< std::vector< double >, 3 > trkX, std::array< std::vector< double >, 3 > trkXErr, std::vector< TVector3 > &TrajPos, std::vector< TVector3 > &TrajDir)
 

Private Member Functions

void ShortTrackTrajectory (std::array< std::vector< geo::WireID >, 3 > trkWID, std::array< std::vector< double >, 3 > trkX, std::array< std::vector< double >, 3 > trkXErr, std::vector< TVector3 > &TrajPos, std::vector< TVector3 > &TrajDir)
 

Private Attributes

art::ServiceHandle< geo::Geometry > geom
 
double minX
 
unsigned short minXPln
 
double maxX
 
unsigned short maxXPln
 
bool prt
 
unsigned short fMaxTrajPoints
 
double fHitWidthFactor
 
TrackLineFitAlg fTrackLineFitAlg
 

Detailed Description

Definition at line 28 of file TrackTrajectoryAlg.h.

Constructor & Destructor Documentation

trkf::TrackTrajectoryAlg::TrackTrajectoryAlg ( )

Definition at line 24 of file TrackTrajectoryAlg.cxx.

24 { }
trkf::TrackTrajectoryAlg::~TrackTrajectoryAlg ( )
virtual

Definition at line 26 of file TrackTrajectoryAlg.cxx.

26 { }

Member Function Documentation

void trkf::TrackTrajectoryAlg::ShortTrackTrajectory ( std::array< std::vector< geo::WireID >, 3 >  trkWID,
std::array< std::vector< double >, 3 >  trkX,
std::array< std::vector< double >, 3 >  trkXErr,
std::vector< TVector3 > &  TrajPos,
std::vector< TVector3 > &  TrajDir 
)
private

Definition at line 289 of file TrackTrajectoryAlg.cxx.

References dir, fHitWidthFactor, geom, geo::GeometryCore::IntersectionPoint(), maxX, maxXPln, minX, minXPln, prt, y, and z.

Referenced by TrackTrajectory().

293  {
294  // Make a short track trajectory composed of a space point at each end.
295 
296  std::vector<unsigned short> usePlns;
297  double y, z, endX, endY, endZ;
298  unsigned short ipl, iht, nend, jpl, iPlane, jPlane;
299  unsigned int iWire, jWire;
300 
301  // do the min X end first
302  float xCut;
303  for(unsigned short nit = 0; nit < 5; ++nit) {
304  xCut = minX + fHitWidthFactor * trkXErr[minXPln][0];
305  for(ipl = 0; ipl < 3; ++ipl) {
306  if(trkX[ipl].size() == 0) continue;
307  if(trkX[ipl][0] < xCut) usePlns.push_back(ipl);
308  } // ipl
309  if(usePlns.size() >= 2) break;
310  fHitWidthFactor += 5;
311  }
312  // Not enough information to find a space point
313  if(prt) mf::LogVerbatim("TTA")<<"ShortTrack minX end "<<xCut<<" usePlns size "<<usePlns.size();
314  if(usePlns.size() < 2) return;
315  endY = 0; endZ = 0; nend = 0;
316  iht = 0;
317  ipl = usePlns[0];
318  endX = trkX[ipl][iht];
319  iPlane = trkWID[ipl][iht].Plane;
320  iWire = trkWID[ipl][iht].Wire;
321  unsigned int cstat = trkWID[ipl][iht].Cryostat;
322  unsigned int tpc = trkWID[ipl][iht].TPC;
323  for(unsigned short jj = 1; jj < usePlns.size(); ++jj) {
324  jpl = usePlns[jj];
325  endX += trkX[jpl][iht];
326  jPlane = trkWID[jpl][iht].Plane;
327  jWire = trkWID[jpl][iht].Wire;
328  geom->IntersectionPoint(iWire, jWire, iPlane, jPlane, cstat, tpc, y, z);
329  endY += y;
330  endZ += z;
331  ++nend;
332  } // ii
333  // nend is the number of y,z values. There are (nend + 1) X values
334  TVector3 xyz;
335  xyz(0) = endX / (float)(nend + 1);
336  xyz(1) = endY / (float)nend;
337  xyz(2) = endZ / (float)nend;
338  if(prt) mf::LogVerbatim("TTA")<<" xyz "<<xyz(0)<<" "<<xyz(1)<<" "<<xyz(2);
339  TrajPos.push_back(xyz);
340 
341  // do the same for the other end
342  fHitWidthFactor = 5;
343 // xCut = maxX - fHitWidthFactor * trkXErr[maxXPln][0];
344  usePlns.clear();
345  for(unsigned short nit = 0; nit < 5; ++nit) {
346  xCut = maxX - fHitWidthFactor * trkXErr[maxXPln][0];
347  for(ipl = 0; ipl < 3; ++ipl) {
348  if(trkX[ipl].size() == 0) continue;
349  iht = trkX[ipl].size() - 1;
350  if(trkX[ipl][iht] > xCut) usePlns.push_back(ipl);
351  } // ipl
352  if(usePlns.size() >= 2) break;
353  fHitWidthFactor += 5;
354  }
355  // Not enough information to find a space point
356  if(prt) mf::LogVerbatim("TTA")<<"ShortTrack maxX end "<<xCut<<" usePlns size "<<usePlns.size();
357  if(usePlns.size() < 2) {
358  TrajPos.clear();
359  TrajDir.clear();
360  return;
361  }
362  endY = 0; endZ = 0; nend = 0;
363  ipl = usePlns[0];
364  iht = trkX[ipl].size() - 1;
365  endX = trkX[ipl][iht];
366  iPlane = trkWID[ipl][iht].Plane;
367  iWire = trkWID[ipl][iht].Wire;
368  for(unsigned short jj = 1; jj < usePlns.size(); ++jj) {
369  jpl = usePlns[jj];
370  iht = trkX[jpl].size() - 1;
371  endX += trkX[jpl][iht];
372  jPlane = trkWID[jpl][iht].Plane;
373  jWire = trkWID[jpl][iht].Wire;
374  geom->IntersectionPoint(iWire, jWire, iPlane, jPlane, cstat, tpc, y, z);
375  endY += y;
376  endZ += z;
377  ++nend;
378  } // ii
379  // nend is the number of y,z values. There are nend + 1 X values
380  xyz(0) = endX / (float)(nend + 1);
381  xyz(1) = endY / (float)nend;
382  xyz(2) = endZ / (float)nend;
383  if(prt) mf::LogVerbatim("TTA")<<" xyz "<<xyz(0)<<" "<<xyz(1)<<" "<<xyz(2);
384  TrajPos.push_back(xyz);
385  TVector3 dir = TrajPos[1] - TrajPos[0];
386  dir = dir.Unit();
387  TrajDir.push_back(dir);
388  TrajDir.push_back(dir);
389 
390  if(prt) mf::LogVerbatim("TTA")<<">>>> Short track ("<<TrajPos[0](0)<<", "<<TrajPos[0](1)<<", "<<TrajPos[0](2)<<") to ("<<TrajPos[1](0)<<", "<<TrajPos[1](1)<<", "<<TrajPos[1](2)<<")";
391  }
trkf::TrackTrajectoryAlg::geom
art::ServiceHandle< geo::Geometry > geom
Definition: TrackTrajectoryAlg.h:42
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:219
y
Float_t y
Definition: compare.C:6
z
Double_t z
Definition: plot.C:279
geo::GeometryCore::IntersectionPoint
bool IntersectionPoint(geo::WireID const &wid1, geo::WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
Definition: GeometryCore.cxx:1749
dir
TDirectory * dir
Definition: macro.C:5
void trkf::TrackTrajectoryAlg::TrackTrajectory ( std::array< std::vector< geo::WireID >, 3 >  trkWID,
std::array< std::vector< double >, 3 >  trkX,
std::array< std::vector< double >, 3 >  trkXErr,
std::vector< TVector3 > &  TrajPos,
std::vector< TVector3 > &  TrajDir 
)

Definition at line 29 of file TrackTrajectoryAlg.cxx.

References evd::details::begin(), dir, evd::details::end(), fHitWidthFactor, fMaxTrajPoints, fTrackLineFitAlg, maxX, maxXPln, minX, minXPln, prt, ShortTrackTrajectory(), and trkf::TrackLineFitAlg::TrkLineFit().

33  {
34 
35  // Make a track trajectory (position, direction) and return it in the TrajPos
36  // and TrajDir vectors. The track hits are received as 3 vectors (one vector per wire plane)
37  // of Wire ID's, hit X and X position errors. The X position errors are used to 1) determine
38  // the significance of the X difference between the beginning and end of the track trajectory
39  // 2) determine the number of trajectory points and 3) weight the track line fit in TrackLineFitAlg.
40  // This code assumes that hits at each end (e.g. trkXW[Plane][0]) of the vectors define the end
41  // points of the trajectory. The ordering of planes in the array is irrelevant since the
42  // plane number is extracted from the WireID. This algorithm will return with a failed condition
43  // (TrajPos, TrajDir size = 0) if there are fewer than 2 planes with hits at each end that are less
44  // than fHitWidthFactor * trkXErr apart. Valid and invalid conditions are shown schematically below
45  // where a . represents hits that are separated by X values > fHitWidthFactor * trkXErr
46  //
47  // minX maxX maxX minX minX maxX
48  // Pln0 .................... Pln0 .................... Pln0 ...................
49  // Pln1 ............... Pln1 .................... Pln1
50  // Pln2 .................... Pln2 ................ Pln2 ...................
51  // VALID VALID VALID - no hits in one plane is OK
52  //
53  // minX maxX
54  // Pln0 .................
55  // Pln1 ...............
56  // Pln2 ....................
57  // NOT VALID - Only one plane has a hit at MaxX
58 
59 
60  TrajPos.clear();
61  TrajDir.clear();
62 
63  prt = false;
64 
65  unsigned short minLen = 9999;
66  unsigned short maxLen = 0;
67  unsigned short nPlnsWithHits = 0;
68  unsigned short ipl, aPlane = 3;
69  for(ipl = 0; ipl < 3; ++ipl) {
70  if(trkX[ipl].size() == 0) continue;
71  ++nPlnsWithHits;
72  if(trkX[ipl].size() < minLen) {
73  minLen = trkX[ipl].size();
74  }
75  if(trkX[ipl].size() > maxLen) {
76  maxLen = trkX[ipl].size();
77  aPlane = ipl;
78  }
79  } // ipl
80  if(prt) mf::LogVerbatim("TTA")<<"trkX sizes "<<trkX[0].size()<<" "<<trkX[1].size()<<" "<<trkX[2].size()<<" "<<" nPlnsWithHits "<<nPlnsWithHits;
81  if(nPlnsWithHits < 2) return;
82  if(aPlane > 2) return;
83 
84  fMaxTrajPoints = 100;
85  fHitWidthFactor = 5.;
86 
87  unsigned short iht;
88 
89  // reverse the order of the hit arrays if necessary so that the minimum X end is at trk[][0] end.
90  // We will use posX to reverse the trajectory later if necessary
91  bool posX = true;
92  iht = trkX[aPlane].size() - 1;
93  if(trkX[aPlane][0] > trkX[aPlane][iht]) {
94  posX = false;
95  for(ipl = 0; ipl < 3; ++ipl) {
96  if(trkX[ipl].size() == 0) continue;
97  std::reverse(trkWID[ipl].begin(), trkWID[ipl].end());
98  std::reverse(trkX[ipl].begin(), trkX[ipl].end());
99  std::reverse(trkXErr[ipl].begin(), trkXErr[ipl].end());
100  } // ipl
101  if(prt) mf::LogVerbatim("TTA")<<"Swapped order";
102  } // posX check
103 
104  if(prt) {
105  mf::LogVerbatim myprt("TTA");
106  myprt<<"TrkXW end0";
107  for(ipl = 0; ipl < 3; ++ipl) {
108  if(trkX[ipl].size() == 0) continue;
109  myprt<<" "<<trkX[ipl][0];
110  }
111  myprt<<"\n";
112  myprt<<"TrkXW end1";
113  for(ipl = 0; ipl < 3; ++ipl) {
114  if(trkX[ipl].size() == 0) continue;
115  iht = trkX[ipl].size() - 1;
116  myprt<<" "<<trkX[ipl][iht];
117  }
118  myprt<<"\n";
119  }
120 
121  // find the min/max X
122  minX = 1E6;
123  minXPln = 0;
124  maxX = -1E6;
125  maxXPln = 0;
126 
127  for(unsigned short ipl = 0; ipl < 3; ++ipl) {
128  if(trkX[ipl].size() == 0) continue;
129  if(trkX[ipl][0] < minX) {
130  minX = trkX[ipl][0];
131  minXPln = ipl;
132  }
133  iht = trkX[ipl].size() - 1;
134  if(trkX[ipl][iht] > maxX) {
135  maxX = trkX[ipl][iht];
136  maxXPln = ipl;
137  }
138  } // ipl
139 
140  if(prt) mf::LogVerbatim("TTA")<<"minX "<<minX<<" in plane "<<minXPln<<" maxX "<<maxX<<" in plane "<<maxXPln;
141 
142  // estimate the number of trajectory points we will want based on the delta T and the
143  // hit errors
144  double aveHitXErr = 0;
145  unsigned short nHit = 0;
146  for(ipl = 0; ipl < 3; ++ipl) {
147  for(iht = 0; iht < trkXErr[ipl].size(); ++iht) {
148  aveHitXErr += trkXErr[ipl][iht];
149  ++nHit;
150  }
151  } // ipl
152  aveHitXErr /= (double)nHit;
153 
154  unsigned short npt = (maxX - minX) / (1 * aveHitXErr);
155  if(npt < 2) npt = 2;
156  if(npt > maxLen) npt = maxLen;
157  if(npt > fMaxTrajPoints) npt = fMaxTrajPoints;
158  if(prt) mf::LogVerbatim("TTA")<<" aveHitXErr "<<aveHitXErr<<" number of traj points ";
159 
160  double maxBinX = (maxX - minX) / (double)(npt - 1);
161  double binX = maxBinX;
162  double xOrigin = minX;
163 
164  TVector3 xyz, dir;
165  // working vectors passed to TrackLineFitAlg
166  std::vector<geo::WireID> hitWID;
167  std::vector<double> hitX;
168  std::vector<double> hitXErr;
169  float ChiDOF;
170 
171  // make a short track trajectory (end points only) to use in case an error occurs
172  std::vector<TVector3> STPos;
173  std::vector<TVector3> STDir;
174  ShortTrackTrajectory(trkWID, trkX, trkXErr, STPos, STDir);
175 
176  if(STPos.size() != 2 || STDir.size() != STPos.size()) {
177  TrajPos.clear();
178  TrajDir.clear();
179  if(posX) {
180  for(ipl = 0; ipl < 3; ++ipl) {
181  if(trkX[ipl].size() == 0) continue;
182  std::reverse(trkWID[ipl].begin(), trkWID[ipl].end());
183  std::reverse(trkX[ipl].begin(), trkX[ipl].end());
184  std::reverse(trkXErr[ipl].begin(), trkXErr[ipl].end());
185  } // ipl
186  } // posX
187  return;
188  } // bad STPos, STDir
189 
190  if(prt) {
191  mf::LogVerbatim("TTA")<<"STPos";
192  for(unsigned short ii = 0; ii < STPos.size(); ++ii) mf::LogVerbatim("TTA")<<ii<<" "<<std::fixed<<std::setprecision(1)<<STPos[ii](0)<<" "<<STPos[ii](1)<<" "<<STPos[ii](2);
193  }
194 
195  if(maxLen < 4 || npt < 2) {
196  TrajPos = STPos;
197  TrajDir = STDir;
198  if(!posX) {
199  // reverse everything
200  std::reverse(TrajPos.begin(), TrajPos.end());
201  std::reverse(TrajDir.begin(), TrajDir.end());
202  for(ipl = 0; ipl < 3; ++ipl) {
203  if(trkX[ipl].size() == 0) continue;
204  std::reverse(trkWID[ipl].begin(), trkWID[ipl].end());
205  std::reverse(trkX[ipl].begin(), trkX[ipl].end());
206  std::reverse(trkXErr[ipl].begin(), trkXErr[ipl].end());
207  } // ipl
208  }
209  return;
210  } // maxLen < 4
211 
212  // Start index of hits to include in the next fit
213  std::array<unsigned short, 3> hStart;
214  for(ipl = 0; ipl < 3; ++ipl) hStart[ipl] = 0;
215 
216  bool gotLastPoint = false;
217  for(unsigned short ipt = 0; ipt < npt + 1; ++ipt) {
218  hitWID.clear();
219  hitX.clear();
220  hitXErr.clear();
221  for(ipl = 0; ipl < 3; ++ipl) {
222  for(iht = hStart[ipl]; iht < trkX[ipl].size(); ++iht) {
223  if(trkX[ipl][iht] < xOrigin - binX) continue;
224  if(trkX[ipl][iht] > xOrigin + binX) break;
225  hitWID.push_back(trkWID[ipl][iht]);
226  hitX.push_back(trkX[ipl][iht]);
227  hitXErr.push_back(trkXErr[ipl][iht]);
228  hStart[ipl] =iht;
229  } // iht
230  } // ipl
231  if(prt) mf::LogVerbatim("TTA")<<"ipt "<<ipt<<" xOrigin "<<xOrigin<<" binX "<<binX<<" hitX size "<<hitX.size();
232  if(hitX.size() > 3) {
233  fTrackLineFitAlg.TrkLineFit(hitWID, hitX, hitXErr, xOrigin, xyz, dir, ChiDOF);
234  if(prt) mf::LogVerbatim("TTA")<<" xyz "<<xyz(0)<<" "<<xyz(1)<<" "<<xyz(2)<<" dir "<<dir(0)<<" "<<dir(1)<<" "<<dir(2)<<" ChiDOF "<<ChiDOF<<" hitX size "<<hitX.size();
235  } else if(ipt == 0 && STPos.size() == 2) {
236  // failure on the first traj point. Use STPos
237  xyz = STPos[0];
238  dir = STDir[0];
239  }
240  if(prt && hitX.size() < 4) mf::LogVerbatim("TTA")<<"\n";
241  if(xOrigin >= maxX) break;
242  // tweak xOrigin if we are close to the end
243  if(maxX - xOrigin < binX) {
244  xOrigin = maxX;
245  } else {
246  xOrigin += binX;
247  }
248  if(ChiDOF < 0 || ChiDOF > 100) continue;
249  TrajPos.push_back(xyz);
250  TrajDir.push_back(dir);
251  if(ipt == npt) gotLastPoint = true;
252  } // ipt
253  if(prt) {
254  mf::LogVerbatim("TTA")<<"gotLastPoint "<<gotLastPoint<<" TTA Traj \n";
255  for(unsigned short ii = 0; ii < TrajPos.size(); ++ii) mf::LogVerbatim("TTA")<<ii<<" "<<std::fixed<<std::setprecision(1)<<TrajPos[ii](0)<<" "<<TrajPos[ii](1)<<" "<<TrajPos[ii](2);
256  }
257 
258  if(!gotLastPoint && STPos.size() == 2) {
259  // failure on the last point. Use STPos, etc
260  TrajPos.push_back(STPos[1]);
261  TrajDir.push_back(STDir[1]);
262  }
263 
264  if(TrajPos.size() < 2) {
265  TrajPos = STPos;
266  TrajDir = STDir;
267  }
268 
269  if(!posX) {
270  // reverse everything
271  std::reverse(TrajPos.begin(), TrajPos.end());
272  std::reverse(TrajDir.begin(), TrajDir.end());
273  for(ipl = 0; ipl < 3; ++ipl) {
274  if(trkX[ipl].size() == 0) continue;
275  std::reverse(trkWID[ipl].begin(), trkWID[ipl].end());
276  std::reverse(trkX[ipl].begin(), trkX[ipl].end());
277  std::reverse(trkXErr[ipl].begin(), trkXErr[ipl].end());
278  } // ipl
279 
280  } // !posX
281  if(prt) {
282  mf::LogVerbatim("TTA")<<"TTA Traj2\n";
283  for(unsigned short ii = 0; ii < TrajPos.size(); ++ii) mf::LogVerbatim("TTA")<<ii<<" "<<std::fixed<<std::setprecision(1)<<TrajPos[ii](0)<<" "<<TrajPos[ii](1)<<" "<<TrajPos[ii](2);
284  }
285 
286  } // TrackTrajectoryAlg
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:219
trkf::TrackTrajectoryAlg::fTrackLineFitAlg
TrackLineFitAlg fTrackLineFitAlg
Definition: TrackTrajectoryAlg.h:54
mf::MaybeLogger_
Definition: MessageLogger.h:114
trkf::TrackLineFitAlg::TrkLineFit
void TrkLineFit(std::vector< geo::WireID > &hitWID, std::vector< double > &hitX, std::vector< double > &hitXErr, double XOrigin, TVector3 &Pos, TVector3 &Dir, float &ChiDOF)
Definition: TrackLineFitAlg.cxx:30
evd::details::begin
std::vector< evd::details::RawDigitInfo_t >::const_iterator begin(RawDigitCacheDataClass const &cache)
Definition: RawDataDrawer.cxx:314
trkf::TrackTrajectoryAlg::ShortTrackTrajectory
void ShortTrackTrajectory(std::array< std::vector< geo::WireID >, 3 > trkWID, std::array< std::vector< double >, 3 > trkX, std::array< std::vector< double >, 3 > trkXErr, std::vector< TVector3 > &TrajPos, std::vector< TVector3 > &TrajDir)
Definition: TrackTrajectoryAlg.cxx:289
dir
TDirectory * dir
Definition: macro.C:5
evd::details::end
std::vector< evd::details::RawDigitInfo_t >::const_iterator end(RawDigitCacheDataClass const &cache)
Definition: RawDataDrawer.cxx:317

Member Data Documentation

double trkf::TrackTrajectoryAlg::fHitWidthFactor
private

Definition at line 51 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

unsigned short trkf::TrackTrajectoryAlg::fMaxTrajPoints
private

Definition at line 50 of file TrackTrajectoryAlg.h.

Referenced by TrackTrajectory().

TrackLineFitAlg trkf::TrackTrajectoryAlg::fTrackLineFitAlg
private

Definition at line 54 of file TrackTrajectoryAlg.h.

Referenced by TrackTrajectory().

art::ServiceHandle<geo::Geometry> trkf::TrackTrajectoryAlg::geom
private

Definition at line 42 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory().

double trkf::TrackTrajectoryAlg::maxX
private

Definition at line 46 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

unsigned short trkf::TrackTrajectoryAlg::maxXPln
private

Definition at line 47 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

double trkf::TrackTrajectoryAlg::minX
private

Definition at line 44 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

unsigned short trkf::TrackTrajectoryAlg::minXPln
private

Definition at line 45 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

bool trkf::TrackTrajectoryAlg::prt
private

Definition at line 48 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

The documentation for this class was generated from the following files: