ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHTpcTrackFollower.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHTpcTrackFollower.cc
1 
8 #include "PHTpcTrackFollower.h"
9 #include "PHTpcConst.h"
10 
11 #include "Fitter.h" // for Fitter
12 #include "PHTpcLookup.h" // for PHTpcLookup
13 #include "SpacepointMeasurement2.h"
14 #include "Track.h"
15 #include "externals/kdfinder.hpp" // for TrackCandidate, Helix, TVector
16 
17 #include <trackbase/TrkrDefs.h> // for cluskey
18 
19 #include <phool/PHLog.h>
20 
21 #include <GenFit/FieldManager.h>
22 #include <GenFit/KalmanFitStatus.h>
23 #include <GenFit/KalmanFitterInfo.h>
24 #include <GenFit/MeasuredStateOnPlane.h>
25 #include <GenFit/RKTrackRep.h>
26 #include <GenFit/Track.h>
27 #include <GenFit/TrackPoint.h>
28 
29 #include <TMatrixDSymfwd.h> // for TMatrixDSym
30 #include <TMatrixTSym.h> // for TMatrixTSym
31 #include <TVector3.h> // for TVector3
32 
33 #include <log4cpp/CategoryStream.hh> // for CategoryStream
34 
35 #include <algorithm>
36 #include <cmath> // for isnan, sqrt
37 #include <cstddef> // for size_t
38 #include <cstdint> // for uint64_t, int64_t
39 #include <limits> // for numeric_limits
40 #include <unordered_set>
41 
42 class PHField;
44 
45 namespace PHGenFit
46 {
47  class Measurement;
48 }
49 namespace genfit
50 {
51  class AbsTrackRep;
52 }
53 
55  : mOptHelix(true)
56  , mOptPreciseFit(true)
57 {
58 }
59 
60 std::vector<PHGenFit2::Track*> PHTpcTrackFollower::followTracks(std::vector<kdfinder::TrackCandidate<double>*>& candidates, PHField* /*B*/, PHTpcLookup* lookup, PHGenFit2::Fitter* fitter)
61 {
62  LOG_DEBUG("tracking.PHTpcTrackFollower.followTracks") << "start";
63 
64  // ----- sort track seeds by n hits desc -----
65  std::sort(candidates.begin(), candidates.end(), [](kdfinder::TrackCandidate<double>* a, kdfinder::TrackCandidate<double>* b) {
66  return a->nhits() > b->nhits();
67  });
68 
69  // ----- process seeds in a loop -----
70  std::vector<PHGenFit2::Track*> gtracks;
71 
72  std::unordered_set<uint64_t> usedClusterKeys;
73  int reused_hits = 0;
74 
75  const int NCANDIDATES = candidates.size();
76  for (int i = 0; i < NCANDIDATES; i++)
77  {
78  LOG_DEBUG("tracking.PHTpcTrackFollower.followTracks") << "+++++ following track " << i << " +++++";
79 
80  // ---- check if candidate hits belong to some already reco-ed track -----
81  const int NHITS = candidates[i]->nhits();
82  reused_hits = 0;
83  for (int j = 0; j < NHITS; j++)
84  {
85  std::vector<double>& hit = candidates[i]->getHit(j);
86  auto it = usedClusterKeys.find(*((int64_t*) &hit[3]));
87  if (it != usedClusterKeys.end())
88  {
89  ++reused_hits;
90  }
91  if (reused_hits > 1)
92  {
93  break;
94  }
95  }
96  // skip seed if it has at least 2 previously used hits
97  if (reused_hits > 1)
98  {
99  LOG_DEBUG("tracking.PHTpcTrackFollower.followTracks") << "more than one hit already used, skipping seed";
100  continue;
101  }
102 
103  // ----- seed looks ok, proceed -----
104  PHGenFit2::Track* gtrack = nullptr;
105  try
106  {
107  gtrack = propagateTrack(candidates[i], lookup, fitter);
108  }
109  catch (...)
110  {
111  LOG_DEBUG("tracking.PHTpcTrackFollower.followTracks") << "exception caught during propagation, skipping track";
112  continue;
113  }
114  LOG_DEBUG("tracking.PHTpcTrackFollower.followTracks") << "track propagated";
115 
116  if (!gtrack || gtrack->getGenFitTrack()->getNumPoints() < 10)
117  {
118  LOG_DEBUG("tracking.PHTpcTrackFollower.followTracks") << "no track or track has less than 10 hit points, skipping";
119  delete gtrack;
120  continue;
121  }
122 
123  genfit::KalmanFitStatus* fstatus = gtrack->getGenFitTrack()->getKalmanFitStatus();
124  double chi2 = fstatus->getChi2(),
125  ndf = fstatus->getNdf();
126  int nfailedpts = fstatus->getNFailedPoints();
127 
128  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "chi2: " << chi2 << ", ndf: " << ndf << ", chi2/ndf: " << chi2 / ndf
129  << ", failed: " << nfailedpts << ", pts: " << gtrack->getGenFitTrack()->getNumPoints();
130 
131  if (nfailedpts > 0)
132  {
133  LOG_DEBUG("tracking.PHTpcTrackFollower.followTracks") << "track has " << nfailedpts << " failed points, skipping";
134  delete gtrack;
135  continue;
136  }
137 
138  gtracks.push_back(gtrack);
139  LOG_DEBUG("tracking.PHTpcTrackFollower.followTracks") << "track added";
140  // ----- mark hits as used in track reco -----
141  const std::vector<TrkrDefs::cluskey>& cluskeys = gtrack->get_cluster_keys();
142  for (TrkrDefs::cluskey cluskey : cluskeys)
143  {
144  usedClusterKeys.insert(cluskey);
145  }
146  }
147 
148  return gtracks;
149 }
150 
153 {
154  //----- candidate hits are pre-sorted, inner R -> outer R -----
155  std::vector<std::vector<double>>& hits = candidate->getHits();
156 
157  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "ctrack nhits: " << hits.size();
158 
159  // TEST: check pos, mom of candidate vs pos, mom of gtrack
160  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "ctrack mom: " << candidate->momentum();
161 
162  std::vector<double> cmom = candidate->getMomForHit(0);
163  std::vector<double> cpos = candidate->getPosForHit(0);
164  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack")
165  << "ctrack pos: " << cpos[0] << ", " << cpos[1] << ", " << cpos[2] << " | "
166  << "mom: " << cmom[0] << ", " << cmom[1] << ", " << cmom[2];
167 
168  // ----- make genfit track out of a track candidate -----
169  PHGenFit2::Track* gtrack = candidate_to_genfit(candidate);
170  try
171  {
172  gtrack->getGenFitTrack()->checkConsistency();
173  }
174  catch (...)
175  {
176  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "Found inconsistent track, skipping";
177  return nullptr;
178  }
179  // ----- make initial fit of the track -----
180  try
181  {
182  if (mOptPreciseFit)
183  {
184  fitter->processTrack5(gtrack);
185  }
186  else
187  {
188  fitter->processTrack1(gtrack);
189  }
190  }
191  catch (...)
192  {
193  // genfit exception thrown
194  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "GenFit exception on track fit, skipping";
195  return nullptr;
196  }
197 
198  if (!gtrack->getGenFitTrack()->hasFitStatus())
199  {
200  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "Track has not been fit, skipping";
201  return nullptr;
202  }
203 
204  genfit::KalmanFitStatus* fstatus = gtrack->getGenFitTrack()->getKalmanFitStatus();
205  double chi2 = fstatus->getChi2(), ndf = fstatus->getNdf();
206  int nfailedpts = fstatus->getNFailedPoints();
207 
208  if (nfailedpts > 0)
209  {
210  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "Track has failed points, skipping";
211  return nullptr;
212  }
213 
214  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "SEED chi2: " << chi2 << ", ndf: " << ndf << ", chi2/ndf: " << chi2 / ndf
215  << ", failed: " << nfailedpts << ", pts: " << gtrack->getGenFitTrack()->getNumPoints();
216 
217  //----- check if track has max possible hits already -----
218  if (hits.size() >= PHTpcConst::TPC_LAYERS_MAX)
219  {
220  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "track has max possible points from seed, no propagation needed";
221  return gtrack;
222  }
223 
224  //----- propagate track outwards using KDTree lookups -----
225  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "follow track outwards";
226  int rc1 = 0;
227  try
228  {
229  rc1 = followTrack(gtrack, lookup, fitter, 1);
230  }
231  catch (...)
232  {
233  // bad track?
234  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "Errors when following track outwards, skipping";
235  return nullptr;
236  }
237 
238  // ---- refit track again ----
239  if (rc1 >= 2)
240  {
241  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "new hits added during outwards following, refitting track";
242  try
243  {
244  // FIXME: 1-pass imprecise fit or partial fit to keep inwards prop precise?
245  fitter->processTrack1(gtrack);
246  }
247  catch (...)
248  {
249  }
250  }
251 
252  //----- propagate track inwards using KDTree lookups ------
253  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "follow track inwards";
254  int rc2 = 0;
255  try
256  {
257  rc2 = followTrack(gtrack, lookup, fitter, -1);
258  }
259  catch (...)
260  {
261  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "Errors when following track inwards";
262  }
263 
264  //----- refit completed track -----
265  if (rc1 || rc2)
266  {
267  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "new hits added during inwards following, refitting track";
268  try
269  {
270  if (mOptPreciseFit)
271  {
272  fitter->processTrack5(gtrack);
273  }
274  else
275  {
276  fitter->processTrack1(gtrack);
277  }
278  }
279  catch (...)
280  {
281  }
282  }
283 
284  try
285  {
286  gtrack->getGenFitTrack()->checkConsistency();
287  }
288  catch (...)
289  {
290  LOG_DEBUG("tracking.PHTpcTrackFollower.propagateTrack") << "got inconsistent track, skipping";
291  return nullptr;
292  }
293 
294  return gtrack;
295 }
296 
298 {
299  // ----- set RungeKutta track rep with the most common particle hypothesis pi- -----
300  genfit::AbsTrackRep* rep = new genfit::RKTrackRep(candidate->sign() > 0 ? 211 : -211 /* Pion PDG ID, as the most common */);
301 
302  // ----- pos, mom from the innermost hit of the track candidate -----
303  std::vector<double> cmom = candidate->getMomForHit(0);
304  std::vector<double> cpos = candidate->getPosForHit(0);
305  TVector3 seed_pos(cpos[0], cpos[1], cpos[2]); // x, y, z
306  TVector3 seed_mom(cmom[0], cmom[1], cmom[2]); // px, py, pz
307 
308  // ----- set up cov matrix -----
309  TMatrixDSym seed_cov(6);
310  seed_cov(0, 0) = 0.1 * 0.1; // dx * dx
311  seed_cov(1, 1) = 0.1 * 0.1; // dy * dy
312  seed_cov(2, 2) = 0.2 * 0.2; // dz * dz
313  seed_cov(3, 3) = cmom[0] * 0.3 * cmom[0] * 0.3; // dpx * dpx
314  seed_cov(4, 4) = cmom[1] * 0.3 * cmom[1] * 0.3; // dpy * dpy
315  seed_cov(5, 5) = cmom[2] * 0.5 * cmom[2] * 0.5; // dpz * dpz
316 
317  PHGenFit2::Track* track = new PHGenFit2::Track(rep, seed_pos, seed_mom, seed_cov);
318 
319  // ----- add track candidate hits as PHGenFit measurements -----
320  const int NHITS = candidate->nhits();
321  std::vector<PHGenFit::Measurement*> measurements;
322  measurements.reserve(NHITS);
323  for (int i = 0; i < NHITS; i++)
324  {
325  std::vector<double>& hit = candidate->getHit(i);
327  measurements.push_back(meas);
328  }
329  track->addMeasurements(measurements);
330 
331  return track;
332 }
333 
335 {
336  TVector3 pos(hit[0], hit[1], hit[2]);
337  TVector3 res(0.1, 0.1, 0.2);
339  meas->set_cluster_key(*((int64_t*) &hit[3]));
340  return meas;
341 }
342 
344 {
345  return get_track_layer(track->getGenFitTrack(), dir);
346 }
347 
348 int PHTpcTrackFollower::get_track_layer(genfit::Track* gftrack, int dir)
349 {
350  int layer = -1; // negative = no layer assigned yet
351  unsigned int GFNUMPOINTS = gftrack->getNumPoints();
352  if (GFNUMPOINTS <= 0)
353  {
354  LOG_ERROR("tracking.PHTpcTrackFollower.get_track_layer") << "track has no points";
355  return layer;
356  }
357  if (!gftrack->hasKalmanFitStatus())
358  {
359  LOG_ERROR("tracking.PHTpcTrackFollower.get_track_layer") << "track has not been fitted";
360  return layer;
361  }
362  try
363  {
364  LOG_DEBUG("tracking.PHTpcTrackFollower.get_track_layer") << "starting layer calc";
365  genfit::TrackPoint* pt = gftrack->getPointWithMeasurementAndFitterInfo(dir == 1 ? GFNUMPOINTS - 1 : 0);
366  if (!pt)
367  {
368  LOG_DEBUG("tracking.PHTpcTrackFollower.get_track_layer") << "track does not have a point with measurement and fitter info";
369  return -1;
370  }
371  LOG_DEBUG("tracking.PHTpcTrackFollower.get_track_layer") << "got point";
372  genfit::KalmanFitterInfo* kfinfo = pt->getKalmanFitterInfo();
373  if (!kfinfo)
374  {
375  LOG_DEBUG("tracking.PHTpcTrackFollower.get_track_layer") << "track point does not have kalman fitter info";
376  return -1;
377  }
378  LOG_DEBUG("tracking.PHTpcTrackFollower.get_track_layer") << "got kfinfo";
379  const genfit::MeasuredStateOnPlane& kfstate = kfinfo->getFittedState();
380  LOG_DEBUG("tracking.PHTpcTrackFollower.get_track_layer") << "got kfstate";
381  TVector3 pos = kfstate.getPos();
382  TVector3 mom = kfstate.getMom();
383  double ptradius = pos.Perp();
384  for (int i = 0; i < PHTpcConst::TPC_LAYERS_MAX; i++)
385  {
386  if (ptradius > (PHTpcConst::TPC_LAYERS_RADIUS[i] - PHTpcConst::TPC_LAYERS_DELTAR[i] / 2.0) &&
388  {
389  layer = i + dir;
390  break;
391  }
392  }
393  }
394  catch (...)
395  {
396  LOG_DEBUG("tracking.PHTpcTrackFollower.get_track_layer") << "error during layer calculation";
397  return -1;
398  }
399  return layer;
400 }
401 
402 std::pair<genfit::MeasuredStateOnPlane*, double> PHTpcTrackFollower::get_projected_coordinate(genfit::Track* gftrack, int dir, double radius)
403 {
404  // project from last fitted point to radius
405  // TODO: for constant field is faster to use helix crossing cylinder - analytical
406  // kdfinder::Helix helix( const TVector<T>& p, const TVector<T>& o, T B, T q )
407  double pathlength = 0;
408  genfit::MeasuredStateOnPlane* state = 0;
409  try
410  {
411  genfit::TrackPoint* pt = gftrack->getPointWithMeasurementAndFitterInfo(dir == 1 ? -1 : 0);
412  if (!pt)
413  {
414  LOG_DEBUG("tracking.PHTpcTrackFollower.get_projected_coordinate") << "no trackpoint with measurement and fitter info";
415  return std::make_pair(nullptr, 0);
416  }
417  genfit::KalmanFitterInfo* kfinfo = pt->getKalmanFitterInfo();
418  if (!kfinfo)
419  {
420  LOG_DEBUG("tracking.PHTpcTrackFollower.get_projected_coordinate") << "trackpoint does not have Kalman Fitter info";
421  return std::make_pair(nullptr, 0);
422  }
423  const genfit::MeasuredStateOnPlane& kfstate = kfinfo->getFittedState();
424 
425  if (mOptHelix)
426  {
427  // make const-field helix projection first to save time if loopers present
428  // FIXME: check if tracks are still correct ( i.e. only loopers are split, but not primaries! )
429  TVector3 pos, mom;
430  kfstate.getPosMom(pos, mom);
431  double charge = kfstate.getCharge();
432  double posX = pos.X(), posY = pos.Y(), posZ = pos.Z(), Bx, By, Bz;
433  genfit::FieldManager::getInstance()->getFieldVal(posX, posY, posZ, Bx, By, Bz);
434  kdfinder::Helix<double> h(kdfinder::TVector<double>(mom.X(), mom.Y(), mom.Z()),
435  kdfinder::TVector<double>(posX, posY, posZ), Bz, charge < 0 ? -1 : 1);
436  std::pair<double, double> s = h.pathLength(radius);
437  if ((!std::isnan(s.first) && s.first < 999.9) || (!std::isnan(s.second) && s.second < 999.9))
438  {
439  LOG_DEBUG("tracking.PHTpcTrackFollower.get_projected_coordinate") << "helix projects to cylinder: " << s.first << ", " << s.second;
440  state = kfstate.clone();
441  pathlength = state->extrapolateToCylinder(radius,
442  TVector3(0., 0., 0.),
443  TVector3(0., 0., 1.),
444  false,
445  false);
446  }
447  else
448  {
449  LOG_DEBUG("tracking.PHTpcTrackFollower.get_projected_coordinate") << "helix does not project to cylinder";
450  // does not project to cylinder :(
451  if (state)
452  {
453  delete state;
454  state = 0;
455  }
456  return std::make_pair(nullptr, 0);
457  }
458  }
459  else
460  { // mOptHelix = false
461 
462  // ALT: no-helix projection, painfully slow for loopers that change direction
463  // GenFit does up to x1000 iterations to see if track projects on helix
464  state = kfstate.clone();
465  pathlength = state->extrapolateToCylinder(radius,
466  TVector3(0., 0., 0.),
467  TVector3(0., 0., 1.),
468  false,
469  false);
470  }
471  }
472  catch (...)
473  {
474  // can't project to helix
475  if (state)
476  {
477  delete state;
478  state = 0;
479  }
480  return std::make_pair(nullptr, 0);
481  }
482  LOG_DEBUG("tracking.PHTpcTrackFollower.get_projected_coordinate") << "pathlength is: " << pathlength;
483  return std::make_pair(state, pathlength);
484 }
485 
486 std::pair<genfit::MeasuredStateOnPlane*, double> PHTpcTrackFollower::get_projected_coordinate(genfit::Track* gftrack, int dir, const TVector3& point)
487 {
488  // project from last fitted point to radius
489  double pathlength = 0;
490  genfit::MeasuredStateOnPlane* state = nullptr;
491 
492  try
493  {
494  genfit::TrackPoint* pt = gftrack->getPointWithMeasurementAndFitterInfo(dir == 1 ? -1 : 0);
495  if (!pt)
496  {
497  LOG_DEBUG("tracking.PHTpcTrackFollower.get_projected_coordinate") << "no trackpoint with measurement and fitter info";
498  return std::make_pair(nullptr, 0);
499  }
500  genfit::KalmanFitterInfo* kfinfo = pt->getKalmanFitterInfo();
501  if (!kfinfo)
502  {
503  LOG_DEBUG("tracking.PHTpcTrackFollower.get_projected_coordinate") << "trackpoint does not have Kalman Fitter info";
504  return std::make_pair(nullptr, 0);
505  }
506  const genfit::MeasuredStateOnPlane& kfstate = kfinfo->getFittedState();
507 
508  state = kfstate.clone();
509  pathlength = state->extrapolateToPoint(
510  point,
511  false,
512  false);
513  }
514  catch (...)
515  {
516  // material effects exception suppression
517  }
518 
519  LOG_DEBUG("tracking.PHTpcTrackFollower.get_projected_coordinate") << "pathlength is: " << pathlength;
520  return std::make_pair(state, pathlength);
521 }
522 
524 {
525  if (!track)
526  {
527  LOG_ERROR("tracking.PHTpcTrackFollower.follow_track") << "cannot follow track, empty pointer!";
528  return 0;
529  }
530 
531  genfit::Track* gftrack = track->getGenFitTrack();
532 
533  int nHitsAdded = 0;
534 
535  LOG_DEBUG("tracking.PHTpcTrackFollower.follow_track") << "about to get layer info";
536  int layer = -1;
537  try
538  {
539  layer = get_track_layer(gftrack, dir);
540  }
541  catch (...)
542  {
543  LOG_DEBUG("tracking.PHTpcTrackFollower.follow_track") << "cannot get track layer, skipping propagation";
544  return 0;
545  }
546  LOG_DEBUG("tracking.PHTpcTrackFollower.follow_track") << "got layer info: " << layer;
547 
548  if (!(layer >= 0 && layer < PHTpcConst::TPC_LAYERS_MAX))
549  {
550  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "track goes beyond min/max layer, stop propagation. Layer: " << layer;
551  return 0;
552  }
553  LOG_DEBUG("tracking.PHTpcTrackFollower.follow_track") << "layer info is good";
554 
555  while (1)
556  {
557  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "dir: " << dir << ", propagating track to layer: " << layer << ", " << PHTpcConst::TPC_LAYERS_RADIUS[layer] << " cm";
558  std::pair<genfit::MeasuredStateOnPlane*, double> p = get_projected_coordinate(gftrack, dir, PHTpcConst::TPC_LAYERS_RADIUS[layer]);
559  if (!p.first)
560  {
561  break;
562  } // can't project to cylinder, likely curler track
563  TVector3 pos = p.first->getPos();
564  delete p.first;
565 
566  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "projected position: " << pos.X() << ", " << pos.Y() << ", " << pos.Z() << ", radius: " << pos.Perp();
567 
568  size_t nMatches = 0;
569  std::vector<std::vector<double>*> matches = lookup->find(pos.X(), pos.Y(), pos.Z(), PHTpcConst::TPC_LAYERS_DELTAR[layer] / 2.0, nMatches);
570  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "lookup returned " << nMatches << " hits nearby, size: " << matches.size();
571 
572  // options: no hits found, one hit found, several hits found
573  if (nMatches > 0)
574  {
575  if (nMatches == 1)
576  {
577  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "processing one hit";
578  std::vector<double>* hit = matches[0];
579  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "hit x: " << (*hit)[0] << ", y: " << (*hit)[1] << ", z: " << (*hit)[2];
580  // add hit as a measurement to the track, check chi2, remove hit if chi2 = bad
581  std::pair<genfit::MeasuredStateOnPlane*, double> p2 = get_projected_coordinate(gftrack, dir, TVector3((*hit)[0], (*hit)[1], (*hit)[2]));
582  if (!p2.first)
583  {
584  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "cannot project to point, skipping";
585  break;
586  }
587  TVector3 pos2 = p2.first->getPos();
588  delete p2.first;
589  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "projected point: " << pos2.X() << ", " << pos2.Y() << ", " << pos2.Z() << ", radius: " << pos.Perp();
590  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "distance to hit: " << std::sqrt(std::pow(pos2.X() - (*hit)[0], 2) + std::pow(pos2.Y() - (*hit)[1], 2) + std::pow(pos2.Z() - (*hit)[2], 2));
591 
592  //FIXME: convert hit to measurement and add measurement to track
594  LOG_DEBUG_IF("tracking.PHTpcTrackFollower.followTrack", !meas) << "bad, no measurement";
595  if (dir == 1)
596  {
597  track->addMeasurement(meas, -1);
598  try
599  {
600  fitter->processTrackPartially(gftrack, -2, -1);
601  }
602  catch (...)
603  {
604  }
605  }
606  else
607  {
608  track->addMeasurement(meas, 0);
609  try
610  {
611  fitter->processTrackPartially(gftrack, 1, 0);
612  }
613  catch (...)
614  {
615  }
616  }
617  ++nHitsAdded;
618  }
619  else
620  {
621  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "processing multiple hits";
622  // project track to point, sort matches by pathlength ASC, try hits one by one, choose best hit
623  int ind = -1;
624  double dist = std::numeric_limits<double>::max();
625  for (unsigned int i = 0; i < nMatches; i++)
626  {
627  std::vector<double>* hit = matches[i];
628  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "hit x: " << (*hit)[0] << ", y: " << (*hit)[1] << ", z: " << (*hit)[2];
629  std::pair<genfit::MeasuredStateOnPlane*, double> p2 = get_projected_coordinate(gftrack, dir, TVector3((*hit)[0], (*hit)[1], (*hit)[2]));
630  if (!p2.first)
631  {
632  continue;
633  }
634  TVector3 pos2 = p2.first->getPos();
635  delete p2.first;
636  double dist2 = std::sqrt(std::pow(pos2.X() - (*hit)[0], 2) + std::pow(pos2.Y() - (*hit)[1], 2) + std::pow(pos2.Z() - (*hit)[2], 2));
637  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "projected point: " << pos2.X() << ", " << pos2.Y() << ", " << pos2.Z() << ", radius: " << pos.Perp();
638  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "distance to hit: " << dist2;
639  if (dist2 < dist)
640  {
641  ind = i;
642  dist = dist2;
643  }
644  }
645  if (ind != -1)
646  {
648  if (dir == 1)
649  {
650  track->addMeasurement(meas, -1);
651  try
652  {
653  fitter->processTrackPartially(gftrack, -2, -1);
654  }
655  catch (...)
656  {
657  }
658  }
659  else
660  {
661  track->addMeasurement(meas, 0);
662  try
663  {
664  fitter->processTrackPartially(gftrack, 1, 0);
665  }
666  catch (...)
667  {
668  }
669  }
670  ++nHitsAdded; // FIXME: only if hit is really added
671  }
672  }
673  } // else => no hits, scan next layer
674 
675  layer += dir;
676  if (layer < 0 || layer >= PHTpcConst::TPC_LAYERS_MAX)
677  {
678  LOG_DEBUG("tracking.PHTpcTrackFollower.followTrack") << "done with " << dir << " following, layer: " << layer;
679  break;
680  }
681  }
682 
683  return nHitsAdded; // number of hits added to the track
684 }