LArSoft  v07_13_02
Liquid Argon Software toolkit - http://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 sortOutputHitsMinLength, bool skipNegProp, bool cleanZigzag, bool rejectHighMultHits, bool rejectHitsNegativeGOF, float hitErr2ScaleFact, bool tryNoSkipWhenFails, 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 (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 (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, 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 (const std::vector< art::Ptr< recob::Hit > > &hits, const std::vector< recob::TrajectoryPointFlags > &flags, const KFTrackState &trackState, bool &reverseHits, 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, 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, const bool reverseHits, 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 > geom
 
const detinfo::DetectorProperties * detprop
 
const TrackStatePropagator * propagator
 
bool useRMS_
 
bool sortHitsByPlane_
 
bool sortOutputHitsMinLength_
 
bool skipNegProp_
 
bool cleanZigzag_
 
bool rejectHighMultHits_
 
bool rejectHitsNegativeGOF_
 
float hitErr2ScaleFact_
 
bool tryNoSkipWhenFails_
 
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 46 of file TrackKalmanFitter.h.

Member Typedef Documentation

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

Definition at line 104 of file TrackKalmanFitter.h.

Constructor & Destructor Documentation

trkf::TrackKalmanFitter::TrackKalmanFitter ( const TrackStatePropagator *  prop,
bool  useRMS,
bool  sortHitsByPlane,
bool  sortOutputHitsMinLength,
bool  skipNegProp,
bool  cleanZigzag,
bool  rejectHighMultHits,
bool  rejectHitsNegativeGOF,
float  hitErr2ScaleFact,
bool  tryNoSkipWhenFails,
int  dumpLevel 
)
inline

Constructor from TrackStatePropagator and values of configuration parameters.

Definition at line 107 of file TrackKalmanFitter.h.

108  {
109  propagator=prop;
110  useRMS_=useRMS;
111  sortHitsByPlane_=sortHitsByPlane;
112  sortOutputHitsMinLength_=sortOutputHitsMinLength;
113  skipNegProp_=skipNegProp;
114  cleanZigzag_=cleanZigzag;
115  rejectHighMultHits_=rejectHighMultHits;
116  rejectHitsNegativeGOF_=rejectHitsNegativeGOF;
117  hitErr2ScaleFact_=hitErr2ScaleFact;
118  tryNoSkipWhenFails_=tryNoSkipWhenFails;
119  dumpLevel_=dumpLevel;
120  detprop = art::ServiceHandle<detinfo::DetectorPropertiesService>()->provider();
121  }
trkf::TrackKalmanFitter::detprop
const detinfo::DetectorProperties * detprop
Definition: TrackKalmanFitter.h:167
trkf::TrackKalmanFitter::propagator
const TrackStatePropagator * propagator
Definition: TrackKalmanFitter.h:168
trkf::TrackKalmanFitter::TrackKalmanFitter ( const TrackStatePropagator *  prop,
Parameters const &  p 
)
inlineexplicit

Constructor from TrackStatePropagator and Parameters table.

Definition at line 124 of file TrackKalmanFitter.h.

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

125  : TrackKalmanFitter(prop,p().useRMS(),p().sortHitsByPlane(),p().sortOutputHitsMinLength(),p().skipNegProp(),p().cleanZigzag(),
126  p().rejectHighMultHits(),p().rejectHitsNegativeGOF(),p().hitErr2ScaleFact(),p().tryNoSkipWhenFails(),p().dumpLevel()) {}
trkf::TrackKalmanFitter::TrackKalmanFitter
TrackKalmanFitter(const TrackStatePropagator *prop, bool useRMS, bool sortHitsByPlane, bool sortOutputHitsMinLength, bool skipNegProp, bool cleanZigzag, bool rejectHighMultHits, bool rejectHitsNegativeGOF, float hitErr2ScaleFact, bool tryNoSkipWhenFails, int dumpLevel)
Constructor from TrackStatePropagator and values of configuration parameters.
Definition: TrackKalmanFitter.h:107

Member Function Documentation

bool trkf::TrackKalmanFitter::doFitWork ( KFTrackState &  trackState,
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 135 of file TrackKalmanFitter.cxx.

References trkf::TrackStatePropagator::BACKWARD, trkf::KFTrackState::covariance(), recob::TrajectoryPointFlagTraits::DeltaRay, recob::TrajectoryPointFlagTraits::DetectorIssue, trkf::TrackStatePropagator::distanceToPlane(), trkf::KFTrackState::dump(), trkf::TrackState::dump(), dumpLevel_, recob::TrajectoryPointFlagTraits::ExcludedFromFit, trkf::TrackStatePropagator::FORWARD, geom, recob::TrajectoryPointFlagTraits::HitIgnored, geo::GeometryCore::MaxPlanes(), recob::TrajectoryPointFlagTraits::Merged, min, recob::TrajectoryPointFlagTraits::NoPoint, trkf::TrackState::plane(), trkf::KFTrackState::position(), trkf::TrackStatePropagator::propagateToPlane(), propagator, recob::TrajectoryPointFlagTraits::Rejected, trkf::KFTrackState::setCovariance(), recob::TrajectoryPointFlagTraits::Shared, skipNegProp_, sortHitsByPlane_, sortOutput(), recob::TrajectoryPointFlagTraits::Suspicious, trkf::KFTrackState::trackState(), and trkf::KFTrackState::updateWithHitState().

Referenced by fitTrack().

138  {
139 
140  fwdPrdTkState.clear();
141  fwdUpdTkState.clear();
142  hitstateidx.clear();
143  rejectedhsidx.clear();
144  sortedtksidx.clear();
145  // these three vectors are aligned
146  fwdPrdTkState.reserve(hitstatev.size());
147  fwdUpdTkState.reserve(hitstatev.size());
148  hitstateidx.reserve(hitstatev.size());
149 
150  if (sortHitsByPlane_) {
151  //array of hit indices in planes, keeping the original sorting by plane
152  const unsigned int nplanes = geom->MaxPlanes();
153  std::vector<std::vector<unsigned int> > hitsInPlanes(nplanes);
154  for (unsigned int ihit = 0; ihit<hitstatev.size(); ihit++) {
155  hitsInPlanes[hitstatev[ihit].wireId().Plane].push_back(ihit);
156  }
157  //array of indices, where iterHitsInPlanes[i] is the iterator over hitsInPlanes[i]
158  std::vector<unsigned int> iterHitsInPlanes(nplanes,0);
159  for (unsigned int p = 0; p<hitstatev.size(); ++p) {
160  if (dumpLevel_>1) std::cout << std::endl << "processing hit #" << p << std::endl;
161  if (dumpLevel_>1) {
162  std::cout << "compute distance from state=" << std::endl; trackState.dump();
163  }
164  int min_plane = -1;
165  double min_dist = DBL_MAX;
166  //find the closest hit according to the sorting in each plane
167  for (unsigned int iplane = 0; iplane<nplanes; ++iplane) {
168  //note: ih is a reference, so when 'continue' we increment iterHitsInPlanes[iplane] and the hit is effectively rejected
169  for (unsigned int& ih = iterHitsInPlanes[iplane]; ih<hitsInPlanes[iplane].size(); ++ih) {
170  if (dumpLevel_>1) std::cout << "iplane=" << iplane << " nplanes=" << nplanes << " iterHitsInPlanes[iplane]=" << iterHitsInPlanes[iplane] << " hitsInPlanes[iplane].size()=" << hitsInPlanes[iplane].size() << std::endl;
171  unsigned int ihit = hitsInPlanes[iplane][ih];
172  const auto& hitstate = hitstatev[ihit];
173  const auto& hitflags = hitflagsv[ihit];
174  if (hitflags.isSet(recob::TrajectoryPointFlagTraits::NoPoint ) ||
175  hitflags.isSet(recob::TrajectoryPointFlagTraits::ExcludedFromFit) ||
176  hitflags.isSet(recob::TrajectoryPointFlagTraits::Rejected )) {
177  if (dumpLevel_>1) std::cout << "rejecting hit flags=" << hitflags.isSet(recob::TrajectoryPointFlagTraits::NoPoint ) << ", " << hitflags.isSet(recob::TrajectoryPointFlagTraits::ExcludedFromFit) << ", " << hitflags.isSet(recob::TrajectoryPointFlagTraits::Rejected ) << std::endl;
178  rejectedhsidx.push_back(ihit);
179  continue;
180  }
181  //get distance to measurement surface
182  bool success = true;
183  const double dist = propagator->distanceToPlane(success, trackState.trackState(), hitstate.plane());
184  if (dumpLevel_>1) std::cout << "distance to plane " << iplane << " wire=" << hitstate.wireId().Wire << " = " << dist << ", min_dist=" << std::min(min_dist,999.) << " min_plane=" << min_plane << " success=" << success << std::endl;
185  if (!success) {
186  rejectedhsidx.push_back(ihit);
187  continue;
188  }
189  if (applySkipClean && skipNegProp_ && dist<0.) {
190  rejectedhsidx.push_back(ihit);
191  continue;
192  }
193  if (dist<min_dist) {
194  min_plane = iplane;
195  min_dist = dist;
196  }
197  break;
198  }
199  }
200  if (dumpLevel_>1) std::cout << "pick min_dist=" << min_dist << " min_plane=" << min_plane << std::endl;
201  //now we know which is the closest wire: get the hitstate and increment the iterator
202  if (min_plane<0) continue;
203  const unsigned int ihit = hitsInPlanes[min_plane][iterHitsInPlanes[min_plane]];
204  const auto& hitstate = hitstatev[ihit];
205  if (dumpLevel_>1) hitstate.dump();
206  auto& hitflags = hitflagsv[ihit];
207  iterHitsInPlanes[min_plane]++;
208  //propagate to measurement surface
209  bool propok = true;
210  trackState = propagator->propagateToPlane(propok, trackState.trackState(), hitstate.plane(), true, true, TrackStatePropagator::FORWARD);
211  if (!propok && !(applySkipClean && skipNegProp_)) {
212  if (dumpLevel_>1) std::cout << "attempt backward prop" << std::endl;
213  trackState = propagator->propagateToPlane(propok, trackState.trackState(), hitstate.plane(), true, true, TrackStatePropagator::BACKWARD);
214  }
215  if (dumpLevel_>1) {
216  std::cout << "hit state " << std::endl; hitstate.dump();
217  std::cout << "propagation result=" << propok << std::endl;
218  std::cout << "propagated state " << std::endl; trackState.dump();
219  std::cout << "propagated planarity=" << hitstate.plane().direction().Dot(hitstate.plane().position()-trackState.position()) << std::endl;
220  }
221  if (propok) {
222  hitstateidx.push_back(ihit);
223  fwdPrdTkState.push_back(trackState);
224  //
225  if (hitflags.isSet(recob::TrajectoryPointFlagTraits::HitIgnored ) ||
226  hitflags.isSet(recob::TrajectoryPointFlagTraits::Merged ) ||
227  hitflags.isSet(recob::TrajectoryPointFlagTraits::Shared ) ||
228  hitflags.isSet(recob::TrajectoryPointFlagTraits::DeltaRay ) ||
229  hitflags.isSet(recob::TrajectoryPointFlagTraits::DetectorIssue) ||
230  hitflags.isSet(recob::TrajectoryPointFlagTraits::Suspicious )) {
231  //do not update the hit, mark as excluded from the fit
232  fwdUpdTkState.push_back(trackState);
233  hitflags.set(recob::TrajectoryPointFlagTraits::ExcludedFromFit);
234  continue;
235  }
236  //now update the forward fitted track
237  bool upok = trackState.updateWithHitState(hitstate);
238  if (upok==0) {
239  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__;
240  return false;
241  }
242  if (dumpLevel_>1) {
243  std::cout << "updated state " << std::endl; trackState.dump();
244  }
245  fwdUpdTkState.push_back(trackState);
246  } else {
247  if (dumpLevel_>0) std::cout << "WARNING: forward propagation failed. Skip this hit..." << std::endl;
248  rejectedhsidx.push_back(ihit);
249  continue;
250  }
251  }
252  } else {
253  for (unsigned int ihit=0; ihit<hitstatev.size(); ++ihit) {
254  const auto& hitstate = hitstatev[ihit];
255  if (dumpLevel_>1) {
256  std::cout << std::endl << "processing hit #" << ihit << std::endl;
257  hitstate.dump();
258  }
259  auto& hitflags = hitflagsv[ihit];
260  if (hitflags.isSet(recob::TrajectoryPointFlagTraits::ExcludedFromFit) ||
261  hitflags.isSet(recob::TrajectoryPointFlagTraits::Rejected)) {
262  rejectedhsidx.push_back(ihit);
263  continue;
264  }
265  if (applySkipClean && skipNegProp_) {
266  bool success = true;
267  const double dist = propagator->distanceToPlane(success, trackState.trackState(), hitstate.plane());
268  if (dist<0. || success==false) {
269  if (dumpLevel_>0) std::cout << "WARNING: negative propagation distance. Skip this hit..." << std::endl;;
270  rejectedhsidx.push_back(ihit);
271  continue;
272  }
273  }
274  //propagate to measurement surface
275  bool propok = true;
276  trackState = propagator->propagateToPlane(propok, trackState.trackState(), hitstate.plane(), true, true, TrackStatePropagator::FORWARD);
277  if (!propok && !(applySkipClean && skipNegProp_)) trackState = propagator->propagateToPlane(propok, trackState.trackState(), hitstate.plane(), true, true, TrackStatePropagator::BACKWARD);
278  if (propok) {
279  hitstateidx.push_back(ihit);
280  fwdPrdTkState.push_back(trackState);
281  //
282  if (hitflags.isSet(recob::TrajectoryPointFlagTraits::HitIgnored ) ||
283  hitflags.isSet(recob::TrajectoryPointFlagTraits::Merged ) ||
284  hitflags.isSet(recob::TrajectoryPointFlagTraits::Shared ) ||
285  hitflags.isSet(recob::TrajectoryPointFlagTraits::DeltaRay ) ||
286  hitflags.isSet(recob::TrajectoryPointFlagTraits::DetectorIssue) ||
287  hitflags.isSet(recob::TrajectoryPointFlagTraits::Suspicious )) {
288  //do not update the hit, mark as excluded from the fit
289  fwdUpdTkState.push_back(trackState);
290  hitflags.set(recob::TrajectoryPointFlagTraits::ExcludedFromFit);
291  continue;
292  }
293  //
294  bool upok = trackState.updateWithHitState(hitstate);
295  if (upok==0) {
296  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__;
297  return false;
298  }
299  fwdUpdTkState.push_back(trackState);
300  } else {
301  if (dumpLevel_>0) std::cout << "WARNING: forward propagation failed. Skip this hit..." << std::endl;;
302  rejectedhsidx.push_back(ihit);
303  continue;
304  }
305  }//for (auto hitstate : hitstatev)
306  }
307 
308  if (dumpLevel_>2) assert( rejectedhsidx.size()+hitstateidx.size() == hitstatev.size() );
309  if (dumpLevel_>0) {
310  std::cout << "TRACK AFTER FWD" << std::endl;
311  trackState.dump();
312  }
313 
314  //reinitialize trf for backward fit, scale the error to avoid biasing the backward fit
315  trackState.setCovariance(100.*trackState.covariance());
316 
317  //backward loop over track states and hits in fwdUpdTracks: use hits for backward fit and fwd track states for smoothing
318  int nvalidhits = 0;
319  for (int itk = fwdPrdTkState.size()-1; itk>=0; itk--) {
320  auto& fwdPrdTrackState = fwdPrdTkState[itk];
321  auto& fwdUpdTrackState = fwdUpdTkState[itk];
322  const auto& hitstate = hitstatev[hitstateidx[itk]];
323  auto& hitflags = hitflagsv[hitstateidx[itk]];
324  if (dumpLevel_>1) {
325  std::cout << std::endl << "processing backward hit #" << itk << std::endl;
326  hitstate.dump();
327  }
328  bool propok = true;
329  trackState = propagator->propagateToPlane(propok, trackState.trackState(), hitstate.plane(), true, true, TrackStatePropagator::BACKWARD);
330  if (!propok) trackState = propagator->propagateToPlane(propok, trackState.trackState(), hitstate.plane(), true, true, TrackStatePropagator::FORWARD);//do we want this?
331  //
332  if (dumpLevel_>1) {
333  std::cout << "propagation result=" << propok << std::endl;
334  std::cout << "propagated state " << std::endl; trackState.dump();
335  std::cout << "propagated planarity=" << hitstate.plane().direction().Dot(hitstate.plane().position()-trackState.position()) << std::endl;
336  }
337  if (propok) {
338  //combine forward predicted and backward predicted
339  bool pcombok = fwdPrdTrackState.combineWithTrackState(trackState.trackState());
340  if (pcombok==0) {
341  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__;
342  return false;
343  }
344  if (dumpLevel_>1) {
345  std::cout << "combined prd state " << std::endl; fwdPrdTrackState.dump();
346  }
347  //combine forward updated and backward predicted
348  bool ucombok = fwdUpdTrackState.combineWithTrackState(trackState.trackState());
349  if (ucombok==0) {
350  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__;
351  return false;
352  }
353  if (dumpLevel_>1) {
354  std::cout << "combined upd state " << std::endl; fwdUpdTrackState.dump();
355  }
356  //update backward predicted, only if the hit was not excluded
357  if (hitflags.isSet(recob::TrajectoryPointFlagTraits::ExcludedFromFit)==0) {
358  bool upok = trackState.updateWithHitState(hitstate);
359  if (upok==0) {
360  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__;
361  return false;
362  }
363  if (dumpLevel_>1) {
364  std::cout << "updated state " << std::endl; trackState.dump();
365  }
366  //
367  nvalidhits++;
368  }
369  } else {
370  // ok, if the backward propagation failed we exclude this point from the rest of the fit,
371  // but we can still use its position from the forward fit, so just mark it as ExcludedFromFit
372  hitflags.set(recob::TrajectoryPointFlagTraits::ExcludedFromFit);
373  if (dumpLevel_>0) std::cout << "WARNING: backward propagation failed. Skip this hit..." << std::endl;;
374  continue;
375  }
376  }//for (unsigned int itk = fwdPrdTkState.size()-1; itk>=0; itk--) {
377 
378  if (nvalidhits<4) {
379  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__ << " ";
380  return false;
381  }
382 
383  // sort output states
384  sortOutput(hitstatev, fwdUpdTkState, hitstateidx, rejectedhsidx, sortedtksidx, applySkipClean);
385  if (sortedtksidx.size()<4) {
386  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__ << " ";
387  return false;
388  }
389 
390  if (dumpLevel_>2) assert( rejectedhsidx.size()+sortedtksidx.size() == hitstatev.size() );
391  return true;
392 }
recob::TrajectoryPointFlagTraits::Merged
static constexpr Flag_t Merged
The hit might have contribution from particles other than this.
Definition: TrajectoryPointFlags.h:135
recob::TrajectoryPointFlagTraits::Suspicious
static constexpr Flag_t Suspicious
The point reconstruction is somehow questionable.
Definition: TrajectoryPointFlags.h:132
recob::TrajectoryPointFlagTraits::NoPoint
static constexpr Flag_t NoPoint
The trajectory point is not defined.
Definition: TrajectoryPointFlags.h:122
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:178
geo::GeometryCore::MaxPlanes
unsigned int MaxPlanes() const
Returns the largest number of planes among all TPCs in this detector.
Definition: GeometryCore.cxx:837
trkf::TrackKalmanFitter::propagator
const TrackStatePropagator * propagator
Definition: TrackKalmanFitter.h:168
trkf::TrackKalmanFitter::geom
art::ServiceHandle< geo::Geometry > geom
Definition: TrackKalmanFitter.h:166
recob::TrajectoryPointFlagTraits::HitIgnored
static constexpr Flag_t HitIgnored
Hit was not included for the computation of the trajectory.
Definition: TrajectoryPointFlags.h:119
recob::TrajectoryPointFlagTraits::Rejected
static constexpr Flag_t Rejected
The hit is extraneous to this track.
Definition: TrajectoryPointFlags.h:175
recob::TrajectoryPointFlagTraits::ExcludedFromFit
static constexpr Flag_t ExcludedFromFit
Definition: TrajectoryPointFlags.h:172
min
Int_t min
Definition: plot.C:26
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h: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, bool applySkipClean=true) const
Sort the output states.
Definition: TrackKalmanFitter.cxx:394
trkf::TrackState::dump
std::ostream & dump(std::ostream &out=std::cout) const
Printout information.
Definition: TrackState.h:119
trkf::TrackStatePropagator::propagateToPlane
TrackState propagateToPlane(bool &success, 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:27
recob::TrajectoryPointFlagTraits::DetectorIssue
static constexpr Flag_t DetectorIssue
The hit is associated to a problematic channel.
Definition: TrajectoryPointFlags.h:141
recob::TrajectoryPointFlagTraits::DeltaRay
static constexpr Flag_t DeltaRay
The hit might have contribution from a δ ray.
Definition: TrajectoryPointFlags.h:138
recob::TrajectoryPointFlagTraits::Shared
static constexpr Flag_t Shared
The hit is known to be associated also to another trajectory.
Definition: TrajectoryPointFlags.h:144
bool trkf::TrackKalmanFitter::fillResult ( const std::vector< art::Ptr< recob::Hit > > &  inHits,
const int  tkID,
const int  pdgid,
const bool  reverseHits,
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 477 of file TrackKalmanFitter.cxx.

References trkmkr::TrackCreationBookKeeper::addPoint(), 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().

481  {
482  // fill output trajectory objects with smoothed track and its hits
483  trkmkr::TrackCreationBookKeeper tcbk(outHits, optionals, tkID, pdgid, true);
484  for (unsigned int p : sortedtksidx) {
485  const auto& trackstate = fwdUpdTkState[p];
486  const auto& hitflags = hitflagsv[hitstateidx[p]];
487  const unsigned int originalPos = (reverseHits ? hitstatev.size()-hitstateidx[p]-1 : hitstateidx[p]);
488  if (dumpLevel_>2) assert(originalPos>=0 && originalPos<hitstatev.size());
489  //
490  const auto& prdtrack = fwdPrdTkState[p];
491  const auto& hitstate = hitstatev[hitstateidx[p]];
492  if (dumpLevel_>2) assert(hitstate.wireId().Plane == inHits[originalPos]->WireID().Plane);
493  //
494  trkmkr::OptionalPointElement ope;
495  if (optionals.isTrackFitInfosInit()) {
496  ope.setTrackFitHitInfo(recob::TrackFitHitInfo(hitstate.hitMeas(),hitstate.hitMeasErr2(),prdtrack.parameters(),prdtrack.covariance(),hitstate.wireId()));
497  }
498  tcbk.addPoint(trackstate.position(), trackstate.momentum(), inHits[originalPos],
499  recob::TrajectoryPointFlags(originalPos,hitflags), prdtrack.chi2(hitstate), ope);
500  }
501 
502  // fill also with rejected hits information
503  SMatrixSym55 fakeCov55;
504  for (int i=0;i<5;i++) for (int j=i;j<5;j++) fakeCov55(i,j) = util::kBogusD;
505  for (unsigned int rejidx = 0; rejidx<rejectedhsidx.size(); ++rejidx) {
506  const unsigned int originalPos = (reverseHits ? hitstatev.size()-rejectedhsidx[rejidx]-1 : rejectedhsidx[rejidx]);
507  auto& mask = hitflagsv[rejectedhsidx[rejidx]];
508  mask.set(recob::TrajectoryPointFlagTraits::HitIgnored,recob::TrajectoryPointFlagTraits::NoPoint);
509  if (mask.isSet(recob::TrajectoryPointFlagTraits::Rejected)==0) mask.set(recob::TrajectoryPointFlagTraits::ExcludedFromFit);
510  //
511  const auto& hitstate = hitstatev[rejectedhsidx[rejidx]];
512  if (dumpLevel_>2) assert(hitstate.wireId().Plane == inHits[originalPos]->WireID().Plane);
513  trkmkr::OptionalPointElement ope;
514  if (optionals.isTrackFitInfosInit()) {
515  ope.setTrackFitHitInfo( recob::TrackFitHitInfo(hitstate.hitMeas(),hitstate.hitMeasErr2(),
516  SVector5(util::kBogusD,util::kBogusD,util::kBogusD,util::kBogusD,util::kBogusD),fakeCov55,hitstate.wireId()) );
517  }
518  tcbk.addPoint(Point_t(util::kBogusD,util::kBogusD,util::kBogusD), Vector_t(util::kBogusD,util::kBogusD,util::kBogusD), inHits[originalPos],
519  recob::TrajectoryPointFlags(originalPos,mask), 0, ope);
520  }
521 
522  if (dumpLevel_>1) std::cout << "outHits.size()=" << outHits.size() << " inHits.size()=" << inHits.size() << std::endl;
523  if (dumpLevel_>2) assert(outHits.size()==inHits.size());
524 
525  bool propok = true;
526  KFTrackState resultF = propagator->rotateToPlane(propok, fwdUpdTkState[sortedtksidx.front()].trackState(),
527  Plane(fwdUpdTkState[sortedtksidx.front()].position(),fwdUpdTkState[sortedtksidx.front()].momentum()));
528  KFTrackState resultB = propagator->rotateToPlane(propok, fwdUpdTkState[sortedtksidx.back()].trackState(),
529  Plane(fwdUpdTkState[sortedtksidx.back()].position(),fwdUpdTkState[sortedtksidx.back()].momentum()));
530 
531  outTrack = tcbk.finalizeTrack(SMatrixSym55(resultF.covariance()),SMatrixSym55(resultB.covariance()));
532 
533  //if there are points with zero momentum return false
534  size_t point = 0;
535  while (outTrack.HasValidPoint(point)) {
536  if (outTrack.MomentumAtPoint( outTrack.NextValidPoint(point++) ) <= 1.0e-9) return false;
537  }
538 
539  if (dumpLevel_>0) {
540  std::cout << "outTrack vertex=" << outTrack.Start()
541  << "\ndir=" << outTrack.StartDirection()
542  << "\ncov=\n" << outTrack.StartCovariance()
543  << "\nlength=" << outTrack.Length() //<< " inLenght=" << track.Length()
544  << std::endl;
545  }
546 
547  return true;
548 }
recob::TrajectoryPointFlagTraits::NoPoint
static constexpr Flag_t NoPoint
The trajectory point is not defined.
Definition: TrajectoryPointFlags.h:122
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:156
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:114
recob::Track::MomentumAtPoint
double MomentumAtPoint(unsigned int p) const
Definition: Track.h:144
recob::Track::StartDirection
Vector_t StartDirection() const
Access to track direction at different points.
Definition: Track.h:134
recob::Track::Length
double Length(size_t p=0) const
Access to various track properties.
Definition: Track.h:170
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:126
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:32
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:98
trkf::TrackKalmanFitter::propagator
const TrackStatePropagator * propagator
Definition: TrackKalmanFitter.h:168
recob::TrajectoryPointFlagTraits::HitIgnored
static constexpr Flag_t HitIgnored
Hit was not included for the computation of the trajectory.
Definition: TrajectoryPointFlags.h:119
recob::TrajectoryPointFlagTraits::Rejected
static constexpr Flag_t Rejected
The hit is extraneous to this track.
Definition: TrajectoryPointFlags.h:175
recob::TrajectoryPointFlagTraits::ExcludedFromFit
static constexpr Flag_t ExcludedFromFit
Definition: TrajectoryPointFlags.h:172
SMatrixSym55
recob::tracking::SMatrixSym55 SMatrixSym55
Definition: PMAlgTracking.cxx:19
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:110
trkmkr::TrackCreationBookKeeper
Helper class to aid the creation of a recob::Track, keeping data vectors in sync. ...
Definition: TrackCreationBookKeeper.h:31
trkmkr::OptionalOutputs::isTrackFitInfosInit
bool isTrackFitInfosInit()
check initialization of the output vector of TrackFitHitInfos
Definition: TrackMaker.h:110
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:35
e
Float_t e
Definition: plot.C:34
trkf::Plane
recob::tracking::Plane Plane
Definition: TrackState.h:17
recob::TrajectoryPointFlags
Set of flags pertaining a point of the track.
Definition: TrajectoryPointFlags.h:318
bool trkf::TrackKalmanFitter::fitTrack ( 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 9 of file TrackKalmanFitter.cxx.

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

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

11  {
12  auto position = traj.Vertex();
13  auto direction = traj.VertexDirection();
14 
15  if (flipDirection) {
16  position = traj.End();
17  direction = -traj.EndDirection();
18  }
19 
20  auto trackStateCov = (flipDirection ? covEnd : covVtx );
21 
22  return fitTrack(position, direction, trackStateCov, hits, traj.Flags(), tkID, pval, pdgid, outTrack, outHits, optionals);
23 }
recob::TrackTrajectory::Flags
Flags_t const & Flags() const
Returns all flags.
Definition: TrackTrajectory.h:171
trkf::TrackKalmanFitter::fitTrack
bool fitTrack(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:9
recob::TrackTrajectory::VertexDirection
Vector_t VertexDirection() const
Returns the direction of the trajectory at the first point.
Definition: TrackTrajectory.h:313
recob::TrackTrajectory::Vertex
Point_t const & Vertex() const
Returns the position of the first valid point of the trajectory [cm].
Definition: TrackTrajectory.h:249
recob::TrackTrajectory::EndDirection
Vector_t EndDirection() const
Returns the direction of the trajectory at the last point.
Definition: TrackTrajectory.h:321
recob::TrackTrajectory::End
Point_t const & End() const
Returns the position of the last valid point of the trajectory [cm].
Definition: TrackTrajectory.h:257
bool trkf::TrackKalmanFitter::fitTrack ( 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 25 of file TrackKalmanFitter.cxx.

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

28  {
29  if (hits.size()<4) {
30  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__;
31  return false;
32  }
33 
34  // setup the KFTrackState we'll use throughout the fit
35  KFTrackState trackState = setupInitialTrackState(position, direction, trackStateCov, pval, pdgid);
36 
37  // setup vector of HitStates and flags, with either same or inverse order as input hit vector
38  // this is what we'll loop over during the fit
39  bool reverseHits = false;
40  std::vector<HitState> hitstatev;
41  std::vector<recob::TrajectoryPointFlags::Mask_t> hitflagsv;
42  bool inputok = setupInputStates(hits, flags, trackState, reverseHits, hitstatev, hitflagsv);
43  if (!inputok) return false;
44 
45  // track and index vectors we use to store the fit results
46  std::vector<KFTrackState> fwdPrdTkState;
47  std::vector<KFTrackState> fwdUpdTkState;
48  std::vector<unsigned int> hitstateidx;
49  std::vector<unsigned int> rejectedhsidx;
50  std::vector<unsigned int> sortedtksidx;
51 
52  // do the actual fit
53  bool fitok = doFitWork(trackState, hitstatev, hitflagsv, fwdPrdTkState, fwdUpdTkState, hitstateidx, rejectedhsidx, sortedtksidx);
54  if (!fitok && (skipNegProp_ || cleanZigzag_) && tryNoSkipWhenFails_) {
55  mf::LogWarning("TrackKalmanFitter") << "Trying to recover with skipNegProp = false and cleanZigzag = false\n";
56  fitok = doFitWork(trackState, hitstatev, hitflagsv, fwdPrdTkState, fwdUpdTkState, hitstateidx, rejectedhsidx, sortedtksidx, false);
57  }
58  if (!fitok) {
59  mf::LogWarning("TrackKalmanFitter") << "Fit failed for track with ID=" << tkID << "\n";
60  return false;
61  }
62 
63  // fill the track, the output hit vector, and the optional output products
64  bool fillok = fillResult(hits, tkID, pdgid, reverseHits, hitstatev, hitflagsv, fwdPrdTkState, fwdUpdTkState, hitstateidx, rejectedhsidx, sortedtksidx, outTrack, outHits, optionals);
65  return fillok;
66 }
trkf::TrackKalmanFitter::setupInputStates
bool setupInputStates(const std::vector< art::Ptr< recob::Hit > > &hits, const std::vector< recob::TrajectoryPointFlags > &flags, const KFTrackState &trackState, bool &reverseHits, 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:83
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:68
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:213
trkf::TrackKalmanFitter::doFitWork
bool doFitWork(KFTrackState &trackState, 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:135
trkf::TrackKalmanFitter::fillResult
bool fillResult(const std::vector< art::Ptr< recob::Hit > > &inHits, const int tkID, const int pdgid, const bool reverseHits, 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:477
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 68 of file TrackKalmanFitter.cxx.

Referenced by fitTrack().

69  {
70  //start from large enough covariance matrix so that the fit is not biased
71  if (trackStateCov==SMatrixSym55()) {
72  trackStateCov(0, 0) = 1000.;
73  trackStateCov(1, 1) = 1000.;
74  trackStateCov(2, 2) = 0.25;
75  trackStateCov(3, 3) = 0.25;
76  trackStateCov(4, 4) = 10.;
77  } else trackStateCov*=100.;
78  // 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)
79  SVector5 trackStatePar(0.,0.,0.,0.,1./pval);
80  return KFTrackState(trackStatePar, trackStateCov, Plane(position,direction), true, pdgid);//along direction by definition
81 }
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 ( const std::vector< art::Ptr< recob::Hit > > &  hits,
const std::vector< recob::TrajectoryPointFlags > &  flags,
const KFTrackState &  trackState,
bool &  reverseHits,
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 83 of file TrackKalmanFitter.cxx.

References detinfo::DetectorProperties::ConvertTicksToX(), recob::TrajectoryPointFlags::DefaultFlagsMask(), detprop, trkf::TrackStatePropagator::distanceToPlane(), dumpLevel_, evd::details::end(), recob::TrajectoryPointFlagTraits::ExcludedFromFit, geom, detinfo::DetectorProperties::GetXTicksCoefficient(), hitErr2ScaleFact_, hits(), recob::tracking::makePlane(), recob::TrajectoryPointFlagTraits::Merged, propagator, rejectHighMultHits_, rejectHitsNegativeGOF_, recob::TrajectoryPointFlagTraits::Suspicious, trkf::KFTrackState::trackState(), useRMS_, geo::GeometryCore::WireIDToWireGeo(), and x.

Referenced by fitTrack().

85  {
86  // figure out hit sorting based on minimum distance to first or last hit
87  // (to be removed once there is a clear convention for hit sorting)
88  geo::WireGeo const& wgeomF = geom->WireIDToWireGeo(hits.front()->WireID());
89  geo::WireGeo const& wgeomB = geom->WireIDToWireGeo(hits.back()->WireID());
90  Plane pF = recob::tracking::makePlane(wgeomF);
91  Plane pB = recob::tracking::makePlane(wgeomB);
92  bool success = true;
93  double distF = propagator->distanceToPlane(success, trackState.trackState(), pF);
94  double distB = propagator->distanceToPlane(success, trackState.trackState(), pB);
95  if (dumpLevel_>1) std::cout << "distF=" << distF << " distB=" << distB << std::endl;
96  if (!success) {
97  mf::LogWarning("TrackKalmanFitter") << "Fit failure at " << __FILE__ << " " << __LINE__;
98  return false;
99  }
100  reverseHits = distB<distF;
101 
102  if (dumpLevel_>1) std::cout << "flags.size()=" << flags.size() << " hits.size()=" << hits.size() << " reverseHits=" << reverseHits << std::endl;
103 
104  // setup vector of HitStates and flags, with either same or inverse order as input hit vector
105  // this is what we'll loop over during the fit
106  const int fsize = flags.size();
107  const int beg = (reverseHits ? hits.size()-1 : 0);
108  const int end = (reverseHits ? -1 : hits.size());
109  for (int ihit = beg; ihit!=end; (reverseHits ? ihit-- : ihit++)) {
110  const auto& hit = hits[ihit];
111  double t = hit->PeakTime();
112  double terr = (useRMS_ ? hit->RMS() : hit->SigmaPeakTime() );
113  double x = detprop->ConvertTicksToX(t, hit->WireID().Plane, hit->WireID().TPC, hit->WireID().Cryostat);
114  double xerr = terr * detprop->GetXTicksCoefficient();
115  hitstatev.emplace_back( x,hitErr2ScaleFact_*xerr*xerr,hit->WireID(),geom->WireIDToWireGeo(hit->WireID()) );
116  //
117  if (fsize>0 && ihit<fsize) hitflagsv.push_back( flags[ihit].mask() );
118  else hitflagsv.push_back(recob::TrajectoryPointFlags::DefaultFlagsMask());
119  //
120  if (dumpLevel_>2) std::cout << "pushed flag mask=" << hitflagsv.back() << std::endl;
121  //
122  if (rejectHighMultHits_ && hit->Multiplicity()>1) {
123  hitflagsv.back().set(recob::TrajectoryPointFlagTraits::Merged);
124  hitflagsv.back().set(recob::TrajectoryPointFlagTraits::ExcludedFromFit);
125  }
126  if (rejectHitsNegativeGOF_ && hit->GoodnessOfFit()<0) {
127  hitflagsv.back().set(recob::TrajectoryPointFlagTraits::Suspicious);
128  hitflagsv.back().set(recob::TrajectoryPointFlagTraits::ExcludedFromFit);
129  }
130  }
131  if (dumpLevel_>2) assert(hits.size()==hitstatev.size());
132  return true;
133 }
x
Float_t x
Definition: compare.C:6
geo::WireGeo
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:61
recob::TrajectoryPointFlagTraits::Merged
static constexpr Flag_t Merged
The hit might have contribution from particles other than this.
Definition: TrajectoryPointFlags.h:135
recob::TrajectoryPointFlagTraits::Suspicious
static constexpr Flag_t Suspicious
The point reconstruction is somehow questionable.
Definition: TrajectoryPointFlags.h:132
recob::TrajectoryPointFlags::DefaultFlagsMask
static constexpr Mask_t DefaultFlagsMask()
Flags used in default construction.
Definition: TrajectoryPointFlags.h:792
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:178
trkf::TrackKalmanFitter::detprop
const detinfo::DetectorProperties * detprop
Definition: TrackKalmanFitter.h:167
detinfo::DetectorProperties::GetXTicksCoefficient
virtual double GetXTicksCoefficient(int t, int c) const =0
recob::tracking::makePlane
Plane makePlane(recob::tracking::Point_t const &pos, recob::tracking::Vector_t const &dir)
helper function to construct a recob::tracking::Plane from a Point_t and a Vector_t; the point is on ...
Definition: TrackingPlaneHelper.h:20
trkf::TrackKalmanFitter::propagator
const TrackStatePropagator * propagator
Definition: TrackKalmanFitter.h:168
trkf::TrackKalmanFitter::geom
art::ServiceHandle< geo::Geometry > geom
Definition: TrackKalmanFitter.h:166
hit
Detector simulation of raw signals on wires.
Definition: DumpHits_module.cc:25
recob::TrajectoryPointFlagTraits::ExcludedFromFit
static constexpr Flag_t ExcludedFromFit
Definition: TrajectoryPointFlags.h:172
detinfo::DetectorProperties::ConvertTicksToX
virtual double ConvertTicksToX(double ticks, int p, int t, int c) const =0
mf::LogWarning
MaybeLogger_< ELseverityLevel::ELsev_warning, false > LogWarning
Definition: MessageLogger.h:213
evd::details::end
std::vector< evd::details::RawDigitInfo_t >::const_iterator end(RawDigitCacheDataClass const &cache)
Definition: RawDataDrawer.cxx:317
trkf::Plane
recob::tracking::Plane Plane
Definition: TrackState.h:17
geo::GeometryCore::WireIDToWireGeo
WireGeo const & WireIDToWireGeo(geo::WireID const &wireid) const
Returns the specified wire.
Definition: GeometryCore.h:3372
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,
bool  applySkipClean = true 
) const
private

Sort the output states.

Definition at line 394 of file TrackKalmanFitter.cxx.

References cleanZigzag_, dir, evd::details::end(), geom, geo::GeometryCore::MaxPlanes(), skipNegProp_, and sortOutputHitsMinLength_.

Referenced by doFitWork().

396  {
397  //
398  if (sortOutputHitsMinLength_) {
399  //try to sort fixing wires order on planes and picking the closest next plane
400  const unsigned int nplanes = geom->MaxPlanes();
401  std::vector<std::vector<unsigned int> > tracksInPlanes(nplanes);
402  for (unsigned int p = 0; p<hitstateidx.size(); ++p) {
403  const auto& hitstate = hitstatev[hitstateidx[p]];
404  tracksInPlanes[hitstate.wireId().Plane].push_back(p);
405  }
406  //this assumes that the first hit/state is a good one, may want to check if that's the case
407  std::vector<unsigned int> iterTracksInPlanes;
408  for (auto it : tracksInPlanes) iterTracksInPlanes.push_back(0);
409  auto pos = fwdUpdTkState.front().position();
410  auto dir = fwdUpdTkState.front().momentum();
411  for (unsigned int p = 0; p<fwdUpdTkState.size(); ++p) {
412  int min_plane = -1;
413  double min_dotp = DBL_MAX;
414  for (unsigned int iplane = 0; iplane<iterTracksInPlanes.size(); ++iplane) {
415  for (unsigned int& itk = iterTracksInPlanes[iplane]; itk<tracksInPlanes[iplane].size(); ++itk) {
416  auto& trackstate = fwdUpdTkState[tracksInPlanes[iplane][iterTracksInPlanes[iplane]]];
417  auto& tmppos = trackstate.position();
418  const double dotp = dir.Dot(tmppos-pos);
419  if (dotp<min_dotp) {
420  min_plane = iplane;
421  min_dotp = dotp;
422  }
423  break;
424  }
425  }
426  if (min_plane<0) continue;
427  const unsigned int ihit = tracksInPlanes[min_plane][iterTracksInPlanes[min_plane]];
428  if (applySkipClean && skipNegProp_ && min_dotp<0.) {
429  rejectedhsidx.push_back(hitstateidx[ihit]);
430  iterTracksInPlanes[min_plane]++;
431  continue;
432  }
433  auto& trackstate = fwdUpdTkState[ihit];
434  pos = trackstate.position();
435  dir = trackstate.momentum();
436  //
437  sortedtksidx.push_back(ihit);
438  iterTracksInPlanes[min_plane]++;
439  }
440  } else {
441  for (unsigned int p = 0; p<fwdUpdTkState.size(); ++p) {
442  sortedtksidx.push_back(p);
443  }
444  }
445  //
446  if (applySkipClean && cleanZigzag_) {
447  std::vector<unsigned int> itoerase;
448  bool clean = false;
449  while (!clean) {
450  bool broken = false;
451  auto pos0 = fwdUpdTkState[sortedtksidx[0]].position();
452  unsigned int i=1;
453  unsigned int end=sortedtksidx.size()-1;
454  for (;i<end;++i) {
455  auto dir0 = fwdUpdTkState[sortedtksidx[i]].position()-pos0;
456  auto dir2 = fwdUpdTkState[sortedtksidx[i+1]].position()-fwdUpdTkState[sortedtksidx[i]].position();
457  dir0/=dir0.R();
458  dir2/=dir2.R();
459  if (dir2.Dot(dir0)<0.) {
460  broken = true;
461  end--;
462  break;
463  } else pos0 = fwdUpdTkState[sortedtksidx[i]].position();
464  }
465  if (!broken) {
466  clean = true;
467  } else {
468  rejectedhsidx.push_back(hitstateidx[sortedtksidx[i]]);
469  sortedtksidx.erase(sortedtksidx.begin()+i);
470  }
471  }
472  }
473  //
474 }
geo::GeometryCore::MaxPlanes
unsigned int MaxPlanes() const
Returns the largest number of planes among all TPCs in this detector.
Definition: GeometryCore.cxx:837
trkf::TrackKalmanFitter::geom
art::ServiceHandle< geo::Geometry > geom
Definition: TrackKalmanFitter.h:166
dir
TDirectory * dir
Definition: macro.C:5
evd::details::end
std::vector< evd::details::RawDigitInfo_t >::const_iterator end(RawDigitCacheDataClass const &cache)
Definition: RawDataDrawer.cxx:317

Member Data Documentation

bool trkf::TrackKalmanFitter::cleanZigzag_
private

Definition at line 173 of file TrackKalmanFitter.h.

Referenced by fitTrack(), and sortOutput().

const detinfo::DetectorProperties* trkf::TrackKalmanFitter::detprop
private

Definition at line 167 of file TrackKalmanFitter.h.

Referenced by setupInputStates().

int trkf::TrackKalmanFitter::dumpLevel_
private

Definition at line 178 of file TrackKalmanFitter.h.

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

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

Definition at line 166 of file TrackKalmanFitter.h.

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

float trkf::TrackKalmanFitter::hitErr2ScaleFact_
private

Definition at line 176 of file TrackKalmanFitter.h.

Referenced by setupInputStates().

const TrackStatePropagator* trkf::TrackKalmanFitter::propagator
private

Definition at line 168 of file TrackKalmanFitter.h.

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

bool trkf::TrackKalmanFitter::rejectHighMultHits_
private

Definition at line 174 of file TrackKalmanFitter.h.

Referenced by setupInputStates().

bool trkf::TrackKalmanFitter::rejectHitsNegativeGOF_
private

Definition at line 175 of file TrackKalmanFitter.h.

Referenced by setupInputStates().

bool trkf::TrackKalmanFitter::skipNegProp_
private

Definition at line 172 of file TrackKalmanFitter.h.

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

bool trkf::TrackKalmanFitter::sortHitsByPlane_
private

Definition at line 170 of file TrackKalmanFitter.h.

Referenced by doFitWork().

bool trkf::TrackKalmanFitter::sortOutputHitsMinLength_
private

Definition at line 171 of file TrackKalmanFitter.h.

Referenced by sortOutput().

bool trkf::TrackKalmanFitter::tryNoSkipWhenFails_
private

Definition at line 177 of file TrackKalmanFitter.h.

Referenced by fitTrack().

bool trkf::TrackKalmanFitter::useRMS_
private

Definition at line 169 of file TrackKalmanFitter.h.

Referenced by setupInputStates().

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