Public Member Functions | Private Attributes | List of all members
trkf::Track3DKalmanHitAlg Class Reference

#include <Track3DKalmanHitAlg.h>

Public Member Functions

 Track3DKalmanHitAlg (const fhicl::ParameterSet &pset)
 Constructor. More...
 
std::vector< trkf::KalmanOutputmakeTracks (detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, KalmanInputs &kalman_inputs)
 
void fetchPFParticleSeeds (const art::PtrVector< recob::Seed > &pfseeds, const std::vector< Hits > &pfseedhits, std::vector< recob::Seed > &seeds, std::vector< Hits > &hitsperseed) const
 Fetch Seeds method. More...
 
recob::Seed makeSeed (detinfo::DetectorPropertiesData const &detProp, const Hits &hits) const
 Make seed method. More...
 
void growSeedsIntoTracks (detinfo::DetectorPropertiesData const &detProp, const bool pfseed, const std::vector< recob::Seed > &seeds, const std::vector< Hits > &hitsperseed, Hits &unusedhits, Hits &hits, std::deque< KGTrack > &kalman_tracks)
 Grow Seeds method. More...
 
void growSeedIntoTracks (detinfo::DetectorPropertiesData const &detProp, const bool pfseed, const recob::Seed &seed, const Hits &hpsit, Hits &unusedhits, Hits &hits, std::deque< KGTrack > &kgtracks)
 
void chopHitsOffSeeds (Hits const &hpsit, bool pfseed, Hits &seedhits) const
 Chop hits off of the end of seeds. More...
 
bool testSeedSlope (const double *dir) const
 
std::shared_ptr< SurfacemakeSurface (const recob::Seed &seed, double *dir) const
 method to return a seed to surface. More...
 
bool makeKalmanTracks (detinfo::DetectorPropertiesData const &detProp, const std::shared_ptr< trkf::Surface > psurf, const Surface::TrackDirection trkdir, Hits &seedhits, Hits &hits, std::deque< KGTrack > &kalman_tracks)
 
bool smoothandextendTrack (detinfo::DetectorPropertiesData const &detProp, Propagator const &propagator, KGTrack &trg0, const Hits hits, unsigned int prefplane, std::deque< KGTrack > &kalman_tracks)
 SMooth and extend track. More...
 
bool extendandsmoothLoop (detinfo::DetectorPropertiesData const &detProp, Propagator const &propagator, KGTrack &trg1, unsigned int prefplane, Hits &trackhits) const
 SMooth and extend a track in a loop. More...
 
void filterHitsOnKalmanTrack (const KGTrack &trg, Hits &hits, Hits &seederhits) const
 Filter hits that are on kalman tracks. More...
 
std::unique_ptr< KHitContainerfillHitContainer (detinfo::DetectorPropertiesData const &detProp, const Hits &hits) const
 Fill hit container with either seedhits or filtered hits i.e. recob::Hit. More...
 
bool qualityCutsOnSeedTrack (const KGTrack &trg0) const
 Quality cuts on seed track. More...
 
void fitnupdateMomentum (Propagator const &propagator, KGTrack &trg1, KGTrack &trg2) const
 fit and update method, used twice. More...
 

Private Attributes

bool fDoDedx
 Global dE/dx enable flag. More...
 
bool fSelfSeed
 Self seed flag. More...
 
double fMaxTcut
 Maximum delta ray energy in MeV for restricted dE/dx. More...
 
bool fLineSurface
 Line surface flag. More...
 
size_t fMinSeedHits
 Minimum number of hits per track seed. More...
 
int fMinSeedChopHits
 Potentially chop seeds that exceed this length. More...
 
int fMaxChopHits
 Maximum number of hits to chop from each end of seed. More...
 
double fMaxSeedChiDF
 Maximum seed track chisquare/dof. More...
 
double fMinSeedSlope
 Minimum seed slope (dx/dz). More...
 
double fInitialMomentum
 Initial (or constant) momentum. More...
 
KalmanFilterAlg fKFAlg
 Kalman filter algorithm. More...
 
SeedFinderAlgorithm fSeedFinderAlg
 Seed finder. More...
 
int fNumTrack
 Number of tracks produced. More...
 

Detailed Description

Definition at line 55 of file Track3DKalmanHitAlg.h.

Constructor & Destructor Documentation

trkf::Track3DKalmanHitAlg::Track3DKalmanHitAlg ( const fhicl::ParameterSet pset)
explicit

Constructor.

Definition at line 84 of file Track3DKalmanHitAlg.cxx.

85  : fDoDedx{pset.get<bool>("DoDedx")}
86  , fSelfSeed{pset.get<bool>("SelfSeed")}
87  , fMaxTcut{pset.get<double>("MaxTcut")}
88  , fLineSurface{pset.get<bool>("LineSurface")}
89  , fMinSeedHits{pset.get<size_t>("MinSeedHits")}
90  , fMinSeedChopHits{pset.get<int>("MinSeedChopHits")}
91  , fMaxChopHits{pset.get<int>("MaxChopHits")}
92  , fMaxSeedChiDF{pset.get<double>("MaxSeedChiDF")}
93  , fMinSeedSlope{pset.get<double>("MinSeedSlope")}
94  , fInitialMomentum{pset.get<double>("InitialMomentum")}
95  , fKFAlg(pset.get<fhicl::ParameterSet>("KalmanFilterAlg"))
96  , fSeedFinderAlg(pset.get<fhicl::ParameterSet>("SeedFinderAlg"))
97  , fNumTrack(0)
98 {
99  mf::LogInfo("Track3DKalmanHitAlg") << "Track3DKalmanHitAlg instantiated.";
100 }
size_t fMinSeedHits
Minimum number of hits per track seed.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
bool fDoDedx
Global dE/dx enable flag.
double fMaxSeedChiDF
Maximum seed track chisquare/dof.
int fMinSeedChopHits
Potentially chop seeds that exceed this length.
double fInitialMomentum
Initial (or constant) momentum.
T get(std::string const &key) const
Definition: ParameterSet.h:271
KalmanFilterAlg fKFAlg
Kalman filter algorithm.
int fNumTrack
Number of tracks produced.
bool fLineSurface
Line surface flag.
SeedFinderAlgorithm fSeedFinderAlg
Seed finder.
int fMaxChopHits
Maximum number of hits to chop from each end of seed.
bool fSelfSeed
Self seed flag.
double fMinSeedSlope
Minimum seed slope (dx/dz).
double fMaxTcut
Maximum delta ray energy in MeV for restricted dE/dx.

Member Function Documentation

void trkf::Track3DKalmanHitAlg::chopHitsOffSeeds ( Hits const &  hpsit,
bool  pfseed,
Hits seedhits 
) const

Chop hits off of the end of seeds.

Definition at line 389 of file Track3DKalmanHitAlg.cxx.

390 {
391  const int nchopmax =
392  (pfseed || fSelfSeed) ? 0 : std::max(0, int((hpsit.size() - fMinSeedChopHits) / 2));
393  //SS: FIXME
394 
395  const int nchop = std::min(nchopmax, fMaxChopHits);
396  Hits::const_iterator itb = hpsit.begin();
397  Hits::const_iterator ite = hpsit.end();
398  itb += nchop;
399  ite -= nchop;
400  seedhits.reserve(hpsit.size());
401  for (Hits::const_iterator it = itb; it != ite; ++it)
402  seedhits.push_back(*it);
403 }
int fMinSeedChopHits
Potentially chop seeds that exceed this length.
typename data_t::const_iterator const_iterator
Definition: PtrVector.h:55
static int max(int a, int b)
T min(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:55
int fMaxChopHits
Maximum number of hits to chop from each end of seed.
bool fSelfSeed
Self seed flag.
bool trkf::Track3DKalmanHitAlg::extendandsmoothLoop ( detinfo::DetectorPropertiesData const &  detProp,
Propagator const &  propagator,
KGTrack trg1,
unsigned int  prefplane,
Hits trackhits 
) const

SMooth and extend a track in a loop.

Definition at line 459 of file Track3DKalmanHitAlg.cxx.

464 {
465  bool ok = true;
466  const int niter = 6;
467  int nfail = 0; // Number of consecutive extend fails.
468  for (int ix = 0; ok && ix < niter && nfail < 2; ++ix) {
469 
470  // Fill a collection of hits from the last good track
471  // (initially the seed track).
472  Hits goodhits;
473  std::vector<unsigned int> hittpindex;
474  trg1.fillHits(goodhits, hittpindex);
475 
476  // Filter hits already on the track out of the available hits.
477  filterHits(trackhits, goodhits);
478 
479  // Fill hit container using filtered hits.
480  std::unique_ptr<KHitContainer> ptrackcont = fillHitContainer(detProp, trackhits);
481 
482  // Extend the track. It is not an error for the
483  // extend operation to fail, meaning that no new hits
484  // were added.
485 
486  if (fKFAlg.extendTrack(trg1, propagator, *ptrackcont))
487  nfail = 0;
488  else
489  ++nfail;
490 
491  // Smooth the extended track, and make a new
492  // unidirectionally fit track in the opposite
493  // direction.
494 
495  KGTrack trg2(prefplane);
496  ok = fKFAlg.smoothTrack(trg1, &trg2, propagator);
497  if (ok) {
498  // Skip momentum estimate for constant-momentum tracks.
499  if (fDoDedx) { fitnupdateMomentum(propagator, trg1, trg2); }
500  trg1 = trg2;
501  }
502  }
503  return ok;
504 }
bool fDoDedx
Global dE/dx enable flag.
bool extendTrack(KGTrack &trg, const Propagator &prop, KHitContainer &hits) const
Add hits to existing track.
std::unique_ptr< KHitContainer > fillHitContainer(detinfo::DetectorPropertiesData const &detProp, const Hits &hits) const
Fill hit container with either seedhits or filtered hits i.e. recob::Hit.
bool smoothTrack(KGTrack &trg, KGTrack *trg1, const Propagator &prop) const
Smooth track.
KalmanFilterAlg fKFAlg
Kalman filter algorithm.
art::PtrVector< recob::Hit > Hits
void fitnupdateMomentum(Propagator const &propagator, KGTrack &trg1, KGTrack &trg2) const
fit and update method, used twice.
void trkf::Track3DKalmanHitAlg::fetchPFParticleSeeds ( const art::PtrVector< recob::Seed > &  pfseeds,
const std::vector< Hits > &  pfseedhits,
std::vector< recob::Seed > &  seeds,
std::vector< Hits > &  hitsperseed 
) const

Fetch Seeds method.

Definition at line 168 of file Track3DKalmanHitAlg.cxx.

172 {
173  for (const auto& pseed : pfseeds) {
174  seeds.push_back(*pseed);
175  }
176  hitsperseed.insert(hitsperseed.end(), pfseedhits.begin(), pfseedhits.end());
177 }
std::unique_ptr< trkf::KHitContainer > trkf::Track3DKalmanHitAlg::fillHitContainer ( detinfo::DetectorPropertiesData const &  detProp,
const Hits hits 
) const

Fill hit container with either seedhits or filtered hits i.e. recob::Hit.

Definition at line 353 of file Track3DKalmanHitAlg.cxx.

355 {
356  std::unique_ptr<KHitContainer> hitcont(fLineSurface ?
357  static_cast<KHitContainer*>(new KHitContainerWireLine) :
358  static_cast<KHitContainer*>(new KHitContainerWireX));
359  hitcont->fill(detProp, hits, -1);
360  return hitcont;
361 }
bool fLineSurface
Line surface flag.
void trkf::Track3DKalmanHitAlg::filterHitsOnKalmanTrack ( const KGTrack trg,
Hits hits,
Hits seederhits 
) const

Filter hits that are on kalman tracks.

Definition at line 339 of file Track3DKalmanHitAlg.cxx.

342 {
343  Hits track_used_hits;
344  std::vector<unsigned int> hittpindex;
345  trg.fillHits(track_used_hits, hittpindex);
346  filterHits(hits, track_used_hits);
347  filterHits(seederhits, track_used_hits);
348 }
art::PtrVector< recob::Hit > Hits
void trkf::Track3DKalmanHitAlg::fitnupdateMomentum ( Propagator const &  propagator,
KGTrack trg1,
KGTrack trg2 
) const

fit and update method, used twice.

Definition at line 508 of file Track3DKalmanHitAlg.cxx.

511 {
512  KETrack tremom;
513  if (fKFAlg.fitMomentum(trg1, propagator, tremom)) {
514  fKFAlg.updateMomentum(tremom, propagator, trg2);
515  }
516 }
KalmanFilterAlg fKFAlg
Kalman filter algorithm.
bool fitMomentum(const KGTrack &trg, const Propagator &prop, KETrack &tremom) const
Estimate track momentum using either range or multiple scattering.
bool updateMomentum(const KETrack &tremom, const Propagator &prop, KGTrack &trg) const
Set track momentum at each track surface.
void trkf::Track3DKalmanHitAlg::growSeedIntoTracks ( detinfo::DetectorPropertiesData const &  detProp,
const bool  pfseed,
const recob::Seed seed,
const Hits hpsit,
Hits unusedhits,
Hits hits,
std::deque< KGTrack > &  kgtracks 
)

Definition at line 203 of file Track3DKalmanHitAlg.cxx.

210 {
211  Hits trimmedhits;
212  // Chop a couple of hits off each end of the seed.
213  chopHitsOffSeeds(hpsit, pfseed, trimmedhits);
214 
215  // Filter hits used by (chopped) seed from hits available to make future seeds.
216  // No matter what, we will never use these hits for another seed.
217  // This eliminates the possibility of an infinite loop.
218 
219  size_t initial_unusedhits = unusedhits.size();
220  filterHits(unusedhits, trimmedhits);
221 
222  // Require that this seed be fully disjoint from existing tracks.
223  //SS: replace this test with a method with appropriate name
224  if (!(trimmedhits.size() + unusedhits.size() == initial_unusedhits)) return;
225 
226  // Convert seed into initial KTracks on surface located at seed point,
227  // and normal to seed direction.
228  double dir[3];
229  std::shared_ptr<Surface> psurf = makeSurface(seed, dir);
230 
231  // Cut on the seed slope dx/ds.
232  //SS: replace test name with a reasonable name
233  if (!testSeedSlope(dir)) return;
234 
235  // Make one or two initial KTracks for forward and backward directions.
236  // The build_both flag specifies whether we should attempt to make
237  // tracks from all both FORWARD and BACKWARD initial tracks,
238  // or alternatively, whether we
239  // should declare victory and quit after getting a successful
240  // track from one initial track.
241  const bool build_both = fDoDedx;
242  const int ninit = 2;
243 
244  auto ntracks = kgtracks.size(); // Remember original track count.
245  bool ok = makeKalmanTracks(detProp, psurf, Surface::FORWARD, trimmedhits, hits, kgtracks);
246  if ((!ok || build_both) && ninit == 2) {
247  makeKalmanTracks(detProp, psurf, Surface::BACKWARD, trimmedhits, hits, kgtracks);
248  }
249 
250  // Loop over newly added tracks and remove hits contained on
251  // these tracks from hits available for making additional
252  // tracks or track seeds.
253  for (unsigned int itrk = ntracks; itrk < kgtracks.size(); ++itrk) {
254  const KGTrack& trg = kgtracks[itrk];
255  filterHitsOnKalmanTrack(trg, hits, unusedhits);
256  }
257 }
int ntracks
Definition: tracks.py:246
void filterHitsOnKalmanTrack(const KGTrack &trg, Hits &hits, Hits &seederhits) const
Filter hits that are on kalman tracks.
bool fDoDedx
Global dE/dx enable flag.
string dir
void chopHitsOffSeeds(Hits const &hpsit, bool pfseed, Hits &seedhits) const
Chop hits off of the end of seeds.
bool makeKalmanTracks(detinfo::DetectorPropertiesData const &detProp, const std::shared_ptr< trkf::Surface > psurf, const Surface::TrackDirection trkdir, Hits &seedhits, Hits &hits, std::deque< KGTrack > &kalman_tracks)
art::PtrVector< recob::Hit > Hits
bool testSeedSlope(const double *dir) const
std::shared_ptr< Surface > makeSurface(const recob::Seed &seed, double *dir) const
method to return a seed to surface.
void trkf::Track3DKalmanHitAlg::growSeedsIntoTracks ( detinfo::DetectorPropertiesData const &  detProp,
const bool  pfseed,
const std::vector< recob::Seed > &  seeds,
const std::vector< Hits > &  hitsperseed,
Hits unusedhits,
Hits hits,
std::deque< KGTrack > &  kalman_tracks 
)

Grow Seeds method.

Definition at line 183 of file Track3DKalmanHitAlg.cxx.

190 {
191  // check for size of both containers
192  if (seeds.size() != hitsperseed.size()) {
193  throw cet::exception("Track3DKalmanHitAlg")
194  << "Different size containers for Seeds and Hits/Seed.\n";
195  }
196  for (size_t i = 0; i < seeds.size(); ++i) {
197  growSeedIntoTracks(detProp, pfseed, seeds[i], hitsperseed[i], unusedhits, hits, kgtracks);
198  }
199 }
void growSeedIntoTracks(detinfo::DetectorPropertiesData const &detProp, const bool pfseed, const recob::Seed &seed, const Hits &hpsit, Hits &unusedhits, Hits &hits, std::deque< KGTrack > &kgtracks)
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
bool trkf::Track3DKalmanHitAlg::makeKalmanTracks ( detinfo::DetectorPropertiesData const &  detProp,
const std::shared_ptr< trkf::Surface psurf,
const Surface::TrackDirection  trkdir,
Hits seedhits,
Hits hits,
std::deque< KGTrack > &  kalman_tracks 
)

Definition at line 282 of file Track3DKalmanHitAlg.cxx.

288 {
289  const int pdg = 13; //SS: FIXME another constant?
290  // SS: FIXME
291  // It is a contant value inside a loop so I took it out
292  // revisit the linear algebra stuff used here (use of ublas)
293  // make a lambda here ... const TrackVector
294  //
295  TrackVector vec(5);
296  vec(0) = 0.;
297  vec(1) = 0.;
298  vec(2) = 0.;
299  vec(3) = 0.;
300  vec(4) = (fInitialMomentum != 0. ? 1. / fInitialMomentum : 2.);
301 
302  KTrack initial_track(psurf, vec, trkdir, pdg);
303 
304  // Fill hit container with current seed hits.
305  // KHitContainer may not be cheaply movabale thats why unique_ptr
306  std::unique_ptr<KHitContainer> pseedcont = fillHitContainer(detProp, seedhits);
307 
308  // Set the preferred plane to be the one with the most hits.
309  unsigned int prefplane = pseedcont->getPreferredPlane();
310  fKFAlg.setPlane(prefplane);
311  if (mf::isDebugEnabled())
312  mf::LogDebug("Track3DKalmanHit") << "Preferred plane = " << prefplane << "\n";
313 
314  PropAny const propagator{detProp, fMaxTcut, fDoDedx};
315 
316  // Build and smooth seed track.
317  KGTrack trg0(prefplane);
318  bool ok =
319  fKFAlg.buildTrack(initial_track, trg0, propagator, Propagator::FORWARD, *pseedcont, fSelfSeed);
320  if (ok) ok = smoothandextendTrack(detProp, propagator, trg0, hits, prefplane, kgtracks);
321 
322  if (mf::isDebugEnabled())
323  mf::LogDebug("Track3DKalmanHit")
324  << (ok ? "Find track succeeded." : "Find track failed.") << "\n";
325  return ok;
326 }
bool fDoDedx
Global dE/dx enable flag.
double fInitialMomentum
Initial (or constant) momentum.
bool buildTrack(const KTrack &trk, KGTrack &trg, const Propagator &prop, const Propagator::PropDirection dir, KHitContainer &hits, bool linear) const
Make a new track.
std::unique_ptr< KHitContainer > fillHitContainer(detinfo::DetectorPropertiesData const &detProp, const Hits &hits) const
Fill hit container with either seedhits or filtered hits i.e. recob::Hit.
bool isDebugEnabled()
KVector< 5 >::type TrackVector
Track state vector, dimension 5.
KalmanFilterAlg fKFAlg
Kalman filter algorithm.
bool smoothandextendTrack(detinfo::DetectorPropertiesData const &detProp, Propagator const &propagator, KGTrack &trg0, const Hits hits, unsigned int prefplane, std::deque< KGTrack > &kalman_tracks)
SMooth and extend track.
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
void setPlane(int plane)
Set preferred view plane.
bool fSelfSeed
Self seed flag.
double fMaxTcut
Maximum delta ray energy in MeV for restricted dE/dx.
recob::Seed trkf::Track3DKalmanHitAlg::makeSeed ( detinfo::DetectorPropertiesData const &  detProp,
const Hits hits 
) const

Make seed method.

Definition at line 521 of file Track3DKalmanHitAlg.cxx.

523 {
525 
526  // Do a linear 3D least squares for of y and z vs. x.
527  // y = y0 + ay*(x-x0)
528  // z = z0 + az*(x-x0)
529  // Here x0 is the global average x coordinate of all hits in all planes.
530  // Parameters y0, z0, ay, and az are determined by a simultaneous least squares
531  // fit in all planes.
532 
533  // First, determine x0 by looping over every hit.
534 
535  double x0 = 0.;
536  int n = 0;
537  for (auto const& phit : hits) {
538  const recob::Hit& hit = *phit;
539  geo::WireID wire_id = hit.WireID();
540  double time = hit.PeakTime();
541  double x = detProp.ConvertTicksToX(time, wire_id);
542  x0 += x;
543  ++n;
544  }
545 
546  // If there are no hits, return invalid seed.
547 
548  if (n == 0) return recob::Seed();
549 
550  // Find the average x.
551 
552  x0 /= n;
553 
554  // Now do the least squares fit proper.
555 
556  KSymMatrix<4>::type sm(4);
557  KVector<4>::type sv(4);
558  sm.clear();
559  sv.clear();
560 
561  // Loop over hits (again).
562 
563  for (auto const& phit : hits) {
564  const recob::Hit& hit = *phit;
565 
566  // Extract the angle, w and x coordinates from hit.
567 
568  geo::WireID wire_id = hit.WireID();
569  const geo::WireGeo& wgeom = geom->Wire(wire_id);
570  double xyz[3];
571  wgeom.GetCenter(xyz);
572 
573  // Phi convention is the one documented in SurfYZPlane.h.
574 
575  double phi = TMath::PiOver2() - wgeom.ThetaZ();
576  double sphi = std::sin(phi);
577  double cphi = std::cos(phi);
578  double w = -xyz[1] * sphi + xyz[2] * cphi;
579 
580  double time = hit.PeakTime();
581  double x = detProp.ConvertTicksToX(time, wire_id);
582 
583  // Accumulate data for least squares fit.
584 
585  double dx = x - x0;
586 
587  sm(0, 0) += sphi * sphi;
588  sm(1, 0) -= sphi * cphi;
589  sm(1, 1) += cphi * cphi;
590  sm(2, 0) += sphi * sphi * dx;
591  sm(2, 1) -= sphi * cphi * dx;
592  sm(2, 2) += sphi * sphi * dx * dx;
593  sm(3, 0) -= sphi * cphi * dx;
594  sm(3, 1) += cphi * cphi * dx;
595  sm(3, 2) -= sphi * cphi * dx * dx;
596  sm(3, 3) += cphi * cphi * dx * dx;
597 
598  sv(0) -= sphi * w;
599  sv(1) += cphi * w;
600  sv(2) -= sphi * w * dx;
601  sv(3) += cphi * w * dx;
602  }
603 
604  // Solve.
605 
606  bool ok = syminvert(sm);
607  if (!ok) return recob::Seed();
608  KVector<4>::type par(4);
609  par = prod(sm, sv);
610 
611  double y0 = par(0);
612  double z0 = par(1);
613  double dydx = par(2);
614  double dzdx = par(3);
615 
616  // Make seed.
617 
618  double dsdx = std::hypot(1., std::hypot(dydx, dzdx));
619 
620  double pos[3] = {x0, y0, z0};
621  double dir[3] = {1. / dsdx, dydx / dsdx, dzdx / dsdx};
622  double poserr[3] = {0., 0., 0.};
623  double direrr[3] = {0., 0., 0.};
624 
625  recob::Seed result(pos, dir, poserr, direrr);
626 
627  return result;
628 }
ublas::symmetric_matrix< double, ublas::lower, ublas::row_major, ublas::bounded_array< double, N *(N+1)/2 > > type
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
static QCString result
geo::WireID WireID() const
Definition: Hit.h:233
string dir
double ThetaZ() const
Returns angle of wire with respect to z axis in the Y-Z plane in radians.
Definition: WireGeo.h:250
std::enable_if_t< std::is_arithmetic_v< T >, T > hypot(T x, T y)
Definition: hypot.h:60
ublas::vector< double, ublas::bounded_array< double, N > > type
std::void_t< T > n
bool syminvert(ublas::symmetric_matrix< T, TRI, L, A > &m)
Detector simulation of raw signals on wires.
WireGeo const & Wire(geo::WireID const &wireid) const
Returns the specified wire.
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:218
list x
Definition: train.py:276
2D representation of charge deposited in the TDC/wire plane
Definition: Hit.h:48
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:73
std::shared_ptr< trkf::Surface > trkf::Track3DKalmanHitAlg::makeSurface ( const recob::Seed seed,
double *  dir 
) const

method to return a seed to surface.

Definition at line 263 of file Track3DKalmanHitAlg.cxx.

264 {
265  double xyz[3];
266  double err[3]; // Dummy.
267  seed.GetPoint(xyz, err);
268  seed.GetDirection(dir, err);
269  if (mf::isDebugEnabled()) {
270  mf::LogDebug("Track3DKalmanHit")
271  //<< "Seed found with " << seedhits.size() <<" hits.\n"
272  << "(x,y,z) = " << xyz[0] << ", " << xyz[1] << ", " << xyz[2] << "\n"
273  << "(dx,dy,dz) = " << dir[0] << ", " << dir[1] << ", " << dir[2] << "\n";
274  } // if debug
275 
276  return std::make_shared<SurfXYZPlane>(xyz[0], xyz[1], xyz[2], dir[0], dir[1], dir[2]);
277 }
void GetPoint(double *Pt, double *Err) const
Definition: Seed.cxx:108
string dir
bool isDebugEnabled()
void err(const char *fmt,...)
Definition: message.cpp:226
MaybeLogger_< ELseverityLevel::ELsev_success, false > LogDebug
void GetDirection(double *Dir, double *Err) const
Definition: Seed.cxx:98
std::vector< trkf::KalmanOutput > trkf::Track3DKalmanHitAlg::makeTracks ( detinfo::DetectorClocksData const &  clockData,
detinfo::DetectorPropertiesData const &  detProp,
KalmanInputs kalman_inputs 
)

Definition at line 105 of file Track3DKalmanHitAlg.cxx.

108 {
109  std::vector<KalmanOutput> outputs(inputs.size());
110  // Loop over input vectors.
111  for (size_t i = 0; i < inputs.size(); ++i) {
112  Hits& hits = inputs[i].hits;
113  // The hit collection "hits" (just filled), initially containing all
114  // hits, represents hits available for making tracks. Now we will
115  // fill a second hit collection called "unusedhits", also initially
116  // containing all hits, which will represent hits available for
117  // making track seeds. These collections are not necessarily the
118  // same, since hits that are not suitable for seeds may still be
119  // suitable for tracks.
120  Hits unusedhits = hits;
121  // Start of loop.
122  bool first = true;
123  while (true) {
124  std::vector<Hits> hitsperseed;
125  std::vector<recob::Seed> seeds;
126  // On the first trip through this loop, try to use pfparticle-associated seeds.
127  // Do this, provided the list of pfparticle-associated seeds and associated
128  // hits are not empty.
129  bool pfseed = false;
130  auto const seedsize = inputs[i].seeds.size();
131  if (first && seedsize > 0 && inputs[i].seedhits.size() > 0) {
132  seeds.reserve(seedsize);
133  fetchPFParticleSeeds(inputs[i].seeds, inputs[i].seedhits, seeds, hitsperseed);
134  pfseed = true;
135  }
136  else {
137  // On subsequent trips, or if there were no usable pfparticle-associated seeds,
138  // attempt to generate our own seeds.
139  if (unusedhits.size() > 0) {
140  if (fSelfSeed) {
141  // Self seed - convert all hits into one big seed.
142  seeds.emplace_back(makeSeed(detProp, unusedhits));
143  hitsperseed.emplace_back();
144  hitsperseed.back().insert(
145  hitsperseed.back().end(), unusedhits.begin(), unusedhits.end());
146  }
147  else
148  seeds =
149  fSeedFinderAlg.GetSeedsFromUnSortedHits(clockData, detProp, unusedhits, hitsperseed);
150  }
151  }
152 
153  assert(seeds.size() == hitsperseed.size());
154  first = false;
155 
156  if (empty(seeds)) break;
157 
158  growSeedsIntoTracks(detProp, pfseed, seeds, hitsperseed, unusedhits, hits, outputs[i].tracks);
159  }
160  }
161  return outputs;
162 }
recob::Seed makeSeed(detinfo::DetectorPropertiesData const &detProp, const Hits &hits) const
Make seed method.
void growSeedsIntoTracks(detinfo::DetectorPropertiesData const &detProp, const bool pfseed, const std::vector< recob::Seed > &seeds, const std::vector< Hits > &hitsperseed, Hits &unusedhits, Hits &hits, std::deque< KGTrack > &kalman_tracks)
Grow Seeds method.
void fetchPFParticleSeeds(const art::PtrVector< recob::Seed > &pfseeds, const std::vector< Hits > &pfseedhits, std::vector< recob::Seed > &seeds, std::vector< Hits > &hitsperseed) const
Fetch Seeds method.
std::vector< recob::Seed > GetSeedsFromUnSortedHits(detinfo::DetectorClocksData const &clockData, detinfo::DetectorPropertiesData const &detProp, art::PtrVector< recob::Hit > const &, std::vector< art::PtrVector< recob::Hit >> &, unsigned int StopAfter=0) const
Definition: tracks.py:1
std::vector< TrajPoint > seeds
Definition: DataStructs.cxx:13
art::PtrVector< recob::Hit > Hits
SeedFinderAlgorithm fSeedFinderAlg
Seed finder.
bool fSelfSeed
Self seed flag.
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:97
bool trkf::Track3DKalmanHitAlg::qualityCutsOnSeedTrack ( const KGTrack trg0) const

Quality cuts on seed track.

Definition at line 368 of file Track3DKalmanHitAlg.cxx.

369 {
370  double mom0[3];
371  double mom1[3];
372  trg0.startTrack().getMomentum(mom0);
373  trg0.endTrack().getMomentum(mom1);
374 
375  double dxds0 = mom0[0] / calcMagnitude(mom0);
376  double dxds1 = mom1[0] / calcMagnitude(mom1);
377 
378  return (std::abs(dxds0) > fMinSeedSlope && std::abs(dxds1) > fMinSeedSlope);
379 }
T abs(T value)
double fMinSeedSlope
Minimum seed slope (dx/dz).
bool trkf::Track3DKalmanHitAlg::smoothandextendTrack ( detinfo::DetectorPropertiesData const &  detProp,
Propagator const &  propagator,
KGTrack trg0,
const Hits  hits,
unsigned int  prefplane,
std::deque< KGTrack > &  kalman_tracks 
)

SMooth and extend track.

Definition at line 408 of file Track3DKalmanHitAlg.cxx.

414 {
415  KGTrack trg1(prefplane);
416  bool ok = fKFAlg.smoothTrack(trg0, &trg1, propagator);
417  if (!ok) return ok;
418 
419  // Now we have the seed track in the form of a KGTrack.
420  // Do additional quality cuts.
421 
422  auto const n = trg1.numHits();
423  auto const chisq = n * fMaxSeedChiDF;
424 
425  ok = n >= fMinSeedHits && trg1.startTrack().getChisq() <= chisq &&
426  trg1.endTrack().getChisq() <= chisq && trg0.startTrack().getChisq() <= chisq &&
427  trg0.endTrack().getChisq() <= chisq && qualityCutsOnSeedTrack(trg0);
428 
429  if (!ok) return ok;
430 
431  // Make a copy of the original hit collection of all
432  // available track hits.
433  Hits trackhits = hits;
434 
435  // Do an extend + smooth loop here.
436  // Exit after two consecutive failures to extend (i.e. from each end),
437  // or if the iteration count reaches the maximum.
438  if (ok) ok = extendandsmoothLoop(detProp, propagator, trg1, prefplane, trackhits);
439 
440  // Do a final smooth.
441  if (!ok) return false;
442 
443  ok = fKFAlg.smoothTrack(trg1, 0, propagator);
444  if (!ok) return false;
445 
446  // Skip momentum estimate for constant-momentum tracks.
447 
448  if (fDoDedx) { fitnupdateMomentum(propagator, trg1, trg1); }
449  // Save this track.
450  ++fNumTrack;
451  kalman_tracks.push_back(trg1);
452  return true;
453 }
size_t fMinSeedHits
Minimum number of hits per track seed.
bool fDoDedx
Global dE/dx enable flag.
double fMaxSeedChiDF
Maximum seed track chisquare/dof.
std::void_t< T > n
bool smoothTrack(KGTrack &trg, KGTrack *trg1, const Propagator &prop) const
Smooth track.
KalmanFilterAlg fKFAlg
Kalman filter algorithm.
art::PtrVector< recob::Hit > Hits
void fitnupdateMomentum(Propagator const &propagator, KGTrack &trg1, KGTrack &trg2) const
fit and update method, used twice.
int fNumTrack
Number of tracks produced.
bool qualityCutsOnSeedTrack(const KGTrack &trg0) const
Quality cuts on seed track.
bool extendandsmoothLoop(detinfo::DetectorPropertiesData const &detProp, Propagator const &propagator, KGTrack &trg1, unsigned int prefplane, Hits &trackhits) const
SMooth and extend a track in a loop.
bool trkf::Track3DKalmanHitAlg::testSeedSlope ( const double *  dir) const

Definition at line 331 of file Track3DKalmanHitAlg.cxx.

332 {
333  return std::abs(dir[0]) >= fMinSeedSlope * calcMagnitude(dir);
334 }
string dir
T abs(T value)
double fMinSeedSlope
Minimum seed slope (dx/dz).

Member Data Documentation

bool trkf::Track3DKalmanHitAlg::fDoDedx
private

Global dE/dx enable flag.

Definition at line 111 of file Track3DKalmanHitAlg.h.

double trkf::Track3DKalmanHitAlg::fInitialMomentum
private

Initial (or constant) momentum.

Definition at line 120 of file Track3DKalmanHitAlg.h.

KalmanFilterAlg trkf::Track3DKalmanHitAlg::fKFAlg
private

Kalman filter algorithm.

Definition at line 124 of file Track3DKalmanHitAlg.h.

bool trkf::Track3DKalmanHitAlg::fLineSurface
private

Line surface flag.

Definition at line 114 of file Track3DKalmanHitAlg.h.

int trkf::Track3DKalmanHitAlg::fMaxChopHits
private

Maximum number of hits to chop from each end of seed.

Definition at line 117 of file Track3DKalmanHitAlg.h.

double trkf::Track3DKalmanHitAlg::fMaxSeedChiDF
private

Maximum seed track chisquare/dof.

Definition at line 118 of file Track3DKalmanHitAlg.h.

double trkf::Track3DKalmanHitAlg::fMaxTcut
private

Maximum delta ray energy in MeV for restricted dE/dx.

Definition at line 113 of file Track3DKalmanHitAlg.h.

int trkf::Track3DKalmanHitAlg::fMinSeedChopHits
private

Potentially chop seeds that exceed this length.

Definition at line 116 of file Track3DKalmanHitAlg.h.

size_t trkf::Track3DKalmanHitAlg::fMinSeedHits
private

Minimum number of hits per track seed.

Definition at line 115 of file Track3DKalmanHitAlg.h.

double trkf::Track3DKalmanHitAlg::fMinSeedSlope
private

Minimum seed slope (dx/dz).

Definition at line 119 of file Track3DKalmanHitAlg.h.

int trkf::Track3DKalmanHitAlg::fNumTrack
private

Number of tracks produced.

Definition at line 128 of file Track3DKalmanHitAlg.h.

SeedFinderAlgorithm trkf::Track3DKalmanHitAlg::fSeedFinderAlg
private

Seed finder.

Definition at line 125 of file Track3DKalmanHitAlg.h.

bool trkf::Track3DKalmanHitAlg::fSelfSeed
private

Self seed flag.

Definition at line 112 of file Track3DKalmanHitAlg.h.


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