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

#include "TrackTrajectoryAlg.h"

Public Member Functions

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

geo::WireReadoutGeom const * wireReadoutGeom = &art::ServiceHandle<geo::WireReadout>()->Get()
 
double minX
 
unsigned short minXPln
 
double maxX
 
unsigned short maxXPln
 
bool prt
 
unsigned short fMaxTrajPoints
 
double fHitWidthFactor
 
TrackLineFitAlg fTrackLineFitAlg
 

Detailed Description

Definition at line 26 of file TrackTrajectoryAlg.h.

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 297 of file TrackTrajectoryAlg.cxx.

References dir, fHitWidthFactor, geo::WireIDIntersection::invalid(), maxX, maxXPln, minX, minXPln, prt, util::size(), geo::WireReadoutGeom::WireIDsIntersect(), and wireReadoutGeom.

Referenced by TrackTrajectory().

302  {
303  // Make a short track trajectory composed of a space point at each end.
304 
305  std::vector<unsigned short> usePlns;
306 
307  // do the min X end first
308  float xCut = minX + fHitWidthFactor * trkXErr[minXPln][0];
309  for (unsigned short nit = 0; nit < 5; ++nit) {
310  for (unsigned int ipl = 0; ipl < 3; ++ipl) {
311  if (trkX[ipl].size() == 0) continue;
312  if (trkX[ipl][0] < xCut) usePlns.push_back(ipl);
313  } // ipl
314  if (usePlns.size() >= 2) break;
315  fHitWidthFactor += 5;
316  }
317  // Not enough information to find a space point
318  if (prt)
319  mf::LogVerbatim("TTA") << "ShortTrack minX end " << xCut << " usePlns size "
320  << usePlns.size();
321  if (usePlns.size() < 2) return;
322  double endY = 0;
323  double endZ = 0;
324  unsigned short nend = 0;
325  unsigned short iht = 0;
326  unsigned short ipl = usePlns[0];
327  double endX = trkX[ipl][iht];
328  unsigned short iPlane = trkWID[ipl][iht].Plane;
329  unsigned int iWire = trkWID[ipl][iht].Wire;
330  unsigned int cstat = trkWID[ipl][iht].Cryostat;
331  unsigned int tpc = trkWID[ipl][iht].TPC;
332  for (unsigned short jj = 1; jj < usePlns.size(); ++jj) {
333  unsigned short jpl = usePlns[jj];
334  endX += trkX[jpl][iht];
335  unsigned short jPlane = trkWID[jpl][iht].Plane;
336  unsigned int jWire = trkWID[jpl][iht].Wire;
337  auto intersection = wireReadoutGeom
338  ->WireIDsIntersect(geo::WireID{cstat, tpc, iPlane, iWire},
339  geo::WireID{cstat, tpc, jPlane, jWire})
340  .value_or(geo::WireIDIntersection::invalid());
341  endY += intersection.y;
342  endZ += intersection.z;
343  ++nend;
344  } // ii
345  // nend is the number of y,z values. There are (nend + 1) X values
346  TVector3 xyz;
347  xyz(0) = endX / (float)(nend + 1);
348  xyz(1) = endY / (float)nend;
349  xyz(2) = endZ / (float)nend;
350  if (prt) mf::LogVerbatim("TTA") << " xyz " << xyz(0) << " " << xyz(1) << " " << xyz(2);
351  TrajPos.push_back(xyz);
352 
353  // do the same for the other end
354  fHitWidthFactor = 5;
355  usePlns.clear();
356  for (unsigned short nit = 0; nit < 5; ++nit) {
357  xCut = maxX - fHitWidthFactor * trkXErr[maxXPln][0];
358  for (ipl = 0; ipl < 3; ++ipl) {
359  if (trkX[ipl].size() == 0) continue;
360  iht = trkX[ipl].size() - 1;
361  if (trkX[ipl][iht] > xCut) usePlns.push_back(ipl);
362  } // ipl
363  if (usePlns.size() >= 2) break;
364  fHitWidthFactor += 5;
365  }
366  // Not enough information to find a space point
367  if (prt)
368  mf::LogVerbatim("TTA") << "ShortTrack maxX end " << xCut << " usePlns size "
369  << usePlns.size();
370  if (usePlns.size() < 2) {
371  TrajPos.clear();
372  TrajDir.clear();
373  return;
374  }
375  endY = 0;
376  endZ = 0;
377  nend = 0;
378  ipl = usePlns[0];
379  iht = trkX[ipl].size() - 1;
380  endX = trkX[ipl][iht];
381  iPlane = trkWID[ipl][iht].Plane;
382  iWire = trkWID[ipl][iht].Wire;
383  for (unsigned short jj = 1; jj < usePlns.size(); ++jj) {
384  unsigned short jpl = usePlns[jj];
385  iht = trkX[jpl].size() - 1;
386  endX += trkX[jpl][iht];
387  unsigned short jPlane = trkWID[jpl][iht].Plane;
388  unsigned int jWire = trkWID[jpl][iht].Wire;
389  auto intersection = wireReadoutGeom
390  ->WireIDsIntersect(geo::WireID{cstat, tpc, iPlane, iWire},
391  geo::WireID{cstat, tpc, jPlane, jWire})
392  .value_or(geo::WireIDIntersection::invalid());
393  endY += intersection.y;
394  endZ += intersection.z;
395  ++nend;
396  } // ii
397  // nend is the number of y,z values. There are nend + 1 X values
398  xyz(0) = endX / (float)(nend + 1);
399  xyz(1) = endY / (float)nend;
400  xyz(2) = endZ / (float)nend;
401  if (prt) mf::LogVerbatim("TTA") << " xyz " << xyz(0) << " " << xyz(1) << " " << xyz(2);
402  TrajPos.push_back(xyz);
403  TVector3 dir = TrajPos[1] - TrajPos[0];
404  dir = dir.Unit();
405  TrajDir.push_back(dir);
406  TrajDir.push_back(dir);
407 
408  if (prt)
409  mf::LogVerbatim("TTA") << ">>>> Short track (" << TrajPos[0](0) << ", " << TrajPos[0](1)
410  << ", " << TrajPos[0](2) << ") to (" << TrajPos[1](0) << ", "
411  << TrajPos[1](1) << ", " << TrajPos[1](2) << ")";
412  }
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:216
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
geo::WireReadoutGeom::WireIDsIntersect
bool WireIDsIntersect(WireID const &wid1, WireID const &wid2, Point_t &intersection) const
Computes the intersection between two wires.
Definition: WireReadoutGeom.cxx:344
dir
TDirectory * dir
Definition: macro.C:5
trkf::TrackTrajectoryAlg::wireReadoutGeom
geo::WireReadoutGeom const * wireReadoutGeom
Definition: TrackTrajectoryAlg.h:35
geo::WireIDIntersection::invalid
static constexpr WireIDIntersection invalid()
Definition: geo_types.h:594
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 23 of file TrackTrajectoryAlg.cxx.

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

Referenced by trkf::CCTrackMaker::StoreTrack().

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

Member Data Documentation

double trkf::TrackTrajectoryAlg::fHitWidthFactor
private

Definition at line 44 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

unsigned short trkf::TrackTrajectoryAlg::fMaxTrajPoints
private

Definition at line 43 of file TrackTrajectoryAlg.h.

Referenced by TrackTrajectory().

TrackLineFitAlg trkf::TrackTrajectoryAlg::fTrackLineFitAlg
private

Definition at line 46 of file TrackTrajectoryAlg.h.

Referenced by TrackTrajectory().

double trkf::TrackTrajectoryAlg::maxX
private

Definition at line 39 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

unsigned short trkf::TrackTrajectoryAlg::maxXPln
private

Definition at line 40 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

double trkf::TrackTrajectoryAlg::minX
private

Definition at line 37 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

unsigned short trkf::TrackTrajectoryAlg::minXPln
private

Definition at line 38 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

bool trkf::TrackTrajectoryAlg::prt
private

Definition at line 41 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

geo::WireReadoutGeom const* trkf::TrackTrajectoryAlg::wireReadoutGeom = &art::ServiceHandle<geo::WireReadout>()->Get()
private

Definition at line 35 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory().

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