All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TrackTrajectoryAlg.cxx
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////
2 ///
3 /// TrackTrajectoryAlg class
4 ///
5 /// Bruce Baller, baller@fnal.gov
6 ///
7 /// Algorithm fitting a 3D trajectory through a set of hits
8 ///
9 ////////////////////////////////////////////////////////////////////////
10 
11 #include <algorithm>
12 #include <iomanip>
13 
14 #include "TVector3.h"
15 #include "messagefacility/MessageLogger/MessageLogger.h"
16 
17 #include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
18 #include "larreco/RecoAlg/TrackTrajectoryAlg.h"
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, std::vector<TVector3>& TrajDir)
27  {
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 
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) {
67  minLen = trkX[ipl].size();
68  }
69  if(trkX[ipl].size() > maxLen) {
70  maxLen = trkX[ipl].size();
71  aPlane = ipl;
72  }
73  } // ipl
74  if(prt) mf::LogVerbatim("TTA")<<"trkX sizes "<<trkX[0].size()<<" "<<trkX[1].size()<<" "<<trkX[2].size()<<" "<<" 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) mf::LogVerbatim("TTA")<<"minX "<<minX<<" in plane "<<minXPln<<" maxX "<<maxX<<" in plane "<<maxXPln;
135 
136  // estimate the number of trajectory points we will want based on the delta T and the
137  // hit errors
138  double aveHitXErr = 0;
139  unsigned short nHit = 0;
140  for(ipl = 0; ipl < 3; ++ipl) {
141  for(iht = 0; iht < trkXErr[ipl].size(); ++iht) {
142  aveHitXErr += trkXErr[ipl][iht];
143  ++nHit;
144  }
145  } // ipl
146  aveHitXErr /= (double)nHit;
147 
148  unsigned short npt = (maxX - minX) / (1 * aveHitXErr);
149  if(npt < 2) npt = 2;
150  if(npt > maxLen) npt = maxLen;
151  if(npt > fMaxTrajPoints) npt = fMaxTrajPoints;
152  if(prt) mf::LogVerbatim("TTA")<<" aveHitXErr "<<aveHitXErr<<" number of traj points ";
153 
154  double maxBinX = (maxX - minX) / (double)(npt - 1);
155  double binX = maxBinX;
156  double xOrigin = minX;
157 
158  TVector3 xyz, dir;
159  // working vectors passed to TrackLineFitAlg
160  std::vector<geo::WireID> hitWID;
161  std::vector<double> hitX;
162  std::vector<double> hitXErr;
163  float ChiDOF;
164 
165  // make a short track trajectory (end points only) to use in case an error occurs
166  std::vector<TVector3> STPos;
167  std::vector<TVector3> STDir;
168  ShortTrackTrajectory(trkWID, trkX, trkXErr, STPos, STDir);
169 
170  if(STPos.size() != 2 || STDir.size() != STPos.size()) {
171  TrajPos.clear();
172  TrajDir.clear();
173  if(posX) {
174  for(ipl = 0; ipl < 3; ++ipl) {
175  if(trkX[ipl].size() == 0) continue;
176  std::reverse(trkWID[ipl].begin(), trkWID[ipl].end());
177  std::reverse(trkX[ipl].begin(), trkX[ipl].end());
178  std::reverse(trkXErr[ipl].begin(), trkXErr[ipl].end());
179  } // ipl
180  } // posX
181  return;
182  } // bad STPos, STDir
183 
184  if(prt) {
185  mf::LogVerbatim("TTA")<<"STPos";
186  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);
187  }
188 
189  if(maxLen < 4 || npt < 2) {
190  TrajPos = STPos;
191  TrajDir = STDir;
192  if(!posX) {
193  // reverse everything
194  std::reverse(TrajPos.begin(), TrajPos.end());
195  std::reverse(TrajDir.begin(), TrajDir.end());
196  for(ipl = 0; ipl < 3; ++ipl) {
197  if(trkX[ipl].size() == 0) continue;
198  std::reverse(trkWID[ipl].begin(), trkWID[ipl].end());
199  std::reverse(trkX[ipl].begin(), trkX[ipl].end());
200  std::reverse(trkXErr[ipl].begin(), trkXErr[ipl].end());
201  } // ipl
202  }
203  return;
204  } // maxLen < 4
205 
206  // Start index of hits to include in the next fit
207  std::array<unsigned short, 3> hStart;
208  for(ipl = 0; ipl < 3; ++ipl) hStart[ipl] = 0;
209 
210  bool gotLastPoint = false;
211  for(unsigned short ipt = 0; ipt < npt + 1; ++ipt) {
212  hitWID.clear();
213  hitX.clear();
214  hitXErr.clear();
215  for(ipl = 0; ipl < 3; ++ipl) {
216  for(iht = hStart[ipl]; iht < trkX[ipl].size(); ++iht) {
217  if(trkX[ipl][iht] < xOrigin - binX) continue;
218  if(trkX[ipl][iht] > xOrigin + binX) break;
219  hitWID.push_back(trkWID[ipl][iht]);
220  hitX.push_back(trkX[ipl][iht]);
221  hitXErr.push_back(trkXErr[ipl][iht]);
222  hStart[ipl] =iht;
223  } // iht
224  } // ipl
225  if(prt) mf::LogVerbatim("TTA")<<"ipt "<<ipt<<" xOrigin "<<xOrigin<<" binX "<<binX<<" hitX size "<<hitX.size();
226  if(hitX.size() > 3) {
227  fTrackLineFitAlg.TrkLineFit(hitWID, hitX, hitXErr, xOrigin, xyz, dir, ChiDOF);
228  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();
229  } else if(ipt == 0 && STPos.size() == 2) {
230  // failure on the first traj point. Use STPos
231  xyz = STPos[0];
232  dir = STDir[0];
233  }
234  if(prt && hitX.size() < 4) mf::LogVerbatim("TTA")<<"\n";
235  if(xOrigin >= maxX) break;
236  // tweak xOrigin if we are close to the end
237  if(maxX - xOrigin < binX) {
238  xOrigin = maxX;
239  } else {
240  xOrigin += binX;
241  }
242  if(ChiDOF < 0 || ChiDOF > 100) continue;
243  TrajPos.push_back(xyz);
244  TrajDir.push_back(dir);
245  if(ipt == npt) gotLastPoint = true;
246  } // ipt
247  if(prt) {
248  mf::LogVerbatim("TTA")<<"gotLastPoint "<<gotLastPoint<<" TTA Traj \n";
249  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);
250  }
251 
252  if(!gotLastPoint && STPos.size() == 2) {
253  // failure on the last point. Use STPos, etc
254  TrajPos.push_back(STPos[1]);
255  TrajDir.push_back(STDir[1]);
256  }
257 
258  if(TrajPos.size() < 2) {
259  TrajPos = STPos;
260  TrajDir = STDir;
261  }
262 
263  if(!posX) {
264  // reverse everything
265  std::reverse(TrajPos.begin(), TrajPos.end());
266  std::reverse(TrajDir.begin(), TrajDir.end());
267  for(ipl = 0; ipl < 3; ++ipl) {
268  if(trkX[ipl].size() == 0) continue;
269  std::reverse(trkWID[ipl].begin(), trkWID[ipl].end());
270  std::reverse(trkX[ipl].begin(), trkX[ipl].end());
271  std::reverse(trkXErr[ipl].begin(), trkXErr[ipl].end());
272  } // ipl
273 
274  } // !posX
275  if(prt) {
276  mf::LogVerbatim("TTA")<<"TTA Traj2\n";
277  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);
278  }
279 
280  } // TrackTrajectoryAlg
281 
282 ////////////////////////////////////////////////////////////////////////////////////////////////////////////
283  void TrackTrajectoryAlg::ShortTrackTrajectory(std::array<std::vector<geo::WireID>,3> trkWID,
284  std::array<std::vector<double>,3> trkX,
285  std::array<std::vector<double>,3> trkXErr,
286  std::vector<TVector3>& TrajPos, std::vector<TVector3>& TrajDir)
287  {
288  // Make a short track trajectory composed of a space point at each end.
289 
290  std::vector<unsigned short> usePlns;
291  double y, z, endX, endY, endZ;
292  unsigned short ipl, iht, nend, jpl, iPlane, jPlane;
293  unsigned int iWire, jWire;
294 
295  // do the min X end first
296  float xCut;
297  for(unsigned short nit = 0; nit < 5; ++nit) {
298  xCut = minX + fHitWidthFactor * trkXErr[minXPln][0];
299  for(ipl = 0; ipl < 3; ++ipl) {
300  if(trkX[ipl].size() == 0) continue;
301  if(trkX[ipl][0] < xCut) usePlns.push_back(ipl);
302  } // ipl
303  if(usePlns.size() >= 2) break;
304  fHitWidthFactor += 5;
305  }
306  // Not enough information to find a space point
307  if(prt) mf::LogVerbatim("TTA")<<"ShortTrack minX end "<<xCut<<" usePlns size "<<usePlns.size();
308  if(usePlns.size() < 2) return;
309  endY = 0; endZ = 0; nend = 0;
310  iht = 0;
311  ipl = usePlns[0];
312  endX = trkX[ipl][iht];
313  iPlane = trkWID[ipl][iht].Plane;
314  iWire = trkWID[ipl][iht].Wire;
315  unsigned int cstat = trkWID[ipl][iht].Cryostat;
316  unsigned int tpc = trkWID[ipl][iht].TPC;
317  for(unsigned short jj = 1; jj < usePlns.size(); ++jj) {
318  jpl = usePlns[jj];
319  endX += trkX[jpl][iht];
320  jPlane = trkWID[jpl][iht].Plane;
321  jWire = trkWID[jpl][iht].Wire;
322  geom->IntersectionPoint(iWire, jWire, iPlane, jPlane, cstat, tpc, y, z);
323  endY += y;
324  endZ += z;
325  ++nend;
326  } // ii
327  // nend is the number of y,z values. There are (nend + 1) X values
328  TVector3 xyz;
329  xyz(0) = endX / (float)(nend + 1);
330  xyz(1) = endY / (float)nend;
331  xyz(2) = endZ / (float)nend;
332  if(prt) mf::LogVerbatim("TTA")<<" xyz "<<xyz(0)<<" "<<xyz(1)<<" "<<xyz(2);
333  TrajPos.push_back(xyz);
334 
335  // do the same for the other end
336  fHitWidthFactor = 5;
337 // xCut = maxX - fHitWidthFactor * trkXErr[maxXPln][0];
338  usePlns.clear();
339  for(unsigned short nit = 0; nit < 5; ++nit) {
340  xCut = maxX - fHitWidthFactor * trkXErr[maxXPln][0];
341  for(ipl = 0; ipl < 3; ++ipl) {
342  if(trkX[ipl].size() == 0) continue;
343  iht = trkX[ipl].size() - 1;
344  if(trkX[ipl][iht] > xCut) usePlns.push_back(ipl);
345  } // ipl
346  if(usePlns.size() >= 2) break;
347  fHitWidthFactor += 5;
348  }
349  // Not enough information to find a space point
350  if(prt) mf::LogVerbatim("TTA")<<"ShortTrack maxX end "<<xCut<<" usePlns size "<<usePlns.size();
351  if(usePlns.size() < 2) {
352  TrajPos.clear();
353  TrajDir.clear();
354  return;
355  }
356  endY = 0; endZ = 0; nend = 0;
357  ipl = usePlns[0];
358  iht = trkX[ipl].size() - 1;
359  endX = trkX[ipl][iht];
360  iPlane = trkWID[ipl][iht].Plane;
361  iWire = trkWID[ipl][iht].Wire;
362  for(unsigned short jj = 1; jj < usePlns.size(); ++jj) {
363  jpl = usePlns[jj];
364  iht = trkX[jpl].size() - 1;
365  endX += trkX[jpl][iht];
366  jPlane = trkWID[jpl][iht].Plane;
367  jWire = trkWID[jpl][iht].Wire;
368  geom->IntersectionPoint(iWire, jWire, iPlane, jPlane, cstat, tpc, y, z);
369  endY += y;
370  endZ += z;
371  ++nend;
372  } // ii
373  // nend is the number of y,z values. There are nend + 1 X values
374  xyz(0) = endX / (float)(nend + 1);
375  xyz(1) = endY / (float)nend;
376  xyz(2) = endZ / (float)nend;
377  if(prt) mf::LogVerbatim("TTA")<<" xyz "<<xyz(0)<<" "<<xyz(1)<<" "<<xyz(2);
378  TrajPos.push_back(xyz);
379  TVector3 dir = TrajPos[1] - TrajPos[0];
380  dir = dir.Unit();
381  TrajDir.push_back(dir);
382  TrajDir.push_back(dir);
383 
384  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)<<")";
385  }
386 
387 
388 } // trkf
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
z
process_name opflash particleana ie ie ie z
Definition: run_opflash_electron.fcl:92
trkf::TrackTrajectoryAlg::fTrackLineFitAlg
TrackLineFitAlg fTrackLineFitAlg
Definition: TrackTrajectoryAlg.h:48
icarus::ns::util::size
std::size_t size(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:561
y
process_name opflash particleana ie ie y
Definition: run_opflash_electron.fcl:91
array
then echo ***************************************echo array
Definition: find_fhicl.sh:28
icarus::ns::util::end
auto end(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:585
trkf::TrackTrajectoryAlg::TrackTrajectory
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: TrackTrajectoryAlg.cxx:23
geo_types.h
Definition of data types for geometry description.
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:283
dropbox.dir
tuple dir
Definition: dropbox.py:28
icarus::ns::util::begin
auto begin(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:573
trkf::TrackTrajectoryAlg::geom
art::ServiceHandle< geo::Geometry const > geom
Definition: TrackTrajectoryAlg.h:36
Generated by   doxygen 1.8.5