LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
Classes | Public Types | Public Member Functions | Private Member Functions | Private Attributes | List of all members
trkf::TrackKalmanFitter Class Reference

Fit tracks using Kalman Filter fit+smooth. More...

#include "TrackKalmanFitter.h"

Classes

struct  Config
 

Public Types

using Parameters = fhicl::Table< Config >
 

Public Member Functions

 TrackKalmanFitter (const TrackStatePropagator *prop, bool useRMS, bool sortHitsByPlane, bool sortHitsByWire, bool sortOutputHitsMinLength, bool skipNegProp, bool cleanZigzag, bool rejectHighMultHits, bool rejectHitsNegativeGOF, float hitErr2ScaleFact, bool tryNoSkipWhenFails, bool tryBothDirs, bool pickBestHitOnWire, float maxResidue, float maxResidueFirstHit, float maxChi2, float maxDist, float negDistTolerance, int dumpLevel)
 Constructor from TrackStatePropagator and values of configuration parameters. More...
 
 TrackKalmanFitter (const TrackStatePropagator *prop, Parameters const &p)
 Constructor from TrackStatePropagator and Parameters table. More...
 
bool fitTrack (detinfo::DetectorPropertiesData const &detProp, const recob::TrackTrajectory &traj, int tkID, const SMatrixSym55 &covVtx, const SMatrixSym55 &covEnd, const std::vector< art::Ptr< recob::Hit >> &hits, const double pval, const int pdgid, const bool flipDirection, recob::Track &outTrack, std::vector< art::Ptr< recob::Hit >> &outHits, trkmkr::OptionalOutputs &optionals) const
 Fit track starting from TrackTrajectory. More...
 
bool fitTrack (detinfo::DetectorPropertiesData const &detProp, const Point_t &position, const Vector_t &direction, SMatrixSym55 &trackStateCov, const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< recob::TrajectoryPointFlags > &flags, const int tkID, const double pval, const int pdgid, recob::Track &outTrack, std::vector< art::Ptr< recob::Hit >> &outHits, trkmkr::OptionalOutputs &optionals) const
 Fit track starting from intial position, direction, and flags. More...
 
bool doFitWork (KFTrackState &trackState, detinfo::DetectorPropertiesData const &detProp, std::vector< HitState > &hitstatev, std::vector< recob::TrajectoryPointFlags::Mask_t > &hitflagsv, std::vector< KFTrackState > &fwdPrdTkState, std::vector< KFTrackState > &fwdUpdTkState, std::vector< unsigned int > &hitstateidx, std::vector< unsigned int > &rejectedhsidx, std::vector< unsigned int > &sortedtksidx, bool applySkipClean=true) const
 Function where the core of the fit is performed. More...
 

Private Member Functions

KFTrackState setupInitialTrackState (const Point_t &position, const Vector_t &direction, SMatrixSym55 &trackStateCov, const double pval, const int pdgid) const
 Return track state from intial position, direction, and covariance. More...
 
bool setupInputStates (detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< recob::TrajectoryPointFlags > &flags, std::vector< HitState > &hitstatev, std::vector< recob::TrajectoryPointFlags::Mask_t > &hitflagsv) const
 Setup vectors of HitState and Masks to be used during the fit. More...
 
void sortOutput (std::vector< HitState > &hitstatev, std::vector< KFTrackState > &fwdUpdTkState, std::vector< unsigned int > &hitstateidx, std::vector< unsigned int > &rejectedhsidx, std::vector< unsigned int > &sortedtksidx, std::vector< recob::TrajectoryPointFlags::Mask_t > &hitflagsv, bool applySkipClean=true) const
 Sort the output states. More...
 
bool fillResult (const std::vector< art::Ptr< recob::Hit >> &inHits, const int tkID, const int pdgid, std::vector< HitState > &hitstatev, std::vector< recob::TrajectoryPointFlags::Mask_t > &hitflagsv, std::vector< KFTrackState > &fwdPrdTkState, std::vector< KFTrackState > &fwdUpdTkState, std::vector< unsigned int > &hitstateidx, std::vector< unsigned int > &rejectedhsidx, std::vector< unsigned int > &sortedtksidx, recob::Track &outTrack, std::vector< art::Ptr< recob::Hit >> &outHits, trkmkr::OptionalOutputs &optionals) const
 Fill the output objects. More...
 

Private Attributes

art::ServiceHandle< geo::Geometry const > geom
 
const TrackStatePropagator * propagator
 
bool useRMS_
 
bool sortHitsByPlane_
 
bool sortHitsByWire_
 
bool sortOutputHitsMinLength_
 
bool skipNegProp_
 
bool cleanZigzag_
 
bool rejectHighMultHits_
 
bool rejectHitsNegativeGOF_
 
float hitErr2ScaleFact_
 
bool tryNoSkipWhenFails_
 
bool tryBothDirs_
 
bool pickBestHitOnWire_
 
float maxResidue_
 
float maxResidueFirstHit_
 
float maxChi2_
 
float maxDist_
 
float negDistTolerance_
 
int dumpLevel_
 

Detailed Description

Fit tracks using Kalman Filter fit+smooth.

This algorithm fits tracks using a Kalman Filter forward fit followed by a backward smoothing. The resulting track will feature covariance matrices at start and end positions. Fundamental components of the fit are trkf::KFTrackState and trkf::TrackStatePropagator.

Inputs are: recob::TrackTrajectory, initial covariance, associated hits, momentum estimate, particle id hypothesis. Alternatively, instead of a TrackTrajectory the fit can be input with initial position and direction, and vector of recob::TrajectoryPointFlags.

Outputs are: resulting recob::Track, associated hits, and trkmkr::OptionalOutputs.

For configuration options see TrackKalmanFitter::Config

Author
G. Cerati (FNAL, MicroBooNE)
Date
2017
Version
1.0

Definition at line 56 of file TrackKalmanFitter.h.

Member Typedef Documentation

using trkf::TrackKalmanFitter::Parameters = fhicl::Table<Config>

Definition at line 151 of file TrackKalmanFitter.h.

Constructor & Destructor Documentation

trkf::TrackKalmanFitter::TrackKalmanFitter ( const TrackStatePropagator *  prop,
bool  useRMS,
bool  sortHitsByPlane,
bool  sortHitsByWire,
bool  sortOutputHitsMinLength,
bool  skipNegProp,
bool  cleanZigzag,
bool  rejectHighMultHits,
bool  rejectHitsNegativeGOF,
float  hitErr2ScaleFact,
bool  tryNoSkipWhenFails,
bool  tryBothDirs,
bool  pickBestHitOnWire,
float  maxResidue,
float  maxResidueFirstHit,
float  maxChi2,
float  maxDist,
float  negDistTolerance,
int  dumpLevel 
)
inline

Constructor from TrackStatePropagator and values of configuration parameters.

Definition at line 154 of file TrackKalmanFitter.h.

References cluster::sortHitsByWire().

173  {
174  propagator = prop;
175  useRMS_ = useRMS;
176  sortHitsByPlane_ = sortHitsByPlane;
177  sortHitsByWire_ = sortHitsByWire;
178  sortOutputHitsMinLength_ = sortOutputHitsMinLength;
179  skipNegProp_ = skipNegProp;
180  cleanZigzag_ = cleanZigzag;
181  rejectHighMultHits_ = rejectHighMultHits;
182  rejectHitsNegativeGOF_ = rejectHitsNegativeGOF;
183  hitErr2ScaleFact_ = hitErr2ScaleFact;
184  tryNoSkipWhenFails_ = tryNoSkipWhenFails;
185  tryBothDirs_ = tryBothDirs;
186  pickBestHitOnWire_ = pickBestHitOnWire;
187  maxResidue_ = (maxResidue > 0 ? maxResidue : std::numeric_limits<float>::max());
188  maxResidueFirstHit_ =
189  (maxResidueFirstHit > 0 ? maxResidueFirstHit : std::numeric_limits<float>::max());
190  maxChi2_ = (maxChi2 > 0 ? maxChi2 : std::numeric_limits<float>::max());
191  maxDist_ = (maxDist > 0 ? maxDist : std::numeric_limits<float>::max());
192  negDistTolerance_ = negDistTolerance;
193  dumpLevel_ = dumpLevel;
194  }
cluster::sortHitsByWire
bool sortHitsByWire(art::Ptr< recob::Hit > a, art::Ptr< recob::Hit > b)
Definition: ClusterCheater_module.cc:70
trkf::TrackKalmanFitter::propagator
const TrackStatePropagator * propagator
Definition: TrackKalmanFitter.h:299
trkf::TrackKalmanFitter::TrackKalmanFitter ( const TrackStatePropagator *  prop,
Parameters const &  p 
)
inlineexplicit

Constructor from TrackStatePropagator and Parameters table.

Definition at line 197 of file TrackKalmanFitter.h.

References hits(), and lar::dump::vector().

198  : TrackKalmanFitter(prop,
199  p().useRMS(),
200  p().sortHitsByPlane(),
201  p().sortHitsByWire(),
202  p().sortOutputHitsMinLength(),
203  p().skipNegProp(),
204  p().cleanZigzag(),
205  p().rejectHighMultHits(),
206  p().rejectHitsNegativeGOF(),
207  p().hitErr2ScaleFact(),
208  p().tryNoSkipWhenFails(),
209  p().tryBothDirs(),
210  p().pickBestHitOnWire(),
211  p().maxResidue(),
212  p().maxResidueFirstHit(),
213  p().maxChi2(),
214  p().maxDist(),
215  p().negDistTolerance(),
216  p().dumpLevel())
217  {}
trkf::TrackKalmanFitter::TrackKalmanFitter
TrackKalmanFitter(const TrackStatePropagator *prop, bool useRMS, bool sortHitsByPlane, bool sortHitsByWire, bool sortOutputHitsMinLength, bool skipNegProp, bool cleanZigzag, bool rejectHighMultHits, bool rejectHitsNegativeGOF, float hitErr2ScaleFact, bool tryNoSkipWhenFails, bool tryBothDirs, bool pickBestHitOnWire, float maxResidue, float maxResidueFirstHit, float maxChi2, float maxDist, float negDistTolerance, int dumpLevel)
Constructor from TrackStatePropagator and values of configuration parameters.
Definition: TrackKalmanFitter.h:154
cluster::sortHitsByWire
bool sortHitsByWire(art::Ptr< recob::Hit > a, art::Ptr< recob::Hit > b)
Definition: ClusterCheater_module.cc:70

Member Function Documentation

bool trkf::TrackKalmanFitter::doFitWork ( KFTrackState &  trackState,
detinfo::DetectorPropertiesData const &  detProp,
std::vector< HitState > &  hitstatev,
std::vector< recob::TrajectoryPointFlags::Mask_t > &  hitflagsv,
std::vector< KFTrackState > &  fwdPrdTkState,
std::vector< KFTrackState > &  fwdUpdTkState,
std::vector< unsigned int > &  hitstateidx,
std::vector< unsigned int > &  rejectedhsidx,
std::vector< unsigned int > &  sortedtksidx,
bool  applySkipClean = true 
) const

Function where the core of the fit is performed.

Definition at line 295 of file TrackKalmanFitter.cxx.

References util::abs(), trkf::TrackStatePropagator::BACKWARD, util::begin(), trkf::KFTrackState::chi2(), trkf::KFTrackState::covariance(), recob::TrajectoryPointFlagTraits::DeltaRay, recob::TrajectoryPointFlagTraits::DetectorIssue, larg4::dist(), trkf::TrackStatePropagator::distanceToPlane(), trkf::KFTrackState::dump(), trkf::TrackState::dump(), dumpLevel_, util::end(), recob::TrajectoryPointFlagTraits::ExcludedFromFit, trkf::TrackStatePropagator::FORWARD, geom, recob::TrajectoryPointFlagTraits::HitIgnored, geo::GeometryCore::Iterate(), maxChi2_, maxDist_, geo::GeometryCore::MaxPlanes(), maxResidue_, maxResidueFirstHit_, recob::TrajectoryPointFlagTraits::Merged, trkf::KFTrackState::momentum(), negDistTolerance_, recob::TrajectoryPointFlagTraits::NoPoint, pickBestHitOnWire_, trkf::TrackState::plane(), trkf::KFTrackState::position(), trkf::TrackStatePropagator::propagateToPlane(), propagator, recob::TrajectoryPointFlagTraits::Rejected, trkf::KFTrackState::residual(), trkf::KFTrackState::setCovariance(), recob::TrajectoryPointFlagTraits::Shared, skipNegProp_, sortHitsByPlane_, sortHitsByWire_, sortOutput(), recob::TrajectoryPointFlagTraits::Suspicious, trkf::KFTrackState::trackState(), and trkf::KFTrackState::updateWithHitState().

Referenced by fitTrack().

305 {
306  fwdPrdTkState.clear();
307  fwdUpdTkState.clear();
308  hitstateidx.clear();
309  rejectedhsidx.clear();
310  sortedtksidx.clear();
311  // these three vectors are aligned
312  fwdPrdTkState.reserve(hitstatev.size());
313  fwdUpdTkState.reserve(hitstatev.size());
314  hitstateidx.reserve(hitstatev.size());
315 
316  // keep a copy in case first propagation fails
317  KFTrackState startState = trackState;
318 
319  if (sortHitsByPlane_) {
320  //array of hit indices in planes, keeping the original sorting by plane
321  const unsigned int nplanes = geom->MaxPlanes();
322  std::vector<std::vector<unsigned int>> hitsInPlanes(nplanes);
323  for (unsigned int ihit = 0; ihit < hitstatev.size(); ihit++) {
324  hitsInPlanes[hitstatev[ihit].wireId().Plane].push_back(ihit);
325  }
326  if (sortHitsByWire_) {
327  constexpr geo::TPCID tpcid{0, 0};
328  for (auto const& plane : geom->Iterate<geo::PlaneGeo>(tpcid)) {
329  auto const iplane = plane.ID().Plane;
330  if (plane.GetIncreasingWireDirection().Dot(trackState.momentum()) > 0) {
331  std::sort(hitsInPlanes[iplane].begin(),
332  hitsInPlanes[iplane].end(),
333  [hitstatev](const unsigned int& a, const unsigned int& b) -> bool {
334  return hitstatev[a].wireId().Wire < hitstatev[b].wireId().Wire;
335  });
336  }
337  else {
338  std::sort(hitsInPlanes[iplane].begin(),
339  hitsInPlanes[iplane].end(),
340  [hitstatev](const unsigned int& a, const unsigned int& b) -> bool {
341  return hitstatev[a].wireId().Wire > hitstatev[b].wireId().Wire;
342  });
343  }
344  }
345  }
346 
347  //dump hits sorted in each plane
348  if (dumpLevel_ > 1) {
349  int ch = 0;
350  for (auto const& p : hitsInPlanes) {
351  for (auto h : p) {
352  std::cout << "hit #/Plane/Wire/x/mask: " << ch++ << " " << hitstatev[h].wireId().Plane
353  << " " << hitstatev[h].wireId().Wire << " " << hitstatev[h].hitMeas() << " "
354  << hitflagsv[h] << std::endl;
355  }
356  }
357  }
358 
359  //array of indices, where iterHitsInPlanes[i] is the iterator over hitsInPlanes[i]
360  std::vector<unsigned int> iterHitsInPlanes(nplanes, 0);
361  for (unsigned int p = 0; p < hitstatev.size(); ++p) {
362  if (dumpLevel_ > 1) std::cout << std::endl << "processing hit #" << p << std::endl;
363  if (dumpLevel_ > 1)
364  std::cout << "hit sizes: rej=" << rejectedhsidx.size() << " good=" << hitstateidx.size()
365  << " input=" << hitstatev.size() << std::endl;
366  if (dumpLevel_ > 1) {
367  std::cout << "compute distance from state=" << std::endl;
368  trackState.dump();
369  }
370  int min_plane = -1;
371  double min_dist = DBL_MAX;
372  //find the closest hit according to the sorting in each plane
373  for (unsigned int iplane = 0; iplane < nplanes; ++iplane) {
374  //note: ih is a reference, so when 'continue' we increment iterHitsInPlanes[iplane] and the hit is effectively rejected
375  if (dumpLevel_ > 1)
376  std::cout << "iplane=" << iplane << " nplanes=" << nplanes
377  << " iterHitsInPlanes[iplane]=" << iterHitsInPlanes[iplane]
378  << " hitsInPlanes[iplane].size()=" << hitsInPlanes[iplane].size() << std::endl;
379  for (unsigned int& ih = iterHitsInPlanes[iplane]; ih < hitsInPlanes[iplane].size(); ++ih) {
380  //
381  unsigned int ihit = hitsInPlanes[iplane][ih];
382  if (dumpLevel_ > 1)
383  std::cout << "ih=" << ih << " ihit=" << ihit << " size=" << hitsInPlanes[iplane].size()
384  << std::endl;
385  const auto& hitstate = hitstatev[ihit];
386  const auto& hitflags = hitflagsv[ihit];
387  if (hitflags.isSet(recob::TrajectoryPointFlagTraits::NoPoint) ||
388  hitflags.isSet(recob::TrajectoryPointFlagTraits::ExcludedFromFit) ||
389  hitflags.isSet(recob::TrajectoryPointFlagTraits::Rejected)) {
390  if (dumpLevel_ > 1)
391  std::cout << "rejecting hit flags="
392  << hitflags.isSet(recob::TrajectoryPointFlagTraits::NoPoint) << ", "
393  << hitflags.isSet(recob::TrajectoryPointFlagTraits::ExcludedFromFit) << ", "
394  << hitflags.isSet(recob::TrajectoryPointFlagTraits::Rejected) << " "
395  << hitflags << std::endl;
396  rejectedhsidx.push_back(ihit);
397  continue;
398  }
399  //get distance to measurement surface
400  bool success = true;
401  const double dist =
402  propagator->distanceToPlane(success, trackState.trackState(), hitstate.plane());
403  if (dumpLevel_ > 1)
404  std::cout << "distance to plane " << iplane << " wire=" << hitstate.wireId().Wire
405  << " = " << dist << ", min_dist=" << std::min(min_dist, 999.)
406  << " min_plane=" << min_plane << " success=" << success
407  << " wirepo=" << hitstate.plane().position() << std::endl;
408  if (!success) {
409  rejectedhsidx.push_back(ihit);
410  continue;
411  }
412  if (applySkipClean && skipNegProp_ && fwdUpdTkState.size() > 0 &&
413  dist < negDistTolerance_) {
414  rejectedhsidx.push_back(ihit);
415  continue;
416  }
417  if (dist < min_dist) {
418  min_plane = iplane;
419  min_dist = dist;
420  }
421  break;
422  }
423  }
424  if (dumpLevel_ > 1)
425  std::cout
426  << "pick min_dist=" << min_dist << " min_plane=" << min_plane << " wire="
427  << (min_plane < 0 ?
428  -1 :
429  hitstatev[hitsInPlanes[min_plane][iterHitsInPlanes[min_plane]]].wireId().Wire)
430  << std::endl;
431  //
432  //now we know which is the closest wire: get the hitstate and increment the iterator
433  if (min_plane < 0) {
434  if (rejectedhsidx.size() + hitstateidx.size() == hitstatev.size())
435  break;
436  else
437  continue;
438  }
439  unsigned int ihit = hitsInPlanes[min_plane][iterHitsInPlanes[min_plane]];
440  const auto* hitstate = &hitstatev[ihit];
441  iterHitsInPlanes[min_plane]++;
442  //
443  //propagate to measurement surface
444  bool propok = true;
445  trackState = propagator->propagateToPlane(propok,
446  detProp,
447  trackState.trackState(),
448  hitstate->plane(),
449  true,
450  true,
451  TrackStatePropagator::FORWARD);
452  if (!propok && !(applySkipClean && fwdUpdTkState.size() > 0 && skipNegProp_)) {
453  if (dumpLevel_ > 1) std::cout << "attempt backward prop" << std::endl;
454  trackState = propagator->propagateToPlane(propok,
455  detProp,
456  trackState.trackState(),
457  hitstate->plane(),
458  true,
459  true,
460  TrackStatePropagator::BACKWARD);
461  }
462  if (dumpLevel_ > 1) {
463  std::cout << "hit state " << std::endl;
464  hitstate->dump();
465  std::cout << "propagation result=" << propok << std::endl;
466  std::cout << "propagated state " << std::endl;
467  trackState.dump();
468  std::cout << "propagated planarity="
469  << hitstate->plane().direction().Dot(hitstate->plane().position() -
470  trackState.position())
471  << std::endl;
472  std::cout << "residual=" << trackState.residual(*hitstate)
473  << " chi2=" << trackState.chi2(*hitstate) << std::endl;
474  }
475  if (propok) {
476 
477  //reject the first hit if its residue is too large (due to e.g. spurious hit or bad hit sorting)
478  if (applySkipClean && fwdUpdTkState.size() == 0 &&
479  std::abs(trackState.residual(*hitstate)) > maxResidueFirstHit_) {
480  if (dumpLevel_ > 1)
481  std::cout << "rejecting first hit with residual=" << trackState.residual(*hitstate)
482  << std::endl;
483  rejectedhsidx.push_back(ihit);
484  trackState = startState;
485  continue;
486  }
487 
488  //if there is >1 consecutive hit on the same wire, choose the one with best chi2 and exclude the others (should this be optional?)
489  if (pickBestHitOnWire_) {
490  auto wire = hitstate->wireId().Wire;
491  float min_chi2_wire = trackState.chi2(*hitstate);
492  for (unsigned int jh = iterHitsInPlanes[min_plane]; jh < hitsInPlanes[min_plane].size();
493  ++jh) {
494  const unsigned int jhit = hitsInPlanes[min_plane][jh];
495  const auto& jhitstate = hitstatev[jhit];
496  if (jhitstate.wireId().Wire != wire) break;
497  if (ihit != jhit) iterHitsInPlanes[min_plane]++;
498  float chi2 = trackState.chi2(jhitstate);
499  if (dumpLevel_ > 1 && ihit != jhit)
500  std::cout << "additional hit on the same wire with jhit=" << jhit << " chi2=" << chi2
501  << " ihit=" << ihit << " min_chi2_wire=" << min_chi2_wire << std::endl;
502  if (chi2 < min_chi2_wire) {
503  min_chi2_wire = chi2;
504  if (dumpLevel_ > 1 && ihit != jhit) {
505  std::cout << "\tnow using ";
506  jhitstate.dump();
507  }
508  if (ihit != jhit) rejectedhsidx.push_back(ihit);
509  ihit = jhit;
510  }
511  else {
512  rejectedhsidx.push_back(jhit);
513  }
514  }
515  }
516 
517  hitstate = &hitstatev[ihit];
518  auto& hitflags = hitflagsv[ihit];
519 
520  //reject hits failing maxResidue, maxChi2, or maxDist cuts
521  if (fwdUpdTkState.size() > 0 && applySkipClean &&
522  (std::abs(trackState.residual(*hitstate)) > maxResidue_ ||
523  trackState.chi2(*hitstate) > maxChi2_ || min_dist > maxDist_)) {
524  //
525  if (dumpLevel_ > 1)
526  std::cout << "rejecting hit with res=" << std::abs(trackState.residual(*hitstate))
527  << " chi2=" << trackState.chi2(*hitstate) << " dist=" << min_dist
528  << std::endl;
529  // reset the current state, do not update the hit, and mark as excluded from the fit
530  if (fwdUpdTkState.size() > 0)
531  trackState = fwdUpdTkState.back();
532  else
533  trackState = startState;
534  rejectedhsidx.push_back(ihit);
535  hitflags.set(recob::TrajectoryPointFlagTraits::ExcludedFromFit);
536  hitflags.set(
537  recob::TrajectoryPointFlagTraits::
538  NoPoint); //fixme: this is only for event display, also needed by current definition of ValidPoint
539  continue;
540  }
541 
542  hitstateidx.push_back(ihit);
543  fwdPrdTkState.push_back(trackState);
544 
545  //hits with problematic flags are kept but not used for update and flagged as excluded from fit
546  if (hitflags.isSet(recob::TrajectoryPointFlagTraits::HitIgnored) ||
547  hitflags.isSet(recob::TrajectoryPointFlagTraits::Merged) ||
548  hitflags.isSet(recob::TrajectoryPointFlagTraits::Shared) ||
549  hitflags.isSet(recob::TrajectoryPointFlagTraits::DeltaRay) ||
550  hitflags.isSet(recob::TrajectoryPointFlagTraits::DetectorIssue) ||
551  hitflags.isSet(recob::TrajectoryPointFlagTraits::Suspicious)) {
552  //
553  // reset the current state, do not update the hit, and mark as excluded from the fit
554  if (fwdUpdTkState.size() > 0)
555  trackState = fwdUpdTkState.back();
556  else
557  trackState = startState;
558  fwdUpdTkState.push_back(trackState);
559  hitflags.set(recob::TrajectoryPointFlagTraits::ExcludedFromFit);
560  hitflags.set(
561  recob::TrajectoryPointFlagTraits::
562  NoPoint); //fixme: this is only for event display, also needed by current definition of ValidPoint
563  continue;
564  }
565  //now update the forward fitted track
566  bool upok = trackState.updateWithHitState(*hitstate);
567  if (upok == 0) {
568  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__;
569  return false;
570  }
571  if (dumpLevel_ > 1) {
572  std::cout << "updated state " << std::endl;
573  trackState.dump();
574  }
575  fwdUpdTkState.push_back(trackState);
576  }
577  else {
578  if (dumpLevel_ > 0)
579  std::cout << "WARNING: forward propagation failed. Skip this hit..." << std::endl;
580  // restore the last successful propagation
581  if (fwdPrdTkState.size() > 0)
582  trackState = fwdPrdTkState.back();
583  else
584  trackState = startState;
585  rejectedhsidx.push_back(ihit);
586  continue;
587  }
588  if (rejectedhsidx.size() + hitstateidx.size() == hitstatev.size()) break;
589  }
590  }
591  else {
592  for (unsigned int ihit = 0; ihit < hitstatev.size(); ++ihit) {
593  const auto& hitstate = hitstatev[ihit];
594  if (dumpLevel_ > 1) {
595  std::cout << std::endl << "processing hit #" << ihit << std::endl;
596  hitstate.dump();
597  }
598  auto& hitflags = hitflagsv[ihit];
599  if (hitflags.isSet(recob::TrajectoryPointFlagTraits::ExcludedFromFit) ||
600  hitflags.isSet(recob::TrajectoryPointFlagTraits::Rejected)) {
601  rejectedhsidx.push_back(ihit);
602  continue;
603  }
604  bool success = true;
605  const double dist =
606  propagator->distanceToPlane(success, trackState.trackState(), hitstate.plane());
607  if (applySkipClean && fwdUpdTkState.size() > 0 && skipNegProp_) {
608  if (dist < 0. || success == false) {
609  if (dumpLevel_ > 0)
610  std::cout << "WARNING: negative propagation distance. Skip this hit..." << std::endl;
611  ;
612  rejectedhsidx.push_back(ihit);
613  continue;
614  }
615  }
616  //propagate to measurement surface
617  bool propok = true;
618  trackState = propagator->propagateToPlane(propok,
619  detProp,
620  trackState.trackState(),
621  hitstate.plane(),
622  true,
623  true,
624  TrackStatePropagator::FORWARD);
625  if (!propok && !(applySkipClean && skipNegProp_))
626  trackState = propagator->propagateToPlane(propok,
627  detProp,
628  trackState.trackState(),
629  hitstate.plane(),
630  true,
631  true,
632  TrackStatePropagator::BACKWARD);
633  if (propok) {
634  hitstateidx.push_back(ihit);
635  fwdPrdTkState.push_back(trackState);
636  //
637  if (hitflags.isSet(recob::TrajectoryPointFlagTraits::HitIgnored) ||
638  hitflags.isSet(recob::TrajectoryPointFlagTraits::Merged) ||
639  hitflags.isSet(recob::TrajectoryPointFlagTraits::Shared) ||
640  hitflags.isSet(recob::TrajectoryPointFlagTraits::DeltaRay) ||
641  hitflags.isSet(recob::TrajectoryPointFlagTraits::DetectorIssue) ||
642  hitflags.isSet(recob::TrajectoryPointFlagTraits::Suspicious) ||
643  (fwdUpdTkState.size() > 0 && applySkipClean &&
644  (std::abs(trackState.residual(hitstate)) > maxResidue_ ||
645  trackState.chi2(hitstate) > maxChi2_ || dist > maxDist_))) {
646  if (dumpLevel_ > 1)
647  std::cout << "rejecting hit with res=" << std::abs(trackState.residual(hitstate))
648  << " chi2=" << trackState.chi2(hitstate) << " dist=" << dist << std::endl;
649  // reset the current state, do not update the hit, mark as excluded from the fit
650  if (fwdUpdTkState.size() > 0)
651  trackState = fwdUpdTkState.back();
652  else
653  trackState = startState;
654  fwdUpdTkState.push_back(trackState);
655  hitflags.set(recob::TrajectoryPointFlagTraits::ExcludedFromFit);
656  hitflags.set(
657  recob::TrajectoryPointFlagTraits::
658  NoPoint); //fixme: this is only for event display, also needed by current definition of ValidPoint
659  continue;
660  }
661  //
662  bool upok = trackState.updateWithHitState(hitstate);
663  if (upok == 0) {
664  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__;
665  return false;
666  }
667  fwdUpdTkState.push_back(trackState);
668  }
669  else {
670  if (dumpLevel_ > 0)
671  std::cout << "WARNING: forward propagation failed. Skip this hit..." << std::endl;
672  ;
673  // restore the last successful propagation
674  if (fwdPrdTkState.size() > 0)
675  trackState = fwdPrdTkState.back();
676  else
677  trackState = startState;
678  rejectedhsidx.push_back(ihit);
679  continue;
680  }
681  } //for (auto hitstate : hitstatev)
682  }
683 
684  if (dumpLevel_ > 2) assert(rejectedhsidx.size() + hitstateidx.size() == hitstatev.size());
685  if (dumpLevel_ > 0) {
686  std::cout << "TRACK AFTER FWD" << std::endl;
687  std::cout << "hit sizes=" << rejectedhsidx.size() << " " << hitstateidx.size() << " "
688  << hitstatev.size() << std::endl;
689  trackState.dump();
690  }
691 
692  //reinitialize trf for backward fit, scale the error to avoid biasing the backward fit
693  trackState.setCovariance(100. * trackState.covariance());
694 
695  startState = trackState;
696 
697  //backward loop over track states and hits in fwdUpdTracks: use hits for backward fit and fwd track states for smoothing
698  int nvalidhits = 0;
699  for (int itk = fwdPrdTkState.size() - 1; itk >= 0; itk--) {
700  auto& fwdPrdTrackState = fwdPrdTkState[itk];
701  auto& fwdUpdTrackState = fwdUpdTkState[itk];
702  const auto& hitstate = hitstatev[hitstateidx[itk]];
703  auto& hitflags = hitflagsv[hitstateidx[itk]];
704  if (dumpLevel_ > 1) {
705  std::cout << std::endl << "processing backward hit #" << itk << std::endl;
706  hitstate.dump();
707  }
708  bool propok = true;
709  trackState = propagator->propagateToPlane(propok,
710  detProp,
711  trackState.trackState(),
712  hitstate.plane(),
713  true,
714  true,
715  TrackStatePropagator::BACKWARD);
716  if (!propok)
717  trackState = propagator->propagateToPlane(propok,
718  detProp,
719  trackState.trackState(),
720  hitstate.plane(),
721  true,
722  true,
723  TrackStatePropagator::FORWARD); //do we want this?
724  //
725  if (dumpLevel_ > 1) {
726  std::cout << "propagation result=" << propok << std::endl;
727  std::cout << "propagated state " << std::endl;
728  trackState.dump();
729  std::cout << "propagated planarity="
730  << hitstate.plane().direction().Dot(hitstate.plane().position() -
731  trackState.position())
732  << std::endl;
733  }
734  if (propok) {
735  //combine forward predicted and backward predicted
736  bool pcombok = fwdPrdTrackState.combineWithTrackState(trackState.trackState());
737  if (pcombok == 0) {
738  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__;
739  return false;
740  }
741  if (dumpLevel_ > 1) {
742  std::cout << "combined prd state " << std::endl;
743  fwdPrdTrackState.dump();
744  }
745  // combine forward updated and backward predicted and update backward predicted, only if the hit was not excluded
746  if (hitflags.isSet(recob::TrajectoryPointFlagTraits::ExcludedFromFit) == 0) {
747  bool ucombok = fwdUpdTrackState.combineWithTrackState(trackState.trackState());
748  if (ucombok == 0) {
749  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__;
750  return false;
751  }
752  if (dumpLevel_ > 1) {
753  std::cout << "combined upd state " << std::endl;
754  fwdUpdTrackState.dump();
755  }
756  bool upok = trackState.updateWithHitState(hitstate);
757  if (upok == 0) {
758  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__;
759  return false;
760  }
761  if (dumpLevel_ > 1) {
762  std::cout << "updated state " << std::endl;
763  trackState.dump();
764  }
765  // Keep a copy in case a future propagation fails
766  startState = trackState;
767  //
768  nvalidhits++;
769  }
770  else {
771  fwdUpdTrackState = fwdPrdTrackState;
772  }
773  }
774  else {
775  // ok, if the backward propagation failed we exclude this point from the rest of the fit,
776  // but we can still use its position from the forward fit, so just mark it as ExcludedFromFit
777  hitflags.set(recob::TrajectoryPointFlagTraits::ExcludedFromFit);
778  hitflags.set(
779  recob::TrajectoryPointFlagTraits::NoPoint); // fixme: this is only for event display, also
780  // needed by current definition of ValidPoint
781  if (dumpLevel_ > 0)
782  std::cout << "WARNING: backward propagation failed. Skip this hit... hitstateidx[itk]="
783  << hitstateidx[itk] << " itk=" << itk << std::endl;
784  ;
785  // restore the last successful propagation
786  trackState = startState;
787  continue;
788  }
789  } //for (unsigned int itk = fwdPrdTkState.size()-1; itk>=0; itk--) {
790 
791  if (nvalidhits < 4) {
792  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__ << " ";
793  return false;
794  }
795 
796  if (dumpLevel_ > 1) std::cout << "sort output with nvalidhits=" << nvalidhits << std::endl;
797 
798  // sort output states
799  sortOutput(
800  hitstatev, fwdUpdTkState, hitstateidx, rejectedhsidx, sortedtksidx, hitflagsv, applySkipClean);
801  size_t nsortvalid = 0;
802  for (auto& idx : sortedtksidx) {
803  auto& hitflags = hitflagsv[hitstateidx[idx]];
804  if (hitflags.isSet(recob::TrajectoryPointFlagTraits::ExcludedFromFit) == 0) nsortvalid++;
805  }
806  if (nsortvalid < 4) {
807  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__ << " ";
808  return false;
809  }
810 
811  if (dumpLevel_ > 2) assert(rejectedhsidx.size() + sortedtksidx.size() == hitstatev.size());
812  return true;
813 }
geo::GeometryCore::Iterate
details::range_type< T > Iterate() const
Initializes the specified ID with the ID of the first cryostat.
Definition: GeometryCore.h:541
recob::TrajectoryPointFlagTraits::Merged
static constexpr Flag_t Merged
The hit might have contribution from particles other than this.
Definition: TrajectoryPointFlags.h:131
recob::TrajectoryPointFlagTraits::Suspicious
static constexpr Flag_t Suspicious
The point reconstruction is somehow questionable.
Definition: TrajectoryPointFlags.h:128
recob::TrajectoryPointFlagTraits::NoPoint
static constexpr Flag_t NoPoint
The trajectory point is not defined.
Definition: TrajectoryPointFlags.h:118
util::abs
constexpr auto abs(T v)
Returns the absolute value of the argument.
Definition: constexpr_math.h:40
recob::TrajectoryPointFlagTraits
Namespace for the trajectory point flags.
Definition: TrajectoryPointFlags.h:57
trkf::TrackKalmanFitter::geom
art::ServiceHandle< geo::Geometry const > geom
Definition: TrackKalmanFitter.h:298
trkf::TrackStatePropagator::propagateToPlane
TrackState propagateToPlane(bool &success, const detinfo::DetectorPropertiesData &detProp, const TrackState &origin, const Plane &target, bool dodedx, bool domcs, PropDirection dir=FORWARD) const
Main function for propagation of a TrackState to a Plane.
Definition: TrackStatePropagator.cxx:29
trkf::TrackStatePropagator::distanceToPlane
double distanceToPlane(bool &success, const Point_t &origpos, const Vector_t &origdir, const Plane &target) const
Distance of a TrackState (Point and Vector) to a Plane, along the TrackState direction.
Definition: TrackStatePropagator.cxx:213
trkf::TrackKalmanFitter::sortOutput
void sortOutput(std::vector< HitState > &hitstatev, std::vector< KFTrackState > &fwdUpdTkState, std::vector< unsigned int > &hitstateidx, std::vector< unsigned int > &rejectedhsidx, std::vector< unsigned int > &sortedtksidx, std::vector< recob::TrajectoryPointFlags::Mask_t > &hitflagsv, bool applySkipClean=true) const
Sort the output states.
Definition: TrackKalmanFitter.cxx:815
util::end
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
geo::GeometryCore::MaxPlanes
unsigned int MaxPlanes() const
Returns the largest number of planes among all TPCs in this detector.
Definition: GeometryCore.cxx:818
geo::PlaneGeo
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:78
trkf::TrackKalmanFitter::propagator
const TrackStatePropagator * propagator
Definition: TrackKalmanFitter.h:299
recob::TrajectoryPointFlagTraits::HitIgnored
static constexpr Flag_t HitIgnored
Hit was not included for the computation of the trajectory.
Definition: TrajectoryPointFlags.h:115
geo::TPCID
The data type to uniquely identify a TPC.
Definition: geo_types.h:381
recob::TrajectoryPointFlagTraits::Rejected
static constexpr Flag_t Rejected
The hit is extraneous to this track.
Definition: TrajectoryPointFlags.h:171
recob::TrajectoryPointFlagTraits::ExcludedFromFit
static constexpr Flag_t ExcludedFromFit
Definition: TrajectoryPointFlags.h:168
larg4::dist
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
Definition: OpFastScintillation.hh:537
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:210
trkf::TrackState::dump
std::ostream & dump(std::ostream &out=std::cout) const
Printout information.
Definition: TrackState.h:144
recob::TrajectoryPointFlagTraits::DetectorIssue
static constexpr Flag_t DetectorIssue
The hit is associated to a problematic channel.
Definition: TrajectoryPointFlags.h:137
util::begin
decltype(auto) constexpr begin(T &&obj)
ADL-aware version of std::begin.
Definition: StdUtils.h:69
recob::TrajectoryPointFlagTraits::DeltaRay
static constexpr Flag_t DeltaRay
The hit might have contribution from a δ ray.
Definition: TrajectoryPointFlags.h:134
recob::TrajectoryPointFlagTraits::Shared
static constexpr Flag_t Shared
The hit is known to be associated also to another trajectory.
Definition: TrajectoryPointFlags.h:140
bool trkf::TrackKalmanFitter::fillResult ( const std::vector< art::Ptr< recob::Hit >> &  inHits,
const int  tkID,
const int  pdgid,
std::vector< HitState > &  hitstatev,
std::vector< recob::TrajectoryPointFlags::Mask_t > &  hitflagsv,
std::vector< KFTrackState > &  fwdPrdTkState,
std::vector< KFTrackState > &  fwdUpdTkState,
std::vector< unsigned int > &  hitstateidx,
std::vector< unsigned int > &  rejectedhsidx,
std::vector< unsigned int > &  sortedtksidx,
recob::Track &  outTrack,
std::vector< art::Ptr< recob::Hit >> &  outHits,
trkmkr::OptionalOutputs &  optionals 
) const
private

Fill the output objects.

Definition at line 940 of file TrackKalmanFitter.cxx.

References trkmkr::TrackCreationBookKeeper::addPoint(), recob::Track::Chi2PerNdof(), trkf::KFTrackState::covariance(), dumpLevel_, e, recob::TrajectoryPointFlagTraits::ExcludedFromFit, trkmkr::TrackCreationBookKeeper::finalizeTrack(), recob::Track::HasValidPoint(), recob::TrajectoryPointFlagTraits::HitIgnored, trkmkr::OptionalOutputs::isTrackFitInfosInit(), util::kBogusD, recob::Track::Length(), recob::Track::MomentumAtPoint(), recob::Track::NextValidPoint(), recob::TrajectoryPointFlagTraits::NoPoint, propagator, recob::TrajectoryPointFlagTraits::Rejected, trkf::TrackStatePropagator::rotateToPlane(), trkmkr::OptionalPointElement::setTrackFitHitInfo(), recob::Track::Start(), recob::Track::StartCovariance(), and recob::Track::StartDirection().

Referenced by fitTrack().

954 {
955  // fill output trajectory objects with smoothed track and its hits
956  int nvalidhits = 0;
957  trkmkr::TrackCreationBookKeeper tcbk(outHits, optionals, tkID, pdgid, true);
958  for (unsigned int p : sortedtksidx) {
959  const auto& trackstate = fwdUpdTkState[p];
960  const auto& hitflags = hitflagsv[hitstateidx[p]];
961  const unsigned int originalPos = hitstateidx[p];
962  if (dumpLevel_ > 2) assert(originalPos < hitstatev.size());
963  //
964  const auto& prdtrack = fwdPrdTkState[p];
965  const auto& hitstate = hitstatev[hitstateidx[p]];
966  if (dumpLevel_ > 2) assert(hitstate.wireId().Plane == inHits[originalPos]->WireID().Plane);
967  //
968  float chi2 =
969  (hitflags.isSet(recob::TrajectoryPointFlagTraits::ExcludedFromFit) ? -1. :
970  prdtrack.chi2(hitstate));
971  if (hitflags.isSet(recob::TrajectoryPointFlagTraits::NoPoint) == false) nvalidhits++;
972  //
973  trkmkr::OptionalPointElement ope;
974  if (optionals.isTrackFitInfosInit()) {
975  ope.setTrackFitHitInfo(recob::TrackFitHitInfo(hitstate.hitMeas(),
976  hitstate.hitMeasErr2(),
977  prdtrack.parameters(),
978  prdtrack.covariance(),
979  hitstate.wireId()));
980  }
981  tcbk.addPoint(trackstate.position(),
982  trackstate.momentum(),
983  inHits[originalPos],
984  recob::TrajectoryPointFlags(originalPos, hitflags),
985  chi2,
986  ope);
987  }
988  if (dumpLevel_ > 0) std::cout << "fillResult nvalidhits=" << nvalidhits << std::endl;
989 
990  // fill also with rejected hits information
991  SMatrixSym55 fakeCov55;
992  for (int i = 0; i < 5; i++)
993  for (int j = i; j < 5; j++)
994  fakeCov55(i, j) = util::kBogusD;
995  for (unsigned int rejidx = 0; rejidx < rejectedhsidx.size(); ++rejidx) {
996  const unsigned int originalPos = rejectedhsidx[rejidx];
997  auto& mask = hitflagsv[rejectedhsidx[rejidx]];
998  mask.set(recob::TrajectoryPointFlagTraits::HitIgnored,
999  recob::TrajectoryPointFlagTraits::NoPoint);
1000  if (mask.isSet(recob::TrajectoryPointFlagTraits::Rejected) == 0)
1001  mask.set(recob::TrajectoryPointFlagTraits::ExcludedFromFit);
1002  //
1003  const auto& hitstate = hitstatev[rejectedhsidx[rejidx]];
1004  if (dumpLevel_ > 2) assert(hitstate.wireId().Plane == inHits[originalPos]->WireID().Plane);
1005  trkmkr::OptionalPointElement ope;
1006  if (optionals.isTrackFitInfosInit()) {
1007  ope.setTrackFitHitInfo(recob::TrackFitHitInfo(
1008  hitstate.hitMeas(),
1009  hitstate.hitMeasErr2(),
1010  SVector5(util::kBogusD, util::kBogusD, util::kBogusD, util::kBogusD, util::kBogusD),
1011  fakeCov55,
1012  hitstate.wireId()));
1013  }
1014  tcbk.addPoint(Point_t(util::kBogusD, util::kBogusD, util::kBogusD),
1015  Vector_t(util::kBogusD, util::kBogusD, util::kBogusD),
1016  inHits[originalPos],
1017  recob::TrajectoryPointFlags(originalPos, mask),
1018  -1.,
1019  ope);
1020  }
1021 
1022  if (dumpLevel_ > 0)
1023  std::cout << "outHits.size()=" << outHits.size() << " inHits.size()=" << inHits.size()
1024  << std::endl;
1025  if (dumpLevel_ > 2) assert(outHits.size() == inHits.size());
1026 
1027  bool propok = true;
1028  KFTrackState resultF =
1029  propagator->rotateToPlane(propok,
1030  fwdUpdTkState[sortedtksidx.front()].trackState(),
1031  Plane(fwdUpdTkState[sortedtksidx.front()].position(),
1032  fwdUpdTkState[sortedtksidx.front()].momentum()));
1033  KFTrackState resultB =
1034  propagator->rotateToPlane(propok,
1035  fwdUpdTkState[sortedtksidx.back()].trackState(),
1036  Plane(fwdUpdTkState[sortedtksidx.back()].position(),
1037  fwdUpdTkState[sortedtksidx.back()].momentum()));
1038 
1039  outTrack =
1040  tcbk.finalizeTrack(SMatrixSym55(resultF.covariance()), SMatrixSym55(resultB.covariance()));
1041 
1042  //if there are points with zero momentum return false
1043  size_t point = 0;
1044  while (outTrack.HasValidPoint(point)) {
1045  if (outTrack.MomentumAtPoint(outTrack.NextValidPoint(point++)) <= 1.0e-9) {
1046  mf::LogWarning("TrackKalmanFitter") << "found point with zero momentum!" << std::endl;
1047  return false;
1048  }
1049  }
1050 
1051  if (dumpLevel_ > 0) {
1052  std::cout << "outTrack vertex=" << outTrack.Start() << "\ndir=" << outTrack.StartDirection()
1053  << "\ncov=\n"
1054  << outTrack.StartCovariance() << "\nlength=" << outTrack.Length()
1055  << "\nchi2/ndof=" << outTrack.Chi2PerNdof() << std::endl;
1056  }
1057 
1058  return true;
1059 }
recob::TrajectoryPointFlagTraits::NoPoint
static constexpr Flag_t NoPoint
The trajectory point is not defined.
Definition: TrajectoryPointFlags.h:118
trkf::Point_t
recob::tracking::Point_t Point_t
Definition: TrackState.h:18
trkf::Vector_t
recob::tracking::Vector_t Vector_t
Definition: TrackState.h:19
recob::Track::StartCovariance
const SMatrixSym55 & StartCovariance() const
Access to covariance matrices.
Definition: Track.h:190
trkf::SMatrixSym55
recob::tracking::SMatrixSym55 SMatrixSym55
Definition: TrackState.h:15
recob::Track::HasValidPoint
bool HasValidPoint(size_t i) const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:145
recob::Track::MomentumAtPoint
double MomentumAtPoint(unsigned int p) const
Definition: Track.h:175
recob::Track::StartDirection
Vector_t StartDirection() const
Access to track direction at different points.
Definition: Track.h:165
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:207
trkf::SVector5
recob::tracking::SVector5 SVector5
Definition: TrackState.h:12
recob::Track::Start
Point_t const & Start() const
Access to track position at different points.
Definition: Track.h:157
recob::TrackFitHitInfo
Object storing per-hit information from a track fit.
Definition: TrackFitHitInfo.h:29
trkmkr::OptionalPointElement
Struct holding point-by-point elements used in OptionalOutputs.
Definition: TrackMaker.h:44
recob::Track::Chi2PerNdof
float Chi2PerNdof() const
Access to various track properties.
Definition: Track.h:209
trkf::TrackStatePropagator::rotateToPlane
TrackState rotateToPlane(bool &success, const TrackState &origin, const Plane &target) const
Rotation of a TrackState to a Plane (zero distance propagation)
Definition: TrackStatePropagator.h:103
trkf::TrackKalmanFitter::propagator
const TrackStatePropagator * propagator
Definition: TrackKalmanFitter.h:299
recob::TrajectoryPointFlagTraits::HitIgnored
static constexpr Flag_t HitIgnored
Hit was not included for the computation of the trajectory.
Definition: TrajectoryPointFlags.h:115
recob::TrajectoryPointFlagTraits::Rejected
static constexpr Flag_t Rejected
The hit is extraneous to this track.
Definition: TrajectoryPointFlags.h:171
recob::TrajectoryPointFlagTraits::ExcludedFromFit
static constexpr Flag_t ExcludedFromFit
Definition: TrajectoryPointFlags.h:168
SMatrixSym55
recob::tracking::SMatrixSym55 SMatrixSym55
Definition: PMAlgTracking.cxx:22
recob::Track::NextValidPoint
size_t NextValidPoint(size_t index) const
Various functions related to the presence and the number of (valid) points.
Definition: Track.h:141
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:210
trkmkr::TrackCreationBookKeeper
Helper class to aid the creation of a recob::Track, keeping data vectors in sync. ...
Definition: TrackCreationBookKeeper.h:32
trkmkr::OptionalOutputs::isTrackFitInfosInit
bool isTrackFitInfosInit()
check initialization of the output vector of TrackFitHitInfos
Definition: TrackMaker.h:147
util::kBogusD
constexpr double kBogusD
obviously bogus double value
Definition: PhysicalConstants.h:71
trkmkr::OptionalPointElement::setTrackFitHitInfo
void setTrackFitHitInfo(recob::TrackFitHitInfo &&aTrackFitHitInfo)
set the recob::TrackFitHitInfo unique_ptr
Definition: TrackMaker.h:47
e
Float_t e
Definition: plot.C:35
trkf::Plane
recob::tracking::Plane Plane
Definition: TrackState.h:17
recob::TrajectoryPointFlags
Set of flags pertaining a point of the track.
Definition: TrajectoryPointFlags.h:304
bool trkf::TrackKalmanFitter::fitTrack ( detinfo::DetectorPropertiesData const &  detProp,
const recob::TrackTrajectory &  traj,
int  tkID,
const SMatrixSym55 &  covVtx,
const SMatrixSym55 &  covEnd,
const std::vector< art::Ptr< recob::Hit >> &  hits,
const double  pval,
const int  pdgid,
const bool  flipDirection,
recob::Track &  outTrack,
std::vector< art::Ptr< recob::Hit >> &  outHits,
trkmkr::OptionalOutputs &  optionals 
) const

Fit track starting from TrackTrajectory.

Definition at line 36 of file TrackKalmanFitter.cxx.

References recob::TrackTrajectory::End(), recob::TrackTrajectory::EndDirection(), recob::TrackTrajectory::Flags(), hits(), tryBothDirs_, recob::TrackTrajectory::Vertex(), and recob::TrackTrajectory::VertexDirection().

Referenced by trkmkr::KalmanFilterFitTrackMaker::makeTrackImpl(), trkf::KalmanFilterTrajectoryFitter::produce(), and trkf::KalmanFilterFinalTrackFitter::produce().

48 {
49  auto position = traj.Vertex();
50  auto direction = traj.VertexDirection();
51 
52  if (tryBothDirs_) {
53  recob::Track fwdTrack;
54  std::vector<art::Ptr<recob::Hit>> fwdHits;
55  trkmkr::OptionalOutputs fwdoptionals;
56  SMatrixSym55 fwdcov = covVtx;
57  bool okfwd = fitTrack(detProp,
58  position,
59  direction,
60  fwdcov,
61  hits,
62  traj.Flags(),
63  tkID,
64  pval,
65  pdgid,
66  fwdTrack,
67  fwdHits,
68  fwdoptionals);
69 
70  recob::Track bwdTrack;
71  std::vector<art::Ptr<recob::Hit>> bwdHits;
72  trkmkr::OptionalOutputs bwdoptionals;
73  SMatrixSym55 bwdcov = covEnd;
74  bool okbwd = fitTrack(detProp,
75  position,
76  -direction,
77  bwdcov,
78  hits,
79  traj.Flags(),
80  tkID,
81  pval,
82  pdgid,
83  bwdTrack,
84  bwdHits,
85  bwdoptionals);
86 
87  if (okfwd == false && okbwd == false) { return false; }
88  else if (okfwd == true && okbwd == true) {
89  if ((fwdTrack.CountValidPoints() / (fwdTrack.Length() * fwdTrack.Chi2PerNdof())) >=
90  (bwdTrack.CountValidPoints() / (bwdTrack.Length() * bwdTrack.Chi2PerNdof()))) {
91  outTrack = fwdTrack;
92  outHits = fwdHits;
93  optionals = std::move(fwdoptionals);
94  }
95  else {
96  outTrack = bwdTrack;
97  outHits = bwdHits;
98  optionals = std::move(bwdoptionals);
99  }
100  }
101  else if (okfwd == true) {
102  outTrack = fwdTrack;
103  outHits = fwdHits;
104  optionals = std::move(fwdoptionals);
105  }
106  else {
107  outTrack = bwdTrack;
108  outHits = bwdHits;
109  optionals = std::move(bwdoptionals);
110  }
111  return true;
112  }
113  else {
114  if (flipDirection) {
115  position = traj.End();
116  direction = -traj.EndDirection();
117  }
118 
119  auto trackStateCov = (flipDirection ? covEnd : covVtx);
120 
121  return fitTrack(detProp,
122  position,
123  direction,
124  trackStateCov,
125  hits,
126  traj.Flags(),
127  tkID,
128  pval,
129  pdgid,
130  outTrack,
131  outHits,
132  optionals);
133  }
134 }
recob::TrackTrajectory::Flags
Flags_t const & Flags() const
Returns all flags.
Definition: TrackTrajectory.h:161
trkf::TrackKalmanFitter::fitTrack
bool fitTrack(detinfo::DetectorPropertiesData const &detProp, const recob::TrackTrajectory &traj, int tkID, const SMatrixSym55 &covVtx, const SMatrixSym55 &covEnd, const std::vector< art::Ptr< recob::Hit >> &hits, const double pval, const int pdgid, const bool flipDirection, recob::Track &outTrack, std::vector< art::Ptr< recob::Hit >> &outHits, trkmkr::OptionalOutputs &optionals) const
Fit track starting from TrackTrajectory.
Definition: TrackKalmanFitter.cxx:36
recob::TrackTrajectory::VertexDirection
Vector_t VertexDirection() const
Returns the direction of the trajectory at the first point.
Definition: TrackTrajectory.h:296
recob::TrackTrajectory::Vertex
Point_t const & Vertex() const
Returns the position of the first valid point of the trajectory [cm].
Definition: TrackTrajectory.h:233
recob::TrackTrajectory::EndDirection
Vector_t EndDirection() const
Returns the direction of the trajectory at the last point.
Definition: TrackTrajectory.h:302
SMatrixSym55
recob::tracking::SMatrixSym55 SMatrixSym55
Definition: PMAlgTracking.cxx:22
recob::TrackTrajectory::End
Point_t const & End() const
Returns the position of the last valid point of the trajectory [cm].
Definition: TrackTrajectory.h:239
trkmkr::OptionalOutputs
Struct holding optional TrackMaker outputs.
Definition: TrackMaker.h:108
recob::Track
Track from a non-cascading particle.A recob::Track consists of a recob::TrackTrajectory, plus additional members relevant for a "fitted" track:
Definition: Track.h:49
bool trkf::TrackKalmanFitter::fitTrack ( detinfo::DetectorPropertiesData const &  detProp,
const Point_t &  position,
const Vector_t &  direction,
SMatrixSym55 &  trackStateCov,
const std::vector< art::Ptr< recob::Hit >> &  hits,
const std::vector< recob::TrajectoryPointFlags > &  flags,
const int  tkID,
const double  pval,
const int  pdgid,
recob::Track &  outTrack,
std::vector< art::Ptr< recob::Hit >> &  outHits,
trkmkr::OptionalOutputs &  optionals 
) const

Fit track starting from intial position, direction, and flags.

Definition at line 136 of file TrackKalmanFitter.cxx.

References cleanZigzag_, doFitWork(), dumpLevel_, fillResult(), hits(), setupInitialTrackState(), setupInputStates(), skipNegProp_, and tryNoSkipWhenFails_.

148 {
149  if (dumpLevel_ > 1)
150  std::cout << "Fitting track with tkID=" << tkID << " start pos=" << position
151  << " dir=" << direction << " nHits=" << hits.size() << " mom=" << pval
152  << " pdg=" << pdgid << std::endl;
153  if (hits.size() < 4) {
154  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__;
155  return false;
156  }
157 
158  // setup the KFTrackState we'll use throughout the fit
159  KFTrackState trackState = setupInitialTrackState(position, direction, trackStateCov, pval, pdgid);
160 
161  // setup vector of HitStates and flags, with either same or inverse order as input hit vector
162  // this is what we'll loop over during the fit
163  std::vector<HitState> hitstatev;
164  std::vector<recob::TrajectoryPointFlags::Mask_t> hitflagsv;
165  bool inputok = setupInputStates(detProp, hits, flags, hitstatev, hitflagsv);
166  if (!inputok) return false;
167 
168  // track and index vectors we use to store the fit results
169  std::vector<KFTrackState> fwdPrdTkState;
170  std::vector<KFTrackState> fwdUpdTkState;
171  std::vector<unsigned int> hitstateidx;
172  std::vector<unsigned int> rejectedhsidx;
173  std::vector<unsigned int> sortedtksidx;
174 
175  // do the actual fit
176  bool fitok = doFitWork(trackState,
177  detProp,
178  hitstatev,
179  hitflagsv,
180  fwdPrdTkState,
181  fwdUpdTkState,
182  hitstateidx,
183  rejectedhsidx,
184  sortedtksidx);
185  if (!fitok && (skipNegProp_ || cleanZigzag_) && tryNoSkipWhenFails_) {
186  mf::LogWarning("TrackKalmanFitter")
187  << "Trying to recover with skipNegProp = false and cleanZigzag = false\n";
188  fitok = doFitWork(trackState,
189  detProp,
190  hitstatev,
191  hitflagsv,
192  fwdPrdTkState,
193  fwdUpdTkState,
194  hitstateidx,
195  rejectedhsidx,
196  sortedtksidx,
197  false);
198  }
199  if (!fitok) {
200  mf::LogWarning("TrackKalmanFitter") << "Fit failed for track with ID=" << tkID << "\n";
201  return false;
202  }
203 
204  // fill the track, the output hit vector, and the optional output products
205  bool fillok = fillResult(hits,
206  tkID,
207  pdgid,
208  hitstatev,
209  hitflagsv,
210  fwdPrdTkState,
211  fwdUpdTkState,
212  hitstateidx,
213  rejectedhsidx,
214  sortedtksidx,
215  outTrack,
216  outHits,
217  optionals);
218  return fillok;
219 }
trkf::TrackKalmanFitter::doFitWork
bool doFitWork(KFTrackState &trackState, detinfo::DetectorPropertiesData const &detProp, std::vector< HitState > &hitstatev, std::vector< recob::TrajectoryPointFlags::Mask_t > &hitflagsv, std::vector< KFTrackState > &fwdPrdTkState, std::vector< KFTrackState > &fwdUpdTkState, std::vector< unsigned int > &hitstateidx, std::vector< unsigned int > &rejectedhsidx, std::vector< unsigned int > &sortedtksidx, bool applySkipClean=true) const
Function where the core of the fit is performed.
Definition: TrackKalmanFitter.cxx:295
trkf::TrackKalmanFitter::fillResult
bool fillResult(const std::vector< art::Ptr< recob::Hit >> &inHits, const int tkID, const int pdgid, std::vector< HitState > &hitstatev, std::vector< recob::TrajectoryPointFlags::Mask_t > &hitflagsv, std::vector< KFTrackState > &fwdPrdTkState, std::vector< KFTrackState > &fwdUpdTkState, std::vector< unsigned int > &hitstateidx, std::vector< unsigned int > &rejectedhsidx, std::vector< unsigned int > &sortedtksidx, recob::Track &outTrack, std::vector< art::Ptr< recob::Hit >> &outHits, trkmkr::OptionalOutputs &optionals) const
Fill the output objects.
Definition: TrackKalmanFitter.cxx:940
trkf::TrackKalmanFitter::setupInitialTrackState
KFTrackState setupInitialTrackState(const Point_t &position, const Vector_t &direction, SMatrixSym55 &trackStateCov, const double pval, const int pdgid) const
Return track state from intial position, direction, and covariance.
Definition: TrackKalmanFitter.cxx:221
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:210
trkf::TrackKalmanFitter::setupInputStates
bool setupInputStates(detinfo::DetectorPropertiesData const &detProp, const std::vector< art::Ptr< recob::Hit >> &hits, const std::vector< recob::TrajectoryPointFlags > &flags, std::vector< HitState > &hitstatev, std::vector< recob::TrajectoryPointFlags::Mask_t > &hitflagsv) const
Setup vectors of HitState and Masks to be used during the fit.
Definition: TrackKalmanFitter.cxx:246
trkf::KFTrackState trkf::TrackKalmanFitter::setupInitialTrackState ( const Point_t &  position,
const Vector_t &  direction,
SMatrixSym55 &  trackStateCov,
const double  pval,
const int  pdgid 
) const
private

Return track state from intial position, direction, and covariance.

Definition at line 221 of file TrackKalmanFitter.cxx.

Referenced by fitTrack().

226 {
227  //start from large enough covariance matrix so that the fit is not biased
228  if (trackStateCov == SMatrixSym55()) {
229  trackStateCov(0, 0) = 1000.;
230  trackStateCov(1, 1) = 1000.;
231  trackStateCov(2, 2) = 0.25;
232  trackStateCov(3, 3) = 0.25;
233  trackStateCov(4, 4) = 10.;
234  }
235  else
236  trackStateCov *= 100.;
237  // build vector of parameters on plane with point on the track and direction normal to the plane parallel to the track (so the first four parameters are zero by construction)
238  SVector5 trackStatePar(0., 0., 0., 0., 1. / pval);
239  return KFTrackState(trackStatePar,
240  trackStateCov,
241  Plane(position, direction),
242  true,
243  pdgid); //along direction by definition
244 }
trkf::SMatrixSym55
recob::tracking::SMatrixSym55 SMatrixSym55
Definition: TrackState.h:15
trkf::SVector5
recob::tracking::SVector5 SVector5
Definition: TrackState.h:12
trkf::Plane
recob::tracking::Plane Plane
Definition: TrackState.h:17
bool trkf::TrackKalmanFitter::setupInputStates ( detinfo::DetectorPropertiesData const &  detProp,
const std::vector< art::Ptr< recob::Hit >> &  hits,
const std::vector< recob::TrajectoryPointFlags > &  flags,
std::vector< HitState > &  hitstatev,
std::vector< recob::TrajectoryPointFlags::Mask_t > &  hitflagsv 
) const
private

Setup vectors of HitState and Masks to be used during the fit.

Definition at line 246 of file TrackKalmanFitter.cxx.

References detinfo::DetectorPropertiesData::ConvertTicksToX(), recob::TrajectoryPointFlags::DefaultFlagsMask(), dumpLevel_, recob::TrajectoryPointFlagTraits::ExcludedFromFit, geom, detinfo::DetectorPropertiesData::GetXTicksCoefficient(), hitErr2ScaleFact_, hits(), recob::TrajectoryPointFlagTraits::Merged, recob::TrajectoryPointFlagTraits::NoPoint, rejectHighMultHits_, rejectHitsNegativeGOF_, recob::TrajectoryPointFlagTraits::Suspicious, useRMS_, geo::GeometryCore::WireIDToWireGeo(), and x.

Referenced by fitTrack().

252 {
253  if (dumpLevel_ > 1)
254  std::cout << "flags.size()=" << flags.size() << " hits.size()=" << hits.size() << std::endl;
255 
256  // setup vector of HitStates and flags, with either same or inverse order as input hit vector
257  // this is what we'll loop over during the fit
258  const size_t fsize = flags.size();
259  for (size_t ihit = 0; ihit != hits.size(); ihit++) {
260  const auto& hit = hits[ihit];
261  double t = hit->PeakTime();
262  double terr = (useRMS_ ? hit->RMS() : hit->SigmaPeakTime());
263  double x =
264  detProp.ConvertTicksToX(t, hit->WireID().Plane, hit->WireID().TPC, hit->WireID().Cryostat);
265  double xerr = terr * detProp.GetXTicksCoefficient();
266  hitstatev.emplace_back(
267  x, hitErr2ScaleFact_ * xerr * xerr, hit->WireID(), geom->WireIDToWireGeo(hit->WireID()));
268  //
269  if (fsize > 0 && ihit < fsize)
270  hitflagsv.push_back(flags[ihit].mask());
271  else
272  hitflagsv.push_back(recob::TrajectoryPointFlags::DefaultFlagsMask());
273  //
274  if (rejectHighMultHits_ && hit->Multiplicity() > 1) {
275  hitflagsv.back().set(recob::TrajectoryPointFlagTraits::Merged);
276  hitflagsv.back().set(recob::TrajectoryPointFlagTraits::ExcludedFromFit);
277  }
278  if (rejectHitsNegativeGOF_ && hit->GoodnessOfFit() < 0) {
279  hitflagsv.back().set(recob::TrajectoryPointFlagTraits::Suspicious);
280  hitflagsv.back().set(recob::TrajectoryPointFlagTraits::ExcludedFromFit);
281  }
282  //
283  if (dumpLevel_ > 2)
284  std::cout << "pushed flag mask=" << hitflagsv.back()
285  << " merged=" << hitflagsv.back().isSet(recob::TrajectoryPointFlagTraits::Merged)
286  << " suspicious="
287  << hitflagsv.back().isSet(recob::TrajectoryPointFlagTraits::Suspicious)
288  << " nopoint=" << hitflagsv.back().isSet(recob::TrajectoryPointFlagTraits::NoPoint)
289  << std::endl;
290  }
291  if (dumpLevel_ > 2) assert(hits.size() == hitstatev.size());
292  return true;
293 }
x
Float_t x
Definition: compare.C:6
recob::TrajectoryPointFlagTraits::Merged
static constexpr Flag_t Merged
The hit might have contribution from particles other than this.
Definition: TrajectoryPointFlags.h:131
geo::GeometryCore::WireIDToWireGeo
WireGeo const & WireIDToWireGeo(WireID const &wireid) const
Returns the specified wire.
Definition: GeometryCore.h:1211
recob::TrajectoryPointFlagTraits::Suspicious
static constexpr Flag_t Suspicious
The point reconstruction is somehow questionable.
Definition: TrajectoryPointFlags.h:128
recob::TrajectoryPointFlagTraits::NoPoint
static constexpr Flag_t NoPoint
The trajectory point is not defined.
Definition: TrajectoryPointFlags.h:118
recob::TrajectoryPointFlags::DefaultFlagsMask
static constexpr Mask_t DefaultFlagsMask()
Flags used in default construction.
Definition: TrajectoryPointFlags.h:731
trkf::TrackKalmanFitter::geom
art::ServiceHandle< geo::Geometry const > geom
Definition: TrackKalmanFitter.h:298
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
recob::TrajectoryPointFlagTraits::ExcludedFromFit
static constexpr Flag_t ExcludedFromFit
Definition: TrajectoryPointFlags.h:168
void trkf::TrackKalmanFitter::sortOutput ( std::vector< HitState > &  hitstatev,
std::vector< KFTrackState > &  fwdUpdTkState,
std::vector< unsigned int > &  hitstateidx,
std::vector< unsigned int > &  rejectedhsidx,
std::vector< unsigned int > &  sortedtksidx,
std::vector< recob::TrajectoryPointFlags::Mask_t > &  hitflagsv,
bool  applySkipClean = true 
) const
private

Sort the output states.

Definition at line 815 of file TrackKalmanFitter.cxx.

References cleanZigzag_, dir, dumpLevel_, util::end(), recob::TrajectoryPointFlagTraits::ExcludedFromFit, geom, geo::GeometryCore::MaxPlanes(), negDistTolerance_, skipNegProp_, and sortOutputHitsMinLength_.

Referenced by doFitWork().

823 {
824  //
825  if (sortOutputHitsMinLength_) {
826  //sort hits keeping fixed the order on planes and picking the closest next plane
827  const unsigned int nplanes = geom->MaxPlanes();
828  std::vector<std::vector<unsigned int>> tracksInPlanes(nplanes);
829  for (unsigned int p = 0; p < hitstateidx.size(); ++p) {
830  const auto& hitstate = hitstatev[hitstateidx[p]];
831  tracksInPlanes[hitstate.wireId().Plane].push_back(p);
832  }
833  if (dumpLevel_ > 2) {
834  for (auto s : fwdUpdTkState) {
835  std::cout << "state pos=" << s.position() << std::endl;
836  }
837  }
838  //find good starting point
839  std::vector<unsigned int> iterTracksInPlanes;
840  for (auto it : tracksInPlanes)
841  iterTracksInPlanes.push_back(0);
842  auto pos = fwdUpdTkState.front().position();
843  auto dir = fwdUpdTkState.front().momentum();
844  unsigned int p = 0;
845  for (; p < fwdUpdTkState.size(); ++p) {
846  if (hitflagsv[hitstateidx[p]].isSet(recob::TrajectoryPointFlagTraits::ExcludedFromFit) == 0) {
847  pos = fwdUpdTkState[p].position();
848  dir = fwdUpdTkState[p].momentum();
849  if (dumpLevel_ > 2)
850  std::cout << "sort output found point not excluded with p=" << p
851  << " hitstateidx[p]=" << hitstateidx[p] << " pos=" << pos << std::endl;
852  break;
853  }
854  else {
855  rejectedhsidx.push_back(hitstateidx[p]);
856  }
857  }
858  if (dumpLevel_ > 1)
859  std::cout << "sort output init with pos=" << pos << " dir=" << dir << std::endl;
860  //pick hits based on minimum dot product with respect to position and direction at previous hit
861  for (; p < fwdUpdTkState.size(); ++p) {
862  int min_plane = -1;
863  double min_dotp = DBL_MAX;
864  for (unsigned int iplane = 0; iplane < iterTracksInPlanes.size(); ++iplane) {
865  for (unsigned int& itk = iterTracksInPlanes[iplane]; itk < tracksInPlanes[iplane].size();
866  ++itk) {
867  auto& trackstate = fwdUpdTkState[tracksInPlanes[iplane][iterTracksInPlanes[iplane]]];
868  auto& tmppos = trackstate.position();
869  const double dotp = dir.Dot(tmppos - pos);
870  if (dumpLevel_ > 2)
871  std::cout << "iplane=" << iplane << " tmppos=" << tmppos << " tmpdir=" << tmppos - pos
872  << " dotp=" << dotp << std::endl;
873  if (dotp < min_dotp) {
874  min_plane = iplane;
875  min_dotp = dotp;
876  }
877  break;
878  }
879  }
880  if (min_plane < 0) continue;
881  const unsigned int ihit = tracksInPlanes[min_plane][iterTracksInPlanes[min_plane]];
882  if (applySkipClean && skipNegProp_ && min_dotp < negDistTolerance_ &&
883  sortedtksidx.size() > 0) {
884  if (dumpLevel_ > 2)
885  std::cout << "sort output rejecting hit #" << ihit << " plane=" << min_plane
886  << " with min_dotp=" << min_dotp << std::endl;
887  rejectedhsidx.push_back(hitstateidx[ihit]);
888  iterTracksInPlanes[min_plane]++;
889  continue;
890  }
891  if (dumpLevel_ > 2)
892  std::cout << "sort output picking hit #" << ihit << " plane=" << min_plane
893  << " with min_dotp=" << min_dotp << std::endl;
894  auto& trackstate = fwdUpdTkState[ihit];
895  pos = trackstate.position();
896  dir = trackstate.momentum();
897  //
898  sortedtksidx.push_back(ihit);
899  iterTracksInPlanes[min_plane]++;
900  }
901  }
902  else {
903  for (unsigned int p = 0; p < fwdUpdTkState.size(); ++p) {
904  sortedtksidx.push_back(p);
905  }
906  }
907  //
908  if (applySkipClean && cleanZigzag_) {
909  std::vector<unsigned int> itoerase;
910  bool clean = false;
911  while (!clean) {
912  bool broken = false;
913  auto pos0 = fwdUpdTkState[sortedtksidx[0]].position();
914  unsigned int i = 1;
915  unsigned int end = sortedtksidx.size() - 1;
916  for (; i < end; ++i) {
917  auto dir0 = fwdUpdTkState[sortedtksidx[i]].position() - pos0;
918  auto dir2 =
919  fwdUpdTkState[sortedtksidx[i + 1]].position() - fwdUpdTkState[sortedtksidx[i]].position();
920  dir0 /= dir0.R();
921  dir2 /= dir2.R();
922  if (dir2.Dot(dir0) < 0.) {
923  broken = true;
924  end--;
925  break;
926  }
927  else
928  pos0 = fwdUpdTkState[sortedtksidx[i]].position();
929  }
930  if (!broken) { clean = true; }
931  else {
932  rejectedhsidx.push_back(hitstateidx[sortedtksidx[i]]);
933  sortedtksidx.erase(sortedtksidx.begin() + i);
934  }
935  }
936  }
937  //
938 }
trkf::TrackKalmanFitter::geom
art::ServiceHandle< geo::Geometry const > geom
Definition: TrackKalmanFitter.h:298
util::end
decltype(auto) constexpr end(T &&obj)
ADL-aware version of std::end.
Definition: StdUtils.h:77
geo::GeometryCore::MaxPlanes
unsigned int MaxPlanes() const
Returns the largest number of planes among all TPCs in this detector.
Definition: GeometryCore.cxx:818
recob::TrajectoryPointFlagTraits::ExcludedFromFit
static constexpr Flag_t ExcludedFromFit
Definition: TrajectoryPointFlags.h:168
dir
TDirectory * dir
Definition: macro.C:5

Member Data Documentation

bool trkf::TrackKalmanFitter::cleanZigzag_
private

Definition at line 305 of file TrackKalmanFitter.h.

Referenced by fitTrack(), and sortOutput().

int trkf::TrackKalmanFitter::dumpLevel_
private

Definition at line 317 of file TrackKalmanFitter.h.

Referenced by doFitWork(), fillResult(), fitTrack(), setupInputStates(), and sortOutput().

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

Definition at line 298 of file TrackKalmanFitter.h.

Referenced by doFitWork(), setupInputStates(), and sortOutput().

float trkf::TrackKalmanFitter::hitErr2ScaleFact_
private

Definition at line 308 of file TrackKalmanFitter.h.

Referenced by setupInputStates().

float trkf::TrackKalmanFitter::maxChi2_
private

Definition at line 314 of file TrackKalmanFitter.h.

Referenced by doFitWork().

float trkf::TrackKalmanFitter::maxDist_
private

Definition at line 315 of file TrackKalmanFitter.h.

Referenced by doFitWork().

float trkf::TrackKalmanFitter::maxResidue_
private

Definition at line 312 of file TrackKalmanFitter.h.

Referenced by doFitWork().

float trkf::TrackKalmanFitter::maxResidueFirstHit_
private

Definition at line 313 of file TrackKalmanFitter.h.

Referenced by doFitWork().

float trkf::TrackKalmanFitter::negDistTolerance_
private

Definition at line 316 of file TrackKalmanFitter.h.

Referenced by doFitWork(), and sortOutput().

bool trkf::TrackKalmanFitter::pickBestHitOnWire_
private

Definition at line 311 of file TrackKalmanFitter.h.

Referenced by doFitWork().

const TrackStatePropagator* trkf::TrackKalmanFitter::propagator
private

Definition at line 299 of file TrackKalmanFitter.h.

Referenced by doFitWork(), and fillResult().

bool trkf::TrackKalmanFitter::rejectHighMultHits_
private

Definition at line 306 of file TrackKalmanFitter.h.

Referenced by setupInputStates().

bool trkf::TrackKalmanFitter::rejectHitsNegativeGOF_
private

Definition at line 307 of file TrackKalmanFitter.h.

Referenced by setupInputStates().

bool trkf::TrackKalmanFitter::skipNegProp_
private

Definition at line 304 of file TrackKalmanFitter.h.

Referenced by doFitWork(), fitTrack(), and sortOutput().

bool trkf::TrackKalmanFitter::sortHitsByPlane_
private

Definition at line 301 of file TrackKalmanFitter.h.

Referenced by doFitWork().

bool trkf::TrackKalmanFitter::sortHitsByWire_
private

Definition at line 302 of file TrackKalmanFitter.h.

Referenced by doFitWork().

bool trkf::TrackKalmanFitter::sortOutputHitsMinLength_
private

Definition at line 303 of file TrackKalmanFitter.h.

Referenced by sortOutput().

bool trkf::TrackKalmanFitter::tryBothDirs_
private

Definition at line 310 of file TrackKalmanFitter.h.

Referenced by fitTrack().

bool trkf::TrackKalmanFitter::tryNoSkipWhenFails_
private

Definition at line 309 of file TrackKalmanFitter.h.

Referenced by fitTrack().

bool trkf::TrackKalmanFitter::useRMS_
private

Definition at line 300 of file TrackKalmanFitter.h.

Referenced by setupInputStates().

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