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>
27 #include <GenFit/TrackPoint.h>
29 #include <TMatrixDSymfwd.h>
30 #include <TMatrixTSym.h>
33 #include <log4cpp/CategoryStream.hh>
40 #include <unordered_set>
56 , mOptPreciseFit(
true)
62 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTracks") <<
"start";
66 return a->
nhits() >
b->nhits();
70 std::vector<PHGenFit2::Track*> gtracks;
72 std::unordered_set<uint64_t> usedClusterKeys;
75 const int NCANDIDATES = candidates.size();
76 for (
int i = 0; i < NCANDIDATES; i++)
78 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTracks") <<
"+++++ following track " << i <<
" +++++";
81 const int NHITS = candidates[i]->nhits();
83 for (
int j = 0; j < NHITS; j++)
85 std::vector<double>& hit = candidates[i]->getHit(j);
86 auto it = usedClusterKeys.find(*((int64_t*) &hit[3]));
87 if (
it != usedClusterKeys.end())
99 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTracks") <<
"more than one hit already used, skipping seed";
111 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTracks") <<
"exception caught during propagation, skipping track";
114 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTracks") <<
"track propagated";
118 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTracks") <<
"no track or track has less than 10 hit points, skipping";
123 genfit::KalmanFitStatus* fstatus = gtrack->
getGenFitTrack()->getKalmanFitStatus();
124 double chi2 = fstatus->getChi2(),
125 ndf = fstatus->getNdf();
126 int nfailedpts = fstatus->getNFailedPoints();
128 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"chi2: " << chi2 <<
", ndf: " << ndf <<
", chi2/ndf: " << chi2 / ndf
129 <<
", failed: " << nfailedpts <<
", pts: " << gtrack->
getGenFitTrack()->getNumPoints();
133 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTracks") <<
"track has " << nfailedpts <<
" failed points, skipping";
138 gtracks.push_back(gtrack);
139 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTracks") <<
"track added";
141 const std::vector<TrkrDefs::cluskey>& cluskeys = gtrack->
get_cluster_keys();
144 usedClusterKeys.insert(
cluskey);
155 std::vector<std::vector<double>>&
hits = candidate->
getHits();
157 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"ctrack nhits: " << hits.size();
160 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"ctrack mom: " << candidate->
momentum();
164 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack")
165 <<
"ctrack pos: " << cpos[0] <<
", " << cpos[1] <<
", " << cpos[2] <<
" | "
166 <<
"mom: " << cmom[0] <<
", " << cmom[1] <<
", " << cmom[2];
176 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"Found inconsistent track, skipping";
194 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"GenFit exception on track fit, skipping";
200 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"Track has not been fit, skipping";
204 genfit::KalmanFitStatus* fstatus = gtrack->
getGenFitTrack()->getKalmanFitStatus();
205 double chi2 = fstatus->getChi2(), ndf = fstatus->getNdf();
206 int nfailedpts = fstatus->getNFailedPoints();
210 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"Track has failed points, skipping";
214 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"SEED chi2: " << chi2 <<
", ndf: " << ndf <<
", chi2/ndf: " << chi2 / ndf
215 <<
", failed: " << nfailedpts <<
", pts: " << gtrack->
getGenFitTrack()->getNumPoints();
220 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"track has max possible points from seed, no propagation needed";
225 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"follow track outwards";
234 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"Errors when following track outwards, skipping";
241 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"new hits added during outwards following, refitting track";
253 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"follow track inwards";
261 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"Errors when following track inwards";
267 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"new hits added during inwards following, refitting track";
290 LOG_DEBUG(
"tracking.PHTpcTrackFollower.propagateTrack") <<
"got inconsistent track, skipping";
300 genfit::AbsTrackRep* rep =
new genfit::RKTrackRep(candidate->
sign() > 0 ? 211 : -211 );
305 TVector3 seed_pos(cpos[0], cpos[1], cpos[2]);
306 TVector3 seed_mom(cmom[0], cmom[1], cmom[2]);
309 TMatrixDSym seed_cov(6);
310 seed_cov(0, 0) = 0.1 * 0.1;
311 seed_cov(1, 1) = 0.1 * 0.1;
312 seed_cov(2, 2) = 0.2 * 0.2;
313 seed_cov(3, 3) = cmom[0] * 0.3 * cmom[0] * 0.3;
314 seed_cov(4, 4) = cmom[1] * 0.3 * cmom[1] * 0.3;
315 seed_cov(5, 5) = cmom[2] * 0.5 * cmom[2] * 0.5;
320 const int NHITS = candidate->
nhits();
321 std::vector<PHGenFit::Measurement*> measurements;
322 measurements.reserve(NHITS);
323 for (
int i = 0; i < NHITS; i++)
325 std::vector<double>& hit = candidate->
getHit(i);
327 measurements.push_back(meas);
336 TVector3
pos(hit[0], hit[1], hit[2]);
337 TVector3 res(0.1, 0.1, 0.2);
351 unsigned int GFNUMPOINTS = gftrack->getNumPoints();
352 if (GFNUMPOINTS <= 0)
354 LOG_ERROR(
"tracking.PHTpcTrackFollower.get_track_layer") <<
"track has no points";
357 if (!gftrack->hasKalmanFitStatus())
359 LOG_ERROR(
"tracking.PHTpcTrackFollower.get_track_layer") <<
"track has not been fitted";
364 LOG_DEBUG(
"tracking.PHTpcTrackFollower.get_track_layer") <<
"starting layer calc";
365 genfit::TrackPoint*
pt = gftrack->getPointWithMeasurementAndFitterInfo(dir == 1 ? GFNUMPOINTS - 1 : 0);
368 LOG_DEBUG(
"tracking.PHTpcTrackFollower.get_track_layer") <<
"track does not have a point with measurement and fitter info";
371 LOG_DEBUG(
"tracking.PHTpcTrackFollower.get_track_layer") <<
"got point";
372 genfit::KalmanFitterInfo* kfinfo = pt->getKalmanFitterInfo();
375 LOG_DEBUG(
"tracking.PHTpcTrackFollower.get_track_layer") <<
"track point does not have kalman fitter info";
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();
396 LOG_DEBUG(
"tracking.PHTpcTrackFollower.get_track_layer") <<
"error during layer calculation";
407 double pathlength = 0;
408 genfit::MeasuredStateOnPlane* state = 0;
411 genfit::TrackPoint*
pt = gftrack->getPointWithMeasurementAndFitterInfo(dir == 1 ? -1 : 0);
414 LOG_DEBUG(
"tracking.PHTpcTrackFollower.get_projected_coordinate") <<
"no trackpoint with measurement and fitter info";
415 return std::make_pair(
nullptr, 0);
417 genfit::KalmanFitterInfo* kfinfo = pt->getKalmanFitterInfo();
420 LOG_DEBUG(
"tracking.PHTpcTrackFollower.get_projected_coordinate") <<
"trackpoint does not have Kalman Fitter info";
421 return std::make_pair(
nullptr, 0);
423 const genfit::MeasuredStateOnPlane& kfstate = kfinfo->getFittedState();
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);
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))
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.),
449 LOG_DEBUG(
"tracking.PHTpcTrackFollower.get_projected_coordinate") <<
"helix does not project to cylinder";
456 return std::make_pair(
nullptr, 0);
464 state = kfstate.clone();
465 pathlength = state->extrapolateToCylinder(radius,
466 TVector3(0., 0., 0.),
467 TVector3(0., 0., 1.),
480 return std::make_pair(
nullptr, 0);
482 LOG_DEBUG(
"tracking.PHTpcTrackFollower.get_projected_coordinate") <<
"pathlength is: " << pathlength;
483 return std::make_pair(state, pathlength);
489 double pathlength = 0;
490 genfit::MeasuredStateOnPlane* state =
nullptr;
494 genfit::TrackPoint*
pt = gftrack->getPointWithMeasurementAndFitterInfo(dir == 1 ? -1 : 0);
497 LOG_DEBUG(
"tracking.PHTpcTrackFollower.get_projected_coordinate") <<
"no trackpoint with measurement and fitter info";
498 return std::make_pair(
nullptr, 0);
500 genfit::KalmanFitterInfo* kfinfo = pt->getKalmanFitterInfo();
503 LOG_DEBUG(
"tracking.PHTpcTrackFollower.get_projected_coordinate") <<
"trackpoint does not have Kalman Fitter info";
504 return std::make_pair(
nullptr, 0);
506 const genfit::MeasuredStateOnPlane& kfstate = kfinfo->getFittedState();
508 state = kfstate.clone();
509 pathlength = state->extrapolateToPoint(
519 LOG_DEBUG(
"tracking.PHTpcTrackFollower.get_projected_coordinate") <<
"pathlength is: " << pathlength;
520 return std::make_pair(state, pathlength);
527 LOG_ERROR(
"tracking.PHTpcTrackFollower.follow_track") <<
"cannot follow track, empty pointer!";
535 LOG_DEBUG(
"tracking.PHTpcTrackFollower.follow_track") <<
"about to get layer info";
543 LOG_DEBUG(
"tracking.PHTpcTrackFollower.follow_track") <<
"cannot get track layer, skipping propagation";
546 LOG_DEBUG(
"tracking.PHTpcTrackFollower.follow_track") <<
"got layer info: " <<
layer;
550 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTrack") <<
"track goes beyond min/max layer, stop propagation. Layer: " <<
layer;
553 LOG_DEBUG(
"tracking.PHTpcTrackFollower.follow_track") <<
"layer info is good";
558 std::pair<genfit::MeasuredStateOnPlane*, double>
p =
get_projected_coordinate(gftrack, dir, PHTpcConst::TPC_LAYERS_RADIUS[layer]);
563 TVector3
pos = p.first->getPos();
566 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTrack") <<
"projected position: " << pos.X() <<
", " << pos.Y() <<
", " << pos.Z() <<
", radius: " << pos.Perp();
570 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTrack") <<
"lookup returned " << nMatches <<
" hits nearby, size: " << matches.size();
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];
581 std::pair<genfit::MeasuredStateOnPlane*, double> p2 =
get_projected_coordinate(gftrack, dir, TVector3((*hit)[0], (*hit)[1], (*hit)[2]));
584 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTrack") <<
"cannot project to point, skipping";
587 TVector3 pos2 = p2.first->getPos();
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));
594 LOG_DEBUG_IF(
"tracking.PHTpcTrackFollower.followTrack", !meas) <<
"bad, no measurement";
621 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTrack") <<
"processing multiple hits";
625 for (
unsigned int i = 0; i < nMatches; i++)
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]));
634 TVector3 pos2 = p2.first->getPos();
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;
678 LOG_DEBUG(
"tracking.PHTpcTrackFollower.followTrack") <<
"done with " << dir <<
" following, layer: " <<
layer;