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

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  }
mf::LogVerbatim
MaybeLogger_< ELseverityLevel::ELsev_info, true > LogVerbatim
Definition: MessageLogger.h:216
reco_momentum_tuples.float
float
Definition: reco_momentum_tuples.py:12
dir
string dir
Definition: INukeNucleonCorr.cxx:61
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
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:1657
trkf::TrackTrajectoryAlg::geom
art::ServiceHandle< geo::Geometry const > geom
Definition: TrackTrajectoryAlg.h:36
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.

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
ValidateOpDetReco.end
end
while True: pbar.update(maxval-len(onlies[E][S])) #print iS, "/", len(onlies[E][S]) found = False for...
Definition: ValidateOpDetReco.py:548
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
dir
string dir
Definition: INukeNucleonCorr.cxx:61
setprecision
Q_EXPORT QTSManip setprecision(int p)
Definition: qtextstream.h:343
util::size
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
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
reverse
static unsigned int reverse(QString &chars, unsigned char *level, unsigned int a, unsigned int b)
Definition: qstring.cpp:11649
util::begin
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:72

Member Data Documentation

double trkf::TrackTrajectoryAlg::fHitWidthFactor
private

Definition at line 45 of file TrackTrajectoryAlg.h.

unsigned short trkf::TrackTrajectoryAlg::fMaxTrajPoints
private

Definition at line 44 of file TrackTrajectoryAlg.h.

TrackLineFitAlg trkf::TrackTrajectoryAlg::fTrackLineFitAlg
private

Definition at line 48 of file TrackTrajectoryAlg.h.

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

Definition at line 36 of file TrackTrajectoryAlg.h.

double trkf::TrackTrajectoryAlg::maxX
private

Definition at line 40 of file TrackTrajectoryAlg.h.

unsigned short trkf::TrackTrajectoryAlg::maxXPln
private

Definition at line 41 of file TrackTrajectoryAlg.h.

double trkf::TrackTrajectoryAlg::minX
private

Definition at line 38 of file TrackTrajectoryAlg.h.

unsigned short trkf::TrackTrajectoryAlg::minXPln
private

Definition at line 39 of file TrackTrajectoryAlg.h.

bool trkf::TrackTrajectoryAlg::prt
private

Definition at line 42 of file TrackTrajectoryAlg.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.11