LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
TrackTrajectoryAlg.cxx
Go to the documentation of this file.
1 
11 #include <algorithm>
12 #include <iomanip>
13 
14 #include "TVector3.h"
16 
19 
20 namespace trkf {
21 
22  //------------------------------------------------------------------------------
23  void TrackTrajectoryAlg::TrackTrajectory(std::array<std::vector<geo::WireID>, 3> trkWID,
24  std::array<std::vector<double>, 3> trkX,
25  std::array<std::vector<double>, 3> trkXErr,
26  std::vector<TVector3>& TrajPos,
27  std::vector<TVector3>& TrajDir)
28  {
29 
30  // Make a track trajectory (position, direction) and return it in the TrajPos
31  // and TrajDir vectors. The track hits are received as 3 vectors (one vector per wire plane)
32  // of Wire ID's, hit X and X position errors. The X position errors are used to 1) determine
33  // the significance of the X difference between the beginning and end of the track trajectory
34  // 2) determine the number of trajectory points and 3) weight the track line fit in TrackLineFitAlg.
35  // This code assumes that hits at each end (e.g. trkXW[Plane][0]) of the vectors define the end
36  // points of the trajectory. The ordering of planes in the array is irrelevant since the
37  // plane number is extracted from the WireID. This algorithm will return with a failed condition
38  // (TrajPos, TrajDir size = 0) if there are fewer than 2 planes with hits at each end that are less
39  // than fHitWidthFactor * trkXErr apart. Valid and invalid conditions are shown schematically below
40  // where a . represents hits that are separated by X values > fHitWidthFactor * trkXErr
41  //
42  // minX maxX maxX minX minX maxX
43  // Pln0 .................... Pln0 .................... Pln0 ...................
44  // Pln1 ............... Pln1 .................... Pln1
45  // Pln2 .................... Pln2 ................ Pln2 ...................
46  // VALID VALID VALID - no hits in one plane is OK
47  //
48  // minX maxX
49  // Pln0 .................
50  // Pln1 ...............
51  // Pln2 ....................
52  // NOT VALID - Only one plane has a hit at MaxX
53 
54  TrajPos.clear();
55  TrajDir.clear();
56 
57  prt = false;
58 
59  unsigned short minLen = 9999;
60  unsigned short maxLen = 0;
61  unsigned short nPlnsWithHits = 0;
62  unsigned short ipl, aPlane = 3;
63  for (ipl = 0; ipl < 3; ++ipl) {
64  if (trkX[ipl].size() == 0) continue;
65  ++nPlnsWithHits;
66  if (trkX[ipl].size() < minLen) { minLen = trkX[ipl].size(); }
67  if (trkX[ipl].size() > maxLen) {
68  maxLen = trkX[ipl].size();
69  aPlane = ipl;
70  }
71  } // ipl
72  if (prt)
73  mf::LogVerbatim("TTA") << "trkX sizes " << trkX[0].size() << " " << trkX[1].size() << " "
74  << trkX[2].size() << " "
75  << " nPlnsWithHits " << nPlnsWithHits;
76  if (nPlnsWithHits < 2) return;
77  if (aPlane > 2) return;
78 
79  fMaxTrajPoints = 100;
80  fHitWidthFactor = 5.;
81 
82  unsigned short iht;
83 
84  // reverse the order of the hit arrays if necessary so that the minimum X end is at trk[][0] end.
85  // We will use posX to reverse the trajectory later if necessary
86  bool posX = true;
87  iht = trkX[aPlane].size() - 1;
88  if (trkX[aPlane][0] > trkX[aPlane][iht]) {
89  posX = false;
90  for (ipl = 0; ipl < 3; ++ipl) {
91  if (trkX[ipl].size() == 0) continue;
92  std::reverse(trkWID[ipl].begin(), trkWID[ipl].end());
93  std::reverse(trkX[ipl].begin(), trkX[ipl].end());
94  std::reverse(trkXErr[ipl].begin(), trkXErr[ipl].end());
95  } // ipl
96  if (prt) mf::LogVerbatim("TTA") << "Swapped order";
97  } // posX check
98 
99  if (prt) {
100  mf::LogVerbatim myprt("TTA");
101  myprt << "TrkXW end0";
102  for (ipl = 0; ipl < 3; ++ipl) {
103  if (trkX[ipl].size() == 0) continue;
104  myprt << " " << trkX[ipl][0];
105  }
106  myprt << "\n";
107  myprt << "TrkXW end1";
108  for (ipl = 0; ipl < 3; ++ipl) {
109  if (trkX[ipl].size() == 0) continue;
110  iht = trkX[ipl].size() - 1;
111  myprt << " " << trkX[ipl][iht];
112  }
113  myprt << "\n";
114  }
115 
116  // find the min/max X
117  minX = 1E6;
118  minXPln = 0;
119  maxX = -1E6;
120  maxXPln = 0;
121 
122  for (unsigned short ipl = 0; ipl < 3; ++ipl) {
123  if (trkX[ipl].size() == 0) continue;
124  if (trkX[ipl][0] < minX) {
125  minX = trkX[ipl][0];
126  minXPln = ipl;
127  }
128  iht = trkX[ipl].size() - 1;
129  if (trkX[ipl][iht] > maxX) {
130  maxX = trkX[ipl][iht];
131  maxXPln = ipl;
132  }
133  } // ipl
134 
135  if (prt)
136  mf::LogVerbatim("TTA") << "minX " << minX << " in plane " << minXPln << " maxX " << maxX
137  << " in plane " << maxXPln;
138 
139  // estimate the number of trajectory points we will want based on the delta T and the
140  // hit errors
141  double aveHitXErr = 0;
142  unsigned short nHit = 0;
143  for (ipl = 0; ipl < 3; ++ipl) {
144  for (iht = 0; iht < trkXErr[ipl].size(); ++iht) {
145  aveHitXErr += trkXErr[ipl][iht];
146  ++nHit;
147  }
148  } // ipl
149  aveHitXErr /= (double)nHit;
150 
151  unsigned short npt = (maxX - minX) / (1 * aveHitXErr);
152  if (npt < 2) npt = 2;
153  if (npt > maxLen) npt = maxLen;
154  if (npt > fMaxTrajPoints) npt = fMaxTrajPoints;
155  if (prt) mf::LogVerbatim("TTA") << " aveHitXErr " << aveHitXErr << " number of traj points ";
156 
157  double maxBinX = (maxX - minX) / (double)(npt - 1);
158  double binX = maxBinX;
159  double xOrigin = minX;
160 
161  TVector3 xyz, dir;
162  // working vectors passed to TrackLineFitAlg
163  std::vector<geo::WireID> hitWID;
164  std::vector<double> hitX;
165  std::vector<double> hitXErr;
166  float ChiDOF;
167 
168  // make a short track trajectory (end points only) to use in case an error occurs
169  std::vector<TVector3> STPos;
170  std::vector<TVector3> STDir;
171  ShortTrackTrajectory(trkWID, trkX, trkXErr, STPos, STDir);
172 
173  if (STPos.size() != 2 || STDir.size() != STPos.size()) {
174  TrajPos.clear();
175  TrajDir.clear();
176  if (posX) {
177  for (ipl = 0; ipl < 3; ++ipl) {
178  if (trkX[ipl].size() == 0) continue;
179  std::reverse(trkWID[ipl].begin(), trkWID[ipl].end());
180  std::reverse(trkX[ipl].begin(), trkX[ipl].end());
181  std::reverse(trkXErr[ipl].begin(), trkXErr[ipl].end());
182  } // ipl
183  } // posX
184  return;
185  } // bad STPos, STDir
186 
187  if (prt) {
188  mf::LogVerbatim("TTA") << "STPos";
189  for (unsigned short ii = 0; ii < STPos.size(); ++ii)
190  mf::LogVerbatim("TTA") << ii << " " << std::fixed << std::setprecision(1) << STPos[ii](0)
191  << " " << STPos[ii](1) << " " << STPos[ii](2);
192  }
193 
194  if (maxLen < 4 || npt < 2) {
195  TrajPos = STPos;
196  TrajDir = STDir;
197  if (!posX) {
198  // reverse everything
199  std::reverse(TrajPos.begin(), TrajPos.end());
200  std::reverse(TrajDir.begin(), TrajDir.end());
201  for (ipl = 0; ipl < 3; ++ipl) {
202  if (trkX[ipl].size() == 0) continue;
203  std::reverse(trkWID[ipl].begin(), trkWID[ipl].end());
204  std::reverse(trkX[ipl].begin(), trkX[ipl].end());
205  std::reverse(trkXErr[ipl].begin(), trkXErr[ipl].end());
206  } // ipl
207  }
208  return;
209  } // maxLen < 4
210 
211  // Start index of hits to include in the next fit
212  std::array<unsigned short, 3> hStart;
213  for (ipl = 0; ipl < 3; ++ipl)
214  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)
232  mf::LogVerbatim("TTA") << "ipt " << ipt << " xOrigin " << xOrigin << " binX " << binX
233  << " hitX size " << hitX.size();
234  if (hitX.size() > 3) {
235  fTrackLineFitAlg.TrkLineFit(hitWID, hitX, hitXErr, xOrigin, xyz, dir, ChiDOF);
236  if (prt)
237  mf::LogVerbatim("TTA") << " xyz " << xyz(0) << " " << xyz(1) << " " << xyz(2) << " dir "
238  << dir(0) << " " << dir(1) << " " << dir(2) << " ChiDOF " << ChiDOF
239  << " hitX size " << hitX.size();
240  }
241  else if (ipt == 0 && STPos.size() == 2) {
242  // failure on the first traj point. Use STPos
243  xyz = STPos[0];
244  dir = STDir[0];
245  }
246  if (prt && hitX.size() < 4) mf::LogVerbatim("TTA") << "\n";
247  if (xOrigin >= maxX) break;
248  // tweak xOrigin if we are close to the end
249  if (maxX - xOrigin < binX) { xOrigin = maxX; }
250  else {
251  xOrigin += binX;
252  }
253  if (ChiDOF < 0 || ChiDOF > 100) continue;
254  TrajPos.push_back(xyz);
255  TrajDir.push_back(dir);
256  if (ipt == npt) gotLastPoint = true;
257  } // ipt
258  if (prt) {
259  mf::LogVerbatim("TTA") << "gotLastPoint " << gotLastPoint << " TTA Traj \n";
260  for (unsigned short ii = 0; ii < TrajPos.size(); ++ii)
261  mf::LogVerbatim("TTA") << ii << " " << std::fixed << std::setprecision(1) << TrajPos[ii](0)
262  << " " << TrajPos[ii](1) << " " << TrajPos[ii](2);
263  }
264 
265  if (!gotLastPoint && STPos.size() == 2) {
266  // failure on the last point. Use STPos, etc
267  TrajPos.push_back(STPos[1]);
268  TrajDir.push_back(STDir[1]);
269  }
270 
271  if (TrajPos.size() < 2) {
272  TrajPos = STPos;
273  TrajDir = STDir;
274  }
275 
276  if (!posX) {
277  // reverse everything
278  std::reverse(TrajPos.begin(), TrajPos.end());
279  std::reverse(TrajDir.begin(), TrajDir.end());
280  for (ipl = 0; ipl < 3; ++ipl) {
281  if (trkX[ipl].size() == 0) continue;
282  std::reverse(trkWID[ipl].begin(), trkWID[ipl].end());
283  std::reverse(trkX[ipl].begin(), trkX[ipl].end());
284  std::reverse(trkXErr[ipl].begin(), trkXErr[ipl].end());
285  } // ipl
286 
287  } // !posX
288  if (prt) {
289  mf::LogVerbatim("TTA") << "TTA Traj2\n";
290  for (unsigned short ii = 0; ii < TrajPos.size(); ++ii)
291  mf::LogVerbatim("TTA") << ii << " " << std::fixed << std::setprecision(1) << TrajPos[ii](0)
292  << " " << TrajPos[ii](1) << " " << TrajPos[ii](2);
293  }
294 
295  } // TrackTrajectoryAlg
296 
298  void TrackTrajectoryAlg::ShortTrackTrajectory(std::array<std::vector<geo::WireID>, 3> trkWID,
299  std::array<std::vector<double>, 3> trkX,
300  std::array<std::vector<double>, 3> trkXErr,
301  std::vector<TVector3>& TrajPos,
302  std::vector<TVector3>& TrajDir)
303  {
304  // Make a short track trajectory composed of a space point at each end.
305 
306  std::vector<unsigned short> usePlns;
307  double y, z, endX, endY, endZ;
308  unsigned short ipl, iht, nend, jpl, iPlane, jPlane;
309  unsigned int iWire, jWire;
310 
311  // do the min X end first
312  float xCut;
313  for (unsigned short nit = 0; nit < 5; ++nit) {
314  xCut = minX + fHitWidthFactor * trkXErr[minXPln][0];
315  for (ipl = 0; ipl < 3; ++ipl) {
316  if (trkX[ipl].size() == 0) continue;
317  if (trkX[ipl][0] < xCut) usePlns.push_back(ipl);
318  } // ipl
319  if (usePlns.size() >= 2) break;
320  fHitWidthFactor += 5;
321  }
322  // Not enough information to find a space point
323  if (prt)
324  mf::LogVerbatim("TTA") << "ShortTrack minX end " << xCut << " usePlns size "
325  << usePlns.size();
326  if (usePlns.size() < 2) return;
327  endY = 0;
328  endZ = 0;
329  nend = 0;
330  iht = 0;
331  ipl = usePlns[0];
332  endX = trkX[ipl][iht];
333  iPlane = trkWID[ipl][iht].Plane;
334  iWire = trkWID[ipl][iht].Wire;
335  unsigned int cstat = trkWID[ipl][iht].Cryostat;
336  unsigned int tpc = trkWID[ipl][iht].TPC;
337  for (unsigned short jj = 1; jj < usePlns.size(); ++jj) {
338  jpl = usePlns[jj];
339  endX += trkX[jpl][iht];
340  jPlane = trkWID[jpl][iht].Plane;
341  jWire = trkWID[jpl][iht].Wire;
343  geo::WireID{cstat, tpc, iPlane, iWire}, geo::WireID{cstat, tpc, jPlane, jWire}, y, z);
344  endY += y;
345  endZ += z;
346  ++nend;
347  } // ii
348  // nend is the number of y,z values. There are (nend + 1) X values
349  TVector3 xyz;
350  xyz(0) = endX / (float)(nend + 1);
351  xyz(1) = endY / (float)nend;
352  xyz(2) = endZ / (float)nend;
353  if (prt) mf::LogVerbatim("TTA") << " xyz " << xyz(0) << " " << xyz(1) << " " << xyz(2);
354  TrajPos.push_back(xyz);
355 
356  // do the same for the other end
357  fHitWidthFactor = 5;
358  // xCut = maxX - fHitWidthFactor * trkXErr[maxXPln][0];
359  usePlns.clear();
360  for (unsigned short nit = 0; nit < 5; ++nit) {
361  xCut = maxX - fHitWidthFactor * trkXErr[maxXPln][0];
362  for (ipl = 0; ipl < 3; ++ipl) {
363  if (trkX[ipl].size() == 0) continue;
364  iht = trkX[ipl].size() - 1;
365  if (trkX[ipl][iht] > xCut) usePlns.push_back(ipl);
366  } // ipl
367  if (usePlns.size() >= 2) break;
368  fHitWidthFactor += 5;
369  }
370  // Not enough information to find a space point
371  if (prt)
372  mf::LogVerbatim("TTA") << "ShortTrack maxX end " << xCut << " usePlns size "
373  << usePlns.size();
374  if (usePlns.size() < 2) {
375  TrajPos.clear();
376  TrajDir.clear();
377  return;
378  }
379  endY = 0;
380  endZ = 0;
381  nend = 0;
382  ipl = usePlns[0];
383  iht = trkX[ipl].size() - 1;
384  endX = trkX[ipl][iht];
385  iPlane = trkWID[ipl][iht].Plane;
386  iWire = trkWID[ipl][iht].Wire;
387  for (unsigned short jj = 1; jj < usePlns.size(); ++jj) {
388  jpl = usePlns[jj];
389  iht = trkX[jpl].size() - 1;
390  endX += trkX[jpl][iht];
391  jPlane = trkWID[jpl][iht].Plane;
392  jWire = trkWID[jpl][iht].Wire;
394  geo::WireID{cstat, tpc, iPlane, iWire}, geo::WireID{cstat, tpc, jPlane, jWire}, y, z);
395  endY += y;
396  endZ += z;
397  ++nend;
398  } // ii
399  // nend is the number of y,z values. There are nend + 1 X values
400  xyz(0) = endX / (float)(nend + 1);
401  xyz(1) = endY / (float)nend;
402  xyz(2) = endZ / (float)nend;
403  if (prt) mf::LogVerbatim("TTA") << " xyz " << xyz(0) << " " << xyz(1) << " " << xyz(2);
404  TrajPos.push_back(xyz);
405  TVector3 dir = TrajPos[1] - TrajPos[0];
406  dir = dir.Unit();
407  TrajDir.push_back(dir);
408  TrajDir.push_back(dir);
409 
410  if (prt)
411  mf::LogVerbatim("TTA") << ">>>> Short track (" << TrajPos[0](0) << ", " << TrajPos[0](1)
412  << ", " << TrajPos[0](2) << ") to (" << TrajPos[1](0) << ", "
413  << TrajPos[1](1) << ", " << TrajPos[1](2) << ")";
414  }
415 
416 } // trkf
void TrkLineFit(std::vector< geo::WireID > &hitWID, std::vector< double > &hitX, std::vector< double > &hitXErr, double XOrigin, TVector3 &Pos, TVector3 &Dir, float &ChiDOF) const
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
TrackLineFitAlg fTrackLineFitAlg
Float_t y
Definition: compare.C:6
Double_t z
Definition: plot.C:276
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
auto array(Array const &a)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:250
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)
Definition of data types for geometry description.
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)
TDirectory * dir
Definition: macro.C:5
art::ServiceHandle< geo::Geometry const > geom
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
bool IntersectionPoint(WireID const &wid1, WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
Float_t posX
Definition: plotDend.C:20