LArSoft  v09_90_00
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

art::ServiceHandle< geo::Geometry const > 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.

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

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

Referenced by TrackTrajectory().

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;
342  geom->IntersectionPoint(
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;
393  geom->IntersectionPoint(
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  }
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:216
y
Float_t y
Definition: compare.C:6
z
Double_t z
Definition: plot.C:276
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:101
dir
TDirectory * dir
Definition: macro.C:5
trkf::TrackTrajectoryAlg::geom
art::ServiceHandle< geo::Geometry const > geom
Definition: TrackTrajectoryAlg.h:37
geo::GeometryCore::IntersectionPoint
bool IntersectionPoint(WireID const &wid1, WireID const &wid2, double &y, double &z) const
Returns the intersection point of two wires.
Definition: GeometryCore.cxx:1451
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 
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
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:28
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:216
trkf::TrackTrajectoryAlg::fTrackLineFitAlg
TrackLineFitAlg fTrackLineFitAlg
Definition: TrackTrajectoryAlg.h:48
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:298
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 46 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

unsigned short trkf::TrackTrajectoryAlg::fMaxTrajPoints
private

Definition at line 45 of file TrackTrajectoryAlg.h.

Referenced by TrackTrajectory().

TrackLineFitAlg trkf::TrackTrajectoryAlg::fTrackLineFitAlg
private

Definition at line 48 of file TrackTrajectoryAlg.h.

Referenced by TrackTrajectory().

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

Definition at line 37 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory().

double trkf::TrackTrajectoryAlg::maxX
private

Definition at line 41 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

unsigned short trkf::TrackTrajectoryAlg::maxXPln
private

Definition at line 42 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

double trkf::TrackTrajectoryAlg::minX
private

Definition at line 39 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

unsigned short trkf::TrackTrajectoryAlg::minXPln
private

Definition at line 40 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

bool trkf::TrackTrajectoryAlg::prt
private

Definition at line 43 of file TrackTrajectoryAlg.h.

Referenced by ShortTrackTrajectory(), and TrackTrajectory().

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