All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | Private Attributes | List of all members
trkf::KalmanFilterAlg Class Reference

#include <KalmanFilterAlg.h>

Public Member Functions

 KalmanFilterAlg (const fhicl::ParameterSet &pset)
 Constructor. More...
 
bool getTrace () const
 Trace config parameters. More...
 
int getPlane () const
 Preferred view plane. More...
 
void setTrace (bool trace)
 Set trace config parameter. More...
 
void setPlane (int plane)
 Set preferred view plane. More...
 
bool buildTrack (const KTrack &trk, KGTrack &trg, const Propagator &prop, const Propagator::PropDirection dir, KHitContainer &hits, bool linear) const
 Make a new track. More...
 
bool smoothTrack (KGTrack &trg, KGTrack *trg1, const Propagator &prop) const
 Smooth track. More...
 
bool extendTrack (KGTrack &trg, const Propagator &prop, KHitContainer &hits) const
 Add hits to existing track. More...
 
bool fitMomentumRange (const KGTrack &trg, KETrack &tremom) const
 Estimate track momentum using range. More...
 
bool fitMomentumMS (const KGTrack &trg, const Propagator &prop, KETrack &tremom) const
 Estimate track momentum using multiple scattering. More...
 
bool fitMomentum (const KGTrack &trg, const Propagator &prop, KETrack &tremom) const
 Estimate track momentum using either range or multiple scattering. More...
 
bool updateMomentum (const KETrack &tremom, const Propagator &prop, KGTrack &trg) const
 Set track momentum at each track surface. More...
 
bool smoothTrackIter (int niter, KGTrack &trg, const Propagator &prop) const
 Iteratively smooth a track. More...
 
void cleanTrack (KGTrack &trg) const
 Clean track by removing noise hits near endpoints. More...
 

Private Attributes

bool fTrace
 Trace flag. More...
 
double fMaxPErr
 Maximum pointing error for free propagation. More...
 
double fGoodPErr
 Pointing error threshold for switching to free propagation. More...
 
double fMaxIncChisq
 Maximum incremental chisquare to accept a hit. More...
 
double fMaxSeedIncChisq
 Maximum incremental chisquare to accept a hit in seed phase. More...
 
double fMaxSmoothIncChisq
 Maximum incremental chisquare to accept a hit in smooth phase. More...
 
double fMaxEndChisq
 Maximum incremental chisquare for endpoint hit. More...
 
int fMinLHits
 Minimum number of hits to turn off linearized propagation. More...
 
double fMaxLDist
 Maximum distance for linearized propagation. More...
 
double fMaxPredDist
 Maximum prediciton distance to accept a hit. More...
 
double fMaxSeedPredDist
 Maximum prediciton distance to accept a hit in seed phase. More...
 
double fMaxPropDist
 Maximum propagation distance to candidate surface. More...
 
double fMinSortDist
 Sort low distance threshold. More...
 
double fMaxSortDist
 Sort high distance threshold. More...
 
int fMaxSamePlane
 Maximum consecutive hits in same plane. More...
 
double fGapDist
 Minimum gap distance. More...
 
int fMaxNoiseHits
 Maximum number of hits in noise cluster. More...
 
double fMinSampleDist
 Minimum sample distance (for momentum measurement). More...
 
bool fFitMomRange
 Fit momentum using range. More...
 
bool fFitMomMS
 Fit momentum using multiple scattering. More...
 
bool fGTrace
 Graphical trace flag. More...
 
double fGTraceWW
 Window width. More...
 
double fGTraceWH
 Window height. More...
 
double fGTraceXMin
 Graphical trace minimum x. More...
 
double fGTraceXMax
 Graphical trace maximum x. More...
 
std::vector< double > fGTraceZMin
 Graphical trace minimum z for each view. More...
 
std::vector< double > fGTraceZMax
 Graphical trace maximum z for each view. More...
 
int fPlane
 Preferred view plane. More...
 
std::vector< std::unique_ptr
< TCanvas > > 
fCanvases
 Graphical trace canvases. More...
 
std::vector< TVirtualPad * > fPads
 View pads in current canvas. More...
 
TVirtualPad * fInfoPad
 Information pad. More...
 
TPaveText * fMessages
 Message box. More...
 
std::map< int, TMarker * > fMarkerMap
 Markers in current canvas. More...
 

Detailed Description

Definition at line 68 of file KalmanFilterAlg.h.

Constructor & Destructor Documentation

trkf::KalmanFilterAlg::KalmanFilterAlg ( const fhicl::ParameterSet &  pset)

Constructor.

Definition at line 226 of file KalmanFilterAlg.cxx.

227  : fTrace(false)
228  , fMaxPErr(0.)
229  , fGoodPErr(0.)
230  , fMaxIncChisq(0.)
231  , fMaxSeedIncChisq(0.)
232  , fMaxSmoothIncChisq(0.)
233  , fMaxEndChisq(0.)
234  , fMinLHits(0)
235  , fMaxLDist(0.)
236  , fMaxPredDist(0.)
237  , fMaxSeedPredDist(0.)
238  , fMaxPropDist(0.)
239  , fMinSortDist(0.)
240  , fMaxSortDist(0.)
241  , fMaxSamePlane(0)
242  , fGapDist(0.)
243  , fMaxNoiseHits(0)
244  , fMinSampleDist(0.)
245  , fFitMomRange(false)
246  , fFitMomMS(false)
247  , fGTrace(false)
248  , fGTraceWW(0)
249  , fGTraceWH(0)
250  , fPlane(-1)
251  , fInfoPad(0)
252  , fMessages(0)
253 {
254  mf::LogInfo("KalmanFilterAlg") << "KalmanFilterAlg instantiated.";
255 
256  fTrace = pset.get<bool>("Trace");
257  fMaxPErr = pset.get<double>("MaxPErr");
258  fGoodPErr = pset.get<double>("GoodPErr");
259  fMaxIncChisq = pset.get<double>("MaxIncChisq");
260  fMaxSeedIncChisq = pset.get<double>("MaxSeedIncChisq");
261  fMaxSmoothIncChisq = pset.get<double>("MaxSmoothIncChisq");
262  fMaxEndChisq = pset.get<double>("MaxEndChisq");
263  fMinLHits = pset.get<int>("MinLHits");
264  fMaxLDist = pset.get<double>("MaxLDist");
265  fMaxPredDist = pset.get<double>("MaxPredDist");
266  fMaxSeedPredDist = pset.get<double>("MaxSeedPredDist");
267  fMaxPropDist = pset.get<double>("MaxPropDist");
268  fMinSortDist = pset.get<double>("MinSortDist");
269  fMaxSortDist = pset.get<double>("MaxSortDist");
270  fMaxSamePlane = pset.get<int>("MaxSamePlane");
271  fGapDist = pset.get<double>("GapDist");
272  fMaxNoiseHits = pset.get<double>("MaxNoiseHits");
273  fMinSampleDist = pset.get<double>("MinSampleDist");
274  fFitMomRange = pset.get<bool>("FitMomRange");
275  fFitMomMS = pset.get<bool>("FitMomMS");
276  fGTrace = pset.get<bool>("GTrace");
277  if (fGTrace) {
278  fGTraceWW = pset.get<double>("GTraceWW");
279  fGTraceWH = pset.get<double>("GTraceWH");
280  fGTraceXMin = pset.get<double>("GTraceXMin");
281  fGTraceXMax = pset.get<double>("GTraceXMax");
282  fGTraceZMin = pset.get<std::vector<double>>("GTraceZMin");
283  fGTraceZMax = pset.get<std::vector<double>>("GTraceZMax");
284  }
285 }
bool fGTrace
Graphical trace flag.
bool fFitMomRange
Fit momentum using range.
double fMaxPropDist
Maximum propagation distance to candidate surface.
double fGapDist
Minimum gap distance.
bool fFitMomMS
Fit momentum using multiple scattering.
double fMaxEndChisq
Maximum incremental chisquare for endpoint hit.
bool fTrace
Trace flag.
std::vector< double > fGTraceZMin
Graphical trace minimum z for each view.
double fMaxLDist
Maximum distance for linearized propagation.
double fMaxSortDist
Sort high distance threshold.
double fMaxSeedPredDist
Maximum prediciton distance to accept a hit in seed phase.
double fMaxPredDist
Maximum prediciton distance to accept a hit.
double fGTraceWH
Window height.
TVirtualPad * fInfoPad
Information pad.
double fGTraceXMin
Graphical trace minimum x.
double fGoodPErr
Pointing error threshold for switching to free propagation.
TPaveText * fMessages
Message box.
double fMaxSeedIncChisq
Maximum incremental chisquare to accept a hit in seed phase.
double fMaxSmoothIncChisq
Maximum incremental chisquare to accept a hit in smooth phase.
double fMaxIncChisq
Maximum incremental chisquare to accept a hit.
double fMinSampleDist
Minimum sample distance (for momentum measurement).
double fMinSortDist
Sort low distance threshold.
double fGTraceXMax
Graphical trace maximum x.
double fGTraceWW
Window width.
std::vector< double > fGTraceZMax
Graphical trace maximum z for each view.
int fMaxSamePlane
Maximum consecutive hits in same plane.
int fMaxNoiseHits
Maximum number of hits in noise cluster.
double fMaxPErr
Maximum pointing error for free propagation.
int fMinLHits
Minimum number of hits to turn off linearized propagation.
int fPlane
Preferred view plane.

Member Function Documentation

bool trkf::KalmanFilterAlg::buildTrack ( const KTrack trk,
KGTrack trg,
const Propagator prop,
const Propagator::PropDirection  dir,
KHitContainer hits,
bool  linear 
) const

Make a new track.

Add hits to track.

Arguments:

trk - Starting track. trg - Result global track. prop - Propagator. dir - Direction. hits - Candidate hits.

Returns: True if success.

This method makes a unidirectional Kalman fit in the specified direction, visiting each surface of the passed KHitContainer and updating the track. In case of multiple measurements on the same surface, keep (at most) the one with the smallest incremental chisquare. Any measurements that fail the incremental chisquare cut are rejected. Resolved hits (accepted or rejected) are moved to the unused list in KHitContainer. The container is sorted at the start of the method, and may be resorted during the progress of the fit.

Definition at line 310 of file KalmanFilterAlg.cxx.

316 {
317  // Direction must be forward or backward (unknown is not allowed).
318 
320  throw cet::exception("KalmanFilterAlg") << "No direction for Kalman fit.\n";
321 
322  // Set up canvas for graphical trace.
323 
324  if (fGTrace) {
325 
326  // Make a new canvas with a unique name.
327 
328  static int cnum = 0;
329  ++cnum;
330  std::ostringstream ostr;
331  ostr << "khit" << cnum;
332  fCanvases.emplace_back(
333  new TCanvas(ostr.str().c_str(), ostr.str().c_str(), fGTraceWW, fGTraceWH));
334  fPads.clear();
335  fMarkerMap.clear();
336  int nview = 3;
337  fCanvases.back()->Divide(2, nview / 2 + 1);
338 
339  // Make subpad for each view.
340 
341  for (int iview = 0; iview < nview; ++iview) {
342 
343  std::ostringstream ostr;
344  ostr << "Plane " << iview;
345 
346  fCanvases.back()->cd(iview + 1);
347  double zmin = 0.06;
348  double zmax = 0.94;
349  double xmin = 0.04;
350  double xmax = 0.95;
351  TPad* p = new TPad(ostr.str().c_str(), ostr.str().c_str(), zmin, xmin, zmax, xmax);
352  p->SetBit(kCanDelete); // Give away ownership.
353  p->Range(fGTraceZMin[iview], fGTraceXMin, fGTraceZMax[iview], fGTraceXMax);
354  p->SetFillStyle(4000); // Transparent.
355  p->Draw();
356  fPads.push_back(p);
357 
358  // Draw label.
359 
360  TText* t = new TText(zmax - 0.02, xmax - 0.03, ostr.str().c_str());
361  t->SetBit(kCanDelete); // Give away ownership.
362  t->SetTextAlign(33);
363  t->SetTextFont(42);
364  t->SetTextSize(0.04);
365  t->Draw();
366 
367  // Draw axes.
368 
369  TGaxis* pz1 =
370  new TGaxis(zmin, xmin, zmax, xmin, fGTraceZMin[iview], fGTraceZMax[iview], 510, "");
371  pz1->SetBit(kCanDelete); // Give away ownership.
372  pz1->Draw();
373 
374  TGaxis* px1 = new TGaxis(zmin, xmin, zmin, xmax, fGTraceXMin, fGTraceXMax, 510, "");
375  px1->SetBit(kCanDelete); // Give away ownership.
376  px1->Draw();
377 
378  TGaxis* pz2 =
379  new TGaxis(zmin, xmax, zmax, xmax, fGTraceZMin[iview], fGTraceZMax[iview], 510, "-");
380  pz2->SetBit(kCanDelete); // Give away ownership.
381  pz2->Draw();
382 
383  TGaxis* px2 = new TGaxis(zmax, xmin, zmax, xmax, fGTraceXMin, fGTraceXMax, 510, "L+");
384  px2->SetBit(kCanDelete); // Give away ownership.
385  px2->Draw();
386  }
387 
388  // Draw legend.
389 
390  fCanvases.back()->cd(nview + 1);
391  fInfoPad = gPad;
392  TLegend* leg = new TLegend(0.6, 0.5, 0.99, 0.99);
393  leg->SetBit(kCanDelete); // Give away ownership.
394 
395  TLegendEntry* entry = 0;
396  TMarker* m = 0;
397 
398  m = new TMarker(0., 0., 20);
399  m->SetBit(kCanDelete);
400  m->SetMarkerSize(1.2);
401  m->SetMarkerColor(kRed);
402  entry = leg->AddEntry(m, "Hits on Track", "P");
403  entry->SetBit(kCanDelete);
404 
405  m = new TMarker(0., 0., 20);
406  m->SetBit(kCanDelete);
407  m->SetMarkerSize(1.2);
408  m->SetMarkerColor(kOrange);
409  entry = leg->AddEntry(m, "Smoothed Hits on Track", "P");
410  entry->SetBit(kCanDelete);
411 
412  m = new TMarker(0., 0., 20);
413  m->SetBit(kCanDelete);
414  m->SetMarkerSize(1.2);
415  m->SetMarkerColor(kBlack);
416  entry = leg->AddEntry(m, "Available Hits", "P");
417  entry->SetBit(kCanDelete);
418 
419  m = new TMarker(0., 0., 20);
420  m->SetBit(kCanDelete);
421  m->SetMarkerSize(1.2);
422  m->SetMarkerColor(kBlue);
423  entry = leg->AddEntry(m, "Rejected Hits", "P");
424  entry->SetBit(kCanDelete);
425 
426  m = new TMarker(0., 0., 20);
427  m->SetBit(kCanDelete);
428  m->SetMarkerSize(1.2);
429  m->SetMarkerColor(kGreen);
430  entry = leg->AddEntry(m, "Removed Hits", "P");
431  entry->SetBit(kCanDelete);
432 
433  m = new TMarker(0., 0., 20);
434  m->SetBit(kCanDelete);
435  m->SetMarkerSize(1.2);
436  m->SetMarkerColor(kMagenta);
437  entry = leg->AddEntry(m, "Seed Hits", "P");
438  entry->SetBit(kCanDelete);
439 
440  m = new TMarker(0., 0., 20);
441  m->SetBit(kCanDelete);
442  m->SetMarkerSize(1.2);
443  m->SetMarkerColor(kCyan);
444  entry = leg->AddEntry(m, "Unreachable Seed Hits", "P");
445  entry->SetBit(kCanDelete);
446 
447  leg->Draw();
448 
449  // Draw message box.
450 
451  fMessages = new TPaveText(0.01, 0.01, 0.55, 0.99, "");
452  fMessages->SetBit(kCanDelete);
453  fMessages->SetTextFont(42);
454  fMessages->SetTextSize(0.04);
455  fMessages->SetTextAlign(12);
456  fMessages->Draw();
457  TText* t = fMessages->AddText("Enter buildTrack");
458  t->SetBit(kCanDelete);
459 
460  fCanvases.back()->Update();
461  }
462 
463  // Sort container using this seed track.
464 
465  hits.sort(trk, true, prop, Propagator::UNKNOWN);
466 
467  // Draw hits and populate hit->marker map.
468 
469  if (fGTrace && fCanvases.size() > 0) {
470 
471  // Loop over sorted KHitGroups.
472  // Paint sorted seed hits magenta.
473 
474  const std::list<KHitGroup>& groups = hits.getSorted();
475  for (auto const& gr : groups) {
476 
477  // Loop over hits in this group.
478 
479  const std::vector<std::shared_ptr<const KHitBase>>& phits = gr.getHits();
480  for (auto const& phit : phits) {
481  const KHitBase& hit = *phit;
482  int pl = hit.getMeasPlane();
483  if (pl >= 0 && pl < int(fPads.size())) {
484  double z = 0.;
485  double x = 0.;
486  hit_position(hit, z, x);
487  TMarker* marker = new TMarker(z, x, 20);
488  fMarkerMap[hit.getID()] = marker;
489  fPads[pl]->cd();
490  marker->SetBit(kCanDelete); // Give away ownership.
491  marker->SetMarkerSize(0.5);
492  marker->SetMarkerColor(kMagenta);
493  marker->Draw();
494  }
495  }
496  }
497 
498  // Loop over unsorted KHitGroups.
499  // Paint unsorted seed hits cyan.
500  // There should be few, if any, unsorted seed hits.
501 
502  const std::list<KHitGroup>& ugroups = hits.getUnsorted();
503  for (auto const& gr : ugroups) {
504 
505  // Loop over hits in this group.
506 
507  const std::vector<std::shared_ptr<const KHitBase>>& phits = gr.getHits();
508  for (auto const& phit : phits) {
509  const KHitBase& hit = *phit;
510  int pl = hit.getMeasPlane();
511  if (pl >= 0 && pl < int(fPads.size())) {
512  double z = 0.;
513  double x = 0.;
514  hit_position(hit, z, x);
515  TMarker* marker = new TMarker(z, x, 20);
516  fMarkerMap[hit.getID()] = marker;
517  fPads[pl]->cd();
518  marker->SetBit(kCanDelete); // Give away ownership.
519  marker->SetMarkerSize(0.5);
520  marker->SetMarkerColor(kCyan);
521  marker->Draw();
522  }
523  }
524  }
525  fCanvases.back()->Update();
526  }
527 
528  // Loop over measurements (KHitGroup) from sorted list.
529 
530  double tchisq = 0.; // Cumulative chisquare.
531  double path = 0.; // Cumulative path distance.
532  int step = 0; // Step count.
533  int nsame = 0; // Number of consecutive measurements in same plane.
534  int last_plane = -1; // Last plane.
535  bool has_pref_plane = false; // Set to true when first preferred plane hit is added to track.
536 
537  // Make a copy of the starting track, in the form of a KFitTrack,
538  // which we will update as we go.
539 
540  TrackError err;
541  trk.getSurface()->getStartingError(err);
542  KETrack tre(trk, err);
543  KFitTrack trf(tre, path, tchisq);
544 
545  // Also use the starting track as the reference track for linearized
546  // propagation, until the track is established with reasonably small
547  // errors.
548 
549  KTrack ref(trk);
550  KTrack* pref = &ref;
551 
552  mf::LogInfo log("KalmanFilterAlg");
553 
554  // Loop over measurement groups (KHitGroups).
555 
556  while (hits.getSorted().size() > 0) {
557  ++step;
558  if (fTrace) {
559  log << "Build Step " << step << "\n";
560  log << "KGTrack has " << trg.numHits() << " hits.\n";
561  log << trf;
562  }
563 
564  // Get an iterator for the next KHitGroup.
565 
566  std::list<KHitGroup>::iterator it;
567  if (dir == Propagator::FORWARD)
568  it = hits.getSorted().begin();
569  else {
570  it = hits.getSorted().end();
571  --it;
572  }
573  const KHitGroup& gr = *it;
574 
575  if (fTrace) {
576  double path_est = gr.getPath();
577  log << "Next surface: " << *(gr.getSurface()) << "\n";
578  log << " Estimated path length of hit group = " << path_est << "\n";
579  }
580 
581  // Get the next prediction surface. If this KHitGroup is on the
582  // preferred plane, use that as the prediction surface.
583  // Otherwise, use the current track surface as the prediction
584  // surface.
585 
586  std::shared_ptr<const Surface> psurf = trf.getSurface();
587  if (gr.getPlane() < 0) throw cet::exception("KalmanFilterAlg") << "negative plane?\n";
588  if (fPlane < 0 || gr.getPlane() < 0 || fPlane == gr.getPlane()) psurf = gr.getSurface();
589 
590  // Propagate track to the prediction surface.
591 
592  std::optional<double> dist = prop.noise_prop(trf, psurf, Propagator::UNKNOWN, true, pref);
593  if (dist && std::abs(*dist) > fMaxPropDist) dist = std::nullopt;
594  double ds = 0.;
595 
596  if (dist) {
597 
598  // Propagation succeeded.
599  // Update cumulative path distance and track status.
600 
601  ds = *dist;
602  path += ds;
603  trf.setPath(path);
604  if (dir == Propagator::FORWARD)
605  trf.setStat(KFitTrack::FORWARD_PREDICTED);
606  else {
607  trf.setStat(KFitTrack::BACKWARD_PREDICTED);
608  }
609  if (fTrace) {
610  log << "After propagation\n";
611  log << " Incremental propagation distance = " << ds << "\n";
612  log << " Path length of prediction surface = " << path << "\n";
613  log << "KGTrack has " << trg.numHits() << " hits.\n";
614  log << trf;
615  }
616 
617  // Loop over measurements in this group.
618 
619  const std::vector<std::shared_ptr<const KHitBase>>& hits = gr.getHits();
620  double best_chisq = 0.;
621  std::shared_ptr<const KHitBase> best_hit;
622  for (std::vector<std::shared_ptr<const KHitBase>>::const_iterator ihit = hits.begin();
623  ihit != hits.end();
624  ++ihit) {
625  const KHitBase& hit = **ihit;
626 
627  // Turn this hit blue.
628 
629  if (fGTrace && fCanvases.size() > 0) {
630  auto marker_it = fMarkerMap.find(hit.getID());
631  if (marker_it != fMarkerMap.end()) {
632  TMarker* marker = marker_it->second;
633  marker->SetMarkerColor(kBlue);
634  }
635  //fCanvases.back()->Update();
636  }
637 
638  // Update predction using current track hypothesis and get
639  // incremental chisquare.
640 
641  if (hit.predict(trf, prop)) {
642  double chisq = hit.getChisq();
643  double preddist = hit.getPredDistance();
644  if (fTrace) {
645  log << "Trying Hit.\n"
646  << hit << "\nchisq = " << chisq << "\n"
647  << "prediction distance = " << preddist << "\n";
648  }
649  if ((!has_pref_plane || abs(preddist) < fMaxSeedPredDist) &&
650  (best_hit.get() == 0 || chisq < best_chisq)) {
651  best_chisq = chisq;
652  if (chisq < fMaxSeedIncChisq) best_hit = *ihit;
653  }
654  }
655  }
656  if (fTrace) {
657  log << "Best hit incremental chisquare = " << best_chisq << "\n";
658  if (best_hit.get() != 0) {
659  log << "Hit after prediction\n";
660  log << *best_hit;
661  }
662  else
663  log << "No hit passed chisquare cut.\n";
664  }
665  if (fGTrace && fCanvases.size() > 0) fCanvases.back()->Update();
666 
667  // If we found a best measurement, and if the incremental
668  // chisquare passes the cut, add it to the track and update
669  // fit information.
670 
671  bool update_ok = false;
672  if (best_hit.get() != 0) {
673  KFitTrack trf0(trf);
674  best_hit->update(trf);
675  update_ok = trf.isValid();
676  if (!update_ok) trf = trf0;
677  }
678  if (update_ok) {
679  ds += best_hit->getPredDistance();
680  tchisq += best_chisq;
681  trf.setChisq(tchisq);
682  if (dir == Propagator::FORWARD)
683  trf.setStat(KFitTrack::FORWARD);
684  else {
685  trf.setStat(KFitTrack::BACKWARD);
686  }
687 
688  // If the pointing error got too large, and there is no
689  // reference track, quit.
690 
691  if (pref == 0 && trf.PointingError() > fMaxPErr) {
692  if (fTrace) log << "Quitting because pointing error got too large.\n";
693  break;
694  }
695 
696  // Test number of consecutive measurements in the same plane.
697 
698  if (gr.getPlane() >= 0) {
699  if (gr.getPlane() == last_plane)
700  ++nsame;
701  else {
702  nsame = 1;
703  last_plane = gr.getPlane();
704  }
705  }
706  else {
707  nsame = 0;
708  last_plane = -1;
709  }
710  if (nsame <= fMaxSamePlane) {
711 
712  // Turn best hit red.
713 
714  if (fGTrace && fCanvases.size() > 0) {
715  int pl = best_hit->getMeasPlane();
716  if (pl >= 0 && pl < int(fPads.size())) {
717  auto marker_it = fMarkerMap.find(best_hit->getID());
718  if (marker_it != fMarkerMap.end()) {
719  TMarker* marker = marker_it->second;
720  marker->SetMarkerColor(kRed);
721 
722  // Redraw marker so that it will be on top.
723 
724  fPads[pl]->cd();
725  marker->Draw();
726  }
727  }
728  fCanvases.back()->Update();
729  }
730 
731  // Make a KHitTrack and add it to the KGTrack.
732 
733  KHitTrack trh(trf, best_hit);
734  trg.addTrack(trh);
735  if (fPlane == gr.getPlane()) has_pref_plane = true;
736 
737  // Decide if we want to kill the reference track.
738 
739  if (!linear && pref != 0 && int(trg.numHits()) >= fMinLHits &&
740  (trf.PointingError() < fGoodPErr || path > fMaxLDist)) {
741  pref = 0;
742  if (fTrace) log << "Killing reference track.\n";
743  }
744 
745  if (fTrace) {
746  log << "After update\n";
747  log << "KGTrack has " << trg.numHits() << " hits.\n";
748  log << trf;
749  }
750  }
751  }
752  }
753 
754  // The current KHitGroup is now resolved.
755  // Move it to unused list.
756 
757  hits.getUnused().splice(hits.getUnused().end(), hits.getSorted(), it);
758 
759  // If the propagation distance was the wrong direction, resort the measurements.
760 
761  if (pref == 0 && dist &&
762  ((dir == Propagator::FORWARD && (ds < fMinSortDist || ds > fMaxSortDist)) ||
763  (dir == Propagator::BACKWARD && (-ds < fMinSortDist || -ds > fMaxSortDist)))) {
764  if (fTrace) log << "Resorting measurements.\n";
765  hits.sort(trf, true, prop, dir);
766  }
767  }
768 
769  // Clean track.
770 
771  cleanTrack(trg);
772 
773  // Set the fit status of the last added KHitTrack to optimal and get
774  // the final chisquare and path length.
775 
776  double fchisq = 0.;
777  path = 0.;
778  if (trg.isValid()) {
779  path = trg.endTrack().getPath() - trg.startTrack().getPath();
780  if (dir == Propagator::FORWARD) {
781  trg.endTrack().setStat(KFitTrack::OPTIMAL);
782  fchisq = trg.endTrack().getChisq();
783  }
784  else {
785  trg.startTrack().setStat(KFitTrack::OPTIMAL);
786  fchisq = trg.startTrack().getChisq();
787  }
788  }
789 
790  // Summary.
791 
792  log << "KalmanFilterAlg build track summary.\n"
793  << "Build direction = " << (dir == Propagator::FORWARD ? "FORWARD" : "BACKWARD") << "\n"
794  << "Track has " << trg.numHits() << " hits.\n"
795  << "Track length = " << path << "\n"
796  << "Track chisquare = " << fchisq << "\n";
797 
798  // Done. Return success if we added at least one measurement.
799 
800  bool ok = trg.numHits() > 0;
801 
802  if (fGTrace && fCanvases.size() > 0) {
803  TText* t = 0;
804  if (ok)
805  t = fMessages->AddText("Exit buildTrack, status success");
806  else
807  t = fMessages->AddText("Exit buildTrack, status fail");
808  t->SetBit(kCanDelete);
809  fInfoPad->Modified();
810  fCanvases.back()->Update();
811  }
812 
813  return ok;
814 }
bool fGTrace
Graphical trace flag.
process_name opflash particleana ie ie ie z
double fMaxPropDist
Maximum propagation distance to candidate surface.
std::map< int, TMarker * > fMarkerMap
Markers in current canvas.
process_name opflash particleana ie x
KSymMatrix< 5 >::type TrackError
Track error matrix, dimension 5x5.
bool fTrace
Trace flag.
EResult err(const char *call)
pdgs p
Definition: selectors.fcl:22
void cleanTrack(KGTrack &trg) const
Clean track by removing noise hits near endpoints.
process_name hit
Definition: cheaterreco.fcl:51
tuple m
now if test mode generate materials, CRT shell, world, gdml header else just generate CRT shell for u...
std::vector< double > fGTraceZMin
Graphical trace minimum z for each view.
double fMaxLDist
Maximum distance for linearized propagation.
double fMaxSortDist
Sort high distance threshold.
double fMaxSeedPredDist
Maximum prediciton distance to accept a hit in seed phase.
process_name pandoraGausCryo1 vertexChargeCryo1 vertexStubCryo1 xmin
BEGIN_PROLOG triggeremu_data_config_icarus settings PMTADCthresholds sequence::icarus_stage0_multiTPC_TPC physics sequence::icarus_stage0_EastHits_TPC physics sequence::icarus_stage0_WestHits_TPC physics producers purityana0 caloskimCalorimetryCryoE physics caloskimCalorimetryCryoW physics path
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
T abs(T value)
double fGTraceWH
Window height.
TVirtualPad * fInfoPad
Information pad.
std::vector< std::unique_ptr< TCanvas > > fCanvases
Graphical trace canvases.
double fGTraceXMin
Graphical trace minimum x.
double fGoodPErr
Pointing error threshold for switching to free propagation.
TPaveText * fMessages
Message box.
double fMaxSeedIncChisq
Maximum incremental chisquare to accept a hit in seed phase.
tuple dir
Definition: dropbox.py:28
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
double fGTraceXMax
Graphical trace maximum x.
std::vector< TVirtualPad * > fPads
View pads in current canvas.
double fGTraceWW
Window width.
std::vector< double > fGTraceZMax
Graphical trace maximum z for each view.
int fMaxSamePlane
Maximum consecutive hits in same plane.
double fMaxPErr
Maximum pointing error for free propagation.
int fMinLHits
Minimum number of hits to turn off linearized propagation.
int fPlane
Preferred view plane.
void trkf::KalmanFilterAlg::cleanTrack ( KGTrack trg) const

Clean track by removing noise hits near endpoints.

Clean track by removing noise hits near endpoints.

Arguments:

trg - Track.

Definition at line 2107 of file KalmanFilterAlg.cxx.

2108 {
2109  // Get hold of a modifiable track map from trg.
2110 
2111  std::multimap<double, KHitTrack>& trackmap = trg.getTrackMap();
2112 
2113  // Do an indefinite loop until we no longer find any dirt.
2114 
2115  bool done = false;
2116  while (!done) {
2117 
2118  // If track map has fewer than fMaxNoiseHits tracks, then this is a
2119  // noise track. Clear the map, making the whole track invalid.
2120 
2121  int ntrack = trackmap.size();
2122  if (ntrack <= fMaxNoiseHits) {
2123  for (auto const& trackmap_entry : trackmap) {
2124  const KHitTrack& trh = trackmap_entry.second;
2125  const KHitBase& hit = *(trh.getHit());
2126  auto marker_it = fMarkerMap.find(hit.getID());
2127  if (marker_it != fMarkerMap.end()) {
2128  TMarker* marker = marker_it->second;
2129  marker->SetMarkerColor(kGreen);
2130  }
2131  }
2132  trackmap.clear();
2133  done = true;
2134  break;
2135  }
2136 
2137  // Make sure the first two and last two tracks belong to different
2138  // views. If not, remove the first or last track.
2139 
2140  if (ntrack >= 2) {
2141 
2142  // Check start.
2143 
2144  std::multimap<double, KHitTrack>::iterator it = trackmap.begin();
2145  const KHitTrack& trh1a = (*it).second;
2146  const KHitBase& hit1a = *(trh1a.getHit());
2147  int plane1 = hit1a.getMeasPlane();
2148  double chisq1 = hit1a.getChisq();
2149  ++it;
2150  const KHitTrack& trh2a = (*it).second;
2151  const KHitBase& hit2a = *(trh2a.getHit());
2152  int plane2 = hit2a.getMeasPlane();
2153  double chisq2 = hit2a.getChisq();
2154  if ((plane1 >= 0 && plane2 >= 0 && plane1 == plane2) || chisq1 > fMaxEndChisq ||
2155  chisq2 > fMaxEndChisq) {
2156  auto marker_it = fMarkerMap.find(hit1a.getID());
2157  if (marker_it != fMarkerMap.end()) {
2158  TMarker* marker = marker_it->second;
2159  marker->SetMarkerColor(kGreen);
2160  }
2161  trackmap.erase(trackmap.begin(), it);
2162  done = false;
2163  continue;
2164  }
2165 
2166  // Check end.
2167 
2168  it = trackmap.end();
2169  --it;
2170  const KHitTrack& trh1b = (*it).second;
2171  const KHitBase& hit1b = *(trh1b.getHit());
2172  plane1 = hit1b.getMeasPlane();
2173  chisq1 = hit1b.getChisq();
2174  --it;
2175  const KHitTrack& trh2b = (*it).second;
2176  const KHitBase& hit2b = *(trh2b.getHit());
2177  plane2 = hit2b.getMeasPlane();
2178  chisq2 = hit2b.getChisq();
2179  if ((plane1 >= 0 && plane2 >= 0 && plane1 == plane2) || chisq1 > fMaxEndChisq ||
2180  chisq2 > fMaxEndChisq) {
2181  ++it;
2182  auto marker_it = fMarkerMap.find(hit1b.getID());
2183  if (marker_it != fMarkerMap.end()) {
2184  TMarker* marker = marker_it->second;
2185  marker->SetMarkerColor(kGreen);
2186  }
2187  trackmap.erase(it, trackmap.end());
2188  done = false;
2189  continue;
2190  }
2191  }
2192 
2193  // Loop over successive pairs of elements of track map. Look for
2194  // adjacent elements with distance separation greater than
2195  // fGapDist.
2196 
2197  std::multimap<double, KHitTrack>::iterator it = trackmap.begin();
2198  std::multimap<double, KHitTrack>::iterator jt = trackmap.end();
2199  int nb = 0; // Number of elements from begin to jt.
2200  int ne = ntrack; // Number of elements it to end.
2201  bool found_noise = false;
2202  for (; it != trackmap.end(); ++it, ++nb, --ne) {
2203  if (jt == trackmap.end())
2204  jt = trackmap.begin();
2205  else {
2206  if (nb < 1)
2207  throw cet::exception("KalmanFilterAlg")
2208  << "KalmanFilterAlg::cleanTrack(): nb not positive\n";
2209  if (ne < 1)
2210  throw cet::exception("KalmanFilterAlg")
2211  << "KalmanFilterAlg::cleanTrack(): ne not positive\n";
2212 
2213  double disti = (*it).first;
2214  double distj = (*jt).first;
2215  double sep = disti - distj;
2216  if (sep < 0.)
2217  throw cet::exception("KalmanFilterAlg")
2218  << "KalmanFilterAlg::fitMomentumRange(): negative separation\n";
2219 
2220  if (sep > fGapDist) {
2221 
2222  // Found a gap. See if we want to trim track.
2223 
2224  if (nb <= fMaxNoiseHits) {
2225 
2226  // Trim front.
2227 
2228  found_noise = true;
2229  for (auto jt = trackmap.begin(); jt != it; ++jt) {
2230  const KHitTrack& trh = jt->second;
2231  const KHitBase& hit = *(trh.getHit());
2232  auto marker_it = fMarkerMap.find(hit.getID());
2233  if (marker_it != fMarkerMap.end()) {
2234  TMarker* marker = marker_it->second;
2235  marker->SetMarkerColor(kGreen);
2236  }
2237  }
2238  trackmap.erase(trackmap.begin(), it);
2239  break;
2240  }
2241  if (ne <= fMaxNoiseHits) {
2242 
2243  // Trim back.
2244 
2245  found_noise = true;
2246  for (auto jt = it; jt != trackmap.end(); ++jt) {
2247  const KHitTrack& trh = jt->second;
2248  const KHitBase& hit = *(trh.getHit());
2249  auto marker_it = fMarkerMap.find(hit.getID());
2250  if (marker_it != fMarkerMap.end()) {
2251  TMarker* marker = marker_it->second;
2252  marker->SetMarkerColor(kGreen);
2253  }
2254  }
2255  trackmap.erase(it, trackmap.end());
2256  break;
2257  }
2258  }
2259  ++jt;
2260  }
2261  }
2262 
2263  // Decide if we are done.
2264 
2265  done = !found_noise;
2266  }
2267  if (fGTrace && fCanvases.size() > 0) fCanvases.back()->Update();
2268 }
bool fGTrace
Graphical trace flag.
double fGapDist
Minimum gap distance.
std::map< int, TMarker * > fMarkerMap
Markers in current canvas.
double fMaxEndChisq
Maximum incremental chisquare for endpoint hit.
process_name hit
Definition: cheaterreco.fcl:51
std::vector< std::unique_ptr< TCanvas > > fCanvases
Graphical trace canvases.
int fMaxNoiseHits
Maximum number of hits in noise cluster.
bool trkf::KalmanFilterAlg::extendTrack ( KGTrack trg,
const Propagator prop,
KHitContainer hits 
) const

Add hits to existing track.

Add hits to existing track.

Arguments:

trg - Track to extend. prop - Propagator. hits - Hit collection to choose hits from.

This method extends a KGTrack by adding hits. The KGTrack must have previously been produced by a unidirectional Kalman fit (it should be optimal at one end). This method finds the optimal end and extends the track in that direction. If any hits are added, the originally optimal end has its status reset to forward or backward, and the new endpoint is optimal. In any case, the final result is unidirectionally fit KGTrack.

Definition at line 1158 of file KalmanFilterAlg.cxx.

1159 {
1160  // Default result failure.
1161 
1162  bool result = false;
1163 
1164  if (fGTrace && fCanvases.size() > 0) {
1165  TText* t = fMessages->AddText("Enter extendTrack");
1166  t->SetBit(kCanDelete);
1167  fInfoPad->Modified();
1168  fCanvases.back()->Update();
1169  }
1170 
1171  // Remember the original number of measurement.
1172 
1173  unsigned int nhits0 = trg.numHits();
1174 
1175  // It is an error if the KGTrack is not valid.
1176 
1177  if (trg.isValid()) {
1178  mf::LogInfo log("KalmanFilterAlg");
1179 
1180  // Examine the track endpoints and figure out which end of the
1181  // track to extend. The track is always extended from the optimal
1182  // end. It is an error if neither end point is optimal, or both
1183  // endoints are optimal. Reset the status of the optimal, and
1184  // make a copy of the starting track fit. Also get starting path
1185  // and chisquare.
1186 
1187  KHitTrack& trh0 = trg.startTrack();
1188  KHitTrack& trh1 = trg.endTrack();
1189  KFitTrack::FitStatus stat0 = trh0.getStat();
1190  KFitTrack::FitStatus stat1 = trh1.getStat();
1191  bool dofit = false;
1193  KFitTrack trf;
1194  double path = 0.;
1195  double tchisq = 0.;
1196 
1197  if (stat0 == KFitTrack::OPTIMAL) {
1198  if (stat1 == KFitTrack::OPTIMAL) {
1199 
1200  // Both ends optimal (do nothing, return failure).
1201 
1202  dofit = false;
1203  result = false;
1204  return result;
1205  }
1206  else {
1207 
1208  // Start optimal. Extend backward.
1209 
1210  dofit = true;
1211  dir = Propagator::BACKWARD;
1212  trh0.setStat(KFitTrack::BACKWARD);
1213  trf = trh0;
1214  path = trh0.getPath();
1215  tchisq = trh0.getChisq();
1216  }
1217  }
1218  else {
1219  if (stat1 == KFitTrack::OPTIMAL) {
1220 
1221  // End optimal. Extend forward.
1222 
1223  dofit = true;
1224  dir = Propagator::FORWARD;
1225  trh1.setStat(KFitTrack::FORWARD);
1226  trf = trh1;
1227  path = trh1.getPath();
1228  tchisq = trh1.getChisq();
1229 
1230  // Make sure forward extend track momentum is over some
1231  // minimum value.
1232 
1233  if (trf.getVector()(4) > 5.) {
1234  trf.getVector()(4) = 5.;
1235  trf.getError()(4, 4) = 5.;
1236  }
1237  }
1238  else {
1239 
1240  // Neither end optimal (do nothing, return failure).
1241 
1242  dofit = false;
1243  result = false;
1244  return result;
1245  }
1246  }
1247  if (dofit) {
1248 
1249  // Sort hit container using starting track.
1250 
1251  hits.sort(trf, true, prop, dir);
1252 
1253  // Draw and add hits in hit->marker map that are not already there.
1254 
1255  if (fGTrace && fCanvases.size() > 0) {
1256 
1257  // Loop over sorted KHitGroups.
1258  // Paint sorted hits black.
1259 
1260  const std::list<KHitGroup>& groups = hits.getSorted();
1261  for (auto const& gr : groups) {
1262 
1263  // Loop over hits in this group.
1264 
1265  const std::vector<std::shared_ptr<const KHitBase>>& phits = gr.getHits();
1266  for (auto const& phit : phits) {
1267  const KHitBase& hit = *phit;
1268  int pl = hit.getMeasPlane();
1269  if (pl >= 0 && pl < int(fPads.size())) {
1270 
1271  // Is this hit already in map?
1272 
1273  TMarker* marker = 0;
1274  auto marker_it = fMarkerMap.find(hit.getID());
1275  if (marker_it == fMarkerMap.end()) {
1276  double z = 0.;
1277  double x = 0.;
1278  hit_position(hit, z, x);
1279  marker = new TMarker(z, x, 20);
1280  fMarkerMap[hit.getID()] = marker;
1281  fPads[pl]->cd();
1282  marker->SetBit(kCanDelete); // Give away ownership.
1283  marker->SetMarkerSize(0.5);
1284  marker->Draw();
1285  }
1286  else
1287  marker = marker_it->second;
1288  marker->SetMarkerColor(kBlack);
1289  }
1290  }
1291  }
1292 
1293  // Loop over unsorted KHitGroups.
1294  // Paint unsorted hits blue.
1295 
1296  const std::list<KHitGroup>& ugroups = hits.getUnsorted();
1297  for (auto const& gr : ugroups) {
1298 
1299  // Loop over hits in this group.
1300 
1301  const std::vector<std::shared_ptr<const KHitBase>>& phits = gr.getHits();
1302  for (auto const& phit : phits) {
1303  const KHitBase& hit = *phit;
1304  int pl = hit.getMeasPlane();
1305  if (pl >= 0 && pl < int(fPads.size())) {
1306 
1307  // Is this hit already in map?
1308 
1309  TMarker* marker = 0;
1310  auto marker_it = fMarkerMap.find(hit.getID());
1311  if (marker_it == fMarkerMap.end()) {
1312  double z = 0.;
1313  double x = 0.;
1314  hit_position(hit, z, x);
1315  marker = new TMarker(z, x, 20);
1316  fMarkerMap[hit.getID()] = marker;
1317  fPads[pl]->cd();
1318  marker->SetBit(kCanDelete); // Give away ownership.
1319  marker->SetMarkerSize(0.5);
1320  marker->Draw();
1321  }
1322  else
1323  marker = marker_it->second;
1324  marker->SetMarkerColor(kBlue);
1325  }
1326  }
1327  }
1328  fCanvases.back()->Update();
1329  }
1330 
1331  //std::cout << "extendTrack: marker map has " << fMarkerMap.size() << " entries." << std::endl;
1332 
1333  // Extend loop starts here.
1334 
1335  int step = 0;
1336  int nsame = 0;
1337  int last_plane = -1;
1338  while (hits.getSorted().size() > 0) {
1339  ++step;
1340  if (fTrace) {
1341  log << "Extend Step " << step << "\n";
1342  log << "KGTrack has " << trg.numHits() << " hits.\n";
1343  log << trf;
1344  }
1345 
1346  // Get an iterator for the next KHitGroup.
1347 
1348  std::list<KHitGroup>::iterator it;
1349  switch (dir) {
1350  case Propagator::FORWARD: it = hits.getSorted().begin(); break;
1351  case Propagator::BACKWARD:
1352  it = hits.getSorted().end();
1353  --it;
1354  break;
1355  default:
1356  throw cet::exception("KalmanFilterAlg")
1357  << "KalmanFilterAlg::extendTrack(): invalid direction\n";
1358  } // switch
1359  const KHitGroup& gr = *it;
1360 
1361  if (fTrace) {
1362  double path_est = gr.getPath();
1363  log << "Next surface: " << *(gr.getSurface()) << "\n";
1364  log << " Estimated path length of hit group = " << path_est << "\n";
1365  }
1366 
1367  // Get the next prediction surface. If this KHitGroup is on the
1368  // preferred plane, use that as the prediction surface.
1369  // Otherwise, use the current track surface as the prediction
1370  // surface.
1371 
1372  std::shared_ptr<const Surface> psurf = trf.getSurface();
1373  if (gr.getPlane() < 0)
1374  throw cet::exception("KalmanFilterAlg")
1375  << "KalmanFilterAlg::extendTrack(): negative plane?\n";
1376  if (fPlane < 0 || gr.getPlane() < 0 || fPlane == gr.getPlane()) psurf = gr.getSurface();
1377 
1378  // Propagate track to the prediction surface.
1379 
1380  std::optional<double> dist = prop.noise_prop(trf, psurf, Propagator::UNKNOWN, true);
1381  if (dist && std::abs(*dist) > fMaxPropDist) dist = std::nullopt;
1382  double ds = 0.;
1383 
1384  if (dist) {
1385 
1386  // Propagation succeeded.
1387  // Update cumulative path distance and track status.
1388 
1389  ds = *dist;
1390  path += ds;
1391  trf.setPath(path);
1392  if (dir == Propagator::FORWARD)
1393  trf.setStat(KFitTrack::FORWARD_PREDICTED);
1394  else {
1395  trf.setStat(KFitTrack::BACKWARD_PREDICTED);
1396  }
1397  if (fTrace) {
1398  log << "After propagation\n";
1399  log << " Incremental distance = " << ds << "\n";
1400  log << " Track path length = " << path << "\n";
1401  log << "KGTrack has " << trg.numHits() << " hits.\n";
1402  log << trf;
1403  }
1404 
1405  // Loop over measurements in this group.
1406 
1407  const std::vector<std::shared_ptr<const KHitBase>>& hits = gr.getHits();
1408  double best_chisq = 0.;
1409  std::shared_ptr<const KHitBase> best_hit;
1410  for (std::vector<std::shared_ptr<const KHitBase>>::const_iterator ihit = hits.begin();
1411  ihit != hits.end();
1412  ++ihit) {
1413  const KHitBase& hit = **ihit;
1414 
1415  // Turn this hit blue.
1416 
1417  if (fGTrace && fCanvases.size() > 0) {
1418  auto marker_it = fMarkerMap.find(hit.getID());
1419  if (marker_it != fMarkerMap.end()) {
1420  TMarker* marker = marker_it->second;
1421  marker->SetMarkerColor(kBlue);
1422  }
1423  //fCanvases.back()->Update();
1424  }
1425 
1426  // Update predction using current track hypothesis and get
1427  // incremental chisquare.
1428 
1429  bool ok = hit.predict(trf, prop);
1430  if (ok) {
1431  double chisq = hit.getChisq();
1432  double preddist = hit.getPredDistance();
1433  if (abs(preddist) < fMaxPredDist && (best_hit.get() == 0 || chisq < best_chisq)) {
1434  best_chisq = chisq;
1435  if (chisq < fMaxIncChisq) best_hit = *ihit;
1436  }
1437  }
1438  }
1439  if (fTrace) {
1440  log << "Best hit incremental chisquare = " << best_chisq << "\n";
1441  if (best_hit.get() != 0) {
1442  log << "Hit after prediction\n";
1443  log << *best_hit;
1444  }
1445  else
1446  log << "No hit passed chisquare cut.\n";
1447  }
1448  if (fGTrace && fCanvases.size() > 0) fCanvases.back()->Update();
1449 
1450  // If we found a best measurement, and if the incremental
1451  // chisquare passes the cut, add it to the track and update
1452  // fit information.
1453 
1454  bool update_ok = false;
1455  if (best_hit.get() != 0) {
1456  KFitTrack trf0(trf);
1457  best_hit->update(trf);
1458  update_ok = trf.isValid();
1459  if (!update_ok) trf = trf0;
1460  }
1461  if (update_ok) {
1462  ds += best_hit->getPredDistance();
1463  tchisq += best_chisq;
1464  trf.setChisq(tchisq);
1465  if (dir == Propagator::FORWARD)
1466  trf.setStat(KFitTrack::FORWARD);
1467  else {
1468  trf.setStat(KFitTrack::BACKWARD);
1469  }
1470 
1471  // If the pointing error got too large, quit.
1472 
1473  if (trf.PointingError() > fMaxPErr) {
1474  if (fTrace) log << "Quitting because pointing error got too large.\n";
1475  break;
1476  }
1477 
1478  // Test number of consecutive measurements in the same plane.
1479 
1480  if (gr.getPlane() >= 0) {
1481  if (gr.getPlane() == last_plane)
1482  ++nsame;
1483  else {
1484  nsame = 1;
1485  last_plane = gr.getPlane();
1486  }
1487  }
1488  else {
1489  nsame = 0;
1490  last_plane = -1;
1491  }
1492  if (nsame <= fMaxSamePlane) {
1493 
1494  // Turn best hit red.
1495 
1496  if (fGTrace && fCanvases.size() > 0) {
1497  int pl = best_hit->getMeasPlane();
1498  if (pl >= 0 && pl < int(fPads.size())) {
1499  auto marker_it = fMarkerMap.find(best_hit->getID());
1500  if (marker_it != fMarkerMap.end()) {
1501  TMarker* marker = marker_it->second;
1502  marker->SetMarkerColor(kRed);
1503 
1504  // Redraw marker so that it will be on top.
1505 
1506  fPads[pl]->cd();
1507  marker->Draw();
1508  }
1509  }
1510  fCanvases.back()->Update();
1511  }
1512 
1513  // Make a KHitTrack and add it to the KGTrack.
1514 
1515  KHitTrack trh(trf, best_hit);
1516  trg.addTrack(trh);
1517 
1518  if (fTrace) {
1519  log << "After update\n";
1520  log << "KGTrack has " << trg.numHits() << " hits.\n";
1521  log << trf;
1522  }
1523  }
1524  }
1525  }
1526 
1527  // The current KHitGroup is now resolved.
1528  // Move it to unused list.
1529 
1530  hits.getUnused().splice(hits.getUnused().end(), hits.getSorted(), it);
1531 
1532  // If the propagation distance was the wrong direction, resort the measurements.
1533 
1534  if (dist && ((dir == Propagator::FORWARD && (ds < fMinSortDist || ds > fMaxSortDist)) ||
1535  (dir == Propagator::BACKWARD && (-ds < fMinSortDist || -ds > fMaxSortDist)))) {
1536  if (fTrace) log << "Resorting measurements.\n";
1537  hits.sort(trf, true, prop, dir);
1538  }
1539  }
1540  }
1541 
1542  // Clean track.
1543 
1544  cleanTrack(trg);
1545 
1546  // Set the fit status of the last added KHitTrack to optimal and
1547  // get the final chisquare and path length.
1548 
1549  double fchisq = 0.;
1550  path = 0.;
1551  if (trg.isValid()) {
1552  path = trg.endTrack().getPath() - trg.startTrack().getPath();
1553  switch (dir) {
1554  case Propagator::FORWARD:
1555  trg.endTrack().setStat(KFitTrack::OPTIMAL);
1556  fchisq = trg.endTrack().getChisq();
1557  break;
1558  case Propagator::BACKWARD:
1559  trg.startTrack().setStat(KFitTrack::OPTIMAL);
1560  fchisq = trg.startTrack().getChisq();
1561  break;
1562  default:
1563  throw cet::exception("KalmanFilterAlg")
1564  << "KalmanFilterAlg::extendTrack(): invalid direction [II]\n";
1565  } // switch
1566  }
1567 
1568  // Summary.
1569 
1570  log << "KalmanFilterAlg extend track summary.\n"
1571  << "Extend direction = " << (dir == Propagator::FORWARD ? "FORWARD" : "BACKWARD") << "\n"
1572  << "Track has " << trg.numHits() << " hits.\n"
1573  << "Track length = " << path << "\n"
1574  << "Track chisquare = " << fchisq << "\n";
1575  // log << trg << "\n";
1576  }
1577 
1578  // Done.
1579 
1580  result = (trg.numHits() > nhits0);
1581 
1582  if (fGTrace && fCanvases.size() > 0) {
1583  TText* t = 0;
1584  if (result)
1585  t = fMessages->AddText("Exit extendTrack, status success");
1586  else
1587  t = fMessages->AddText("Exit extendTrack, status fail");
1588  t->SetBit(kCanDelete);
1589  fInfoPad->Modified();
1590  fCanvases.back()->Update();
1591  }
1592 
1593  return result;
1594 }
bool fGTrace
Graphical trace flag.
process_name opflash particleana ie ie ie z
double fMaxPropDist
Maximum propagation distance to candidate surface.
std::map< int, TMarker * > fMarkerMap
Markers in current canvas.
process_name opflash particleana ie x
bool fTrace
Trace flag.
void cleanTrack(KGTrack &trg) const
Clean track by removing noise hits near endpoints.
FitStatus
Fit status enum.
Definition: KFitTrack.h:41
process_name hit
Definition: cheaterreco.fcl:51
double fMaxSortDist
Sort high distance threshold.
double fMaxPredDist
Maximum prediciton distance to accept a hit.
BEGIN_PROLOG triggeremu_data_config_icarus settings PMTADCthresholds sequence::icarus_stage0_multiTPC_TPC physics sequence::icarus_stage0_EastHits_TPC physics sequence::icarus_stage0_WestHits_TPC physics producers purityana0 caloskimCalorimetryCryoE physics caloskimCalorimetryCryoW physics path
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
T abs(T value)
TVirtualPad * fInfoPad
Information pad.
std::vector< std::unique_ptr< TCanvas > > fCanvases
Graphical trace canvases.
TPaveText * fMessages
Message box.
double fMaxIncChisq
Maximum incremental chisquare to accept a hit.
tuple dir
Definition: dropbox.py:28
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
std::vector< TVirtualPad * > fPads
View pads in current canvas.
int fMaxSamePlane
Maximum consecutive hits in same plane.
double fMaxPErr
Maximum pointing error for free propagation.
PropDirection
Propagation direction enum.
Definition: Propagator.h:94
int fPlane
Preferred view plane.
bool trkf::KalmanFilterAlg::fitMomentum ( const KGTrack trg,
const Propagator prop,
KETrack tremom 
) const

Estimate track momentum using either range or multiple scattering.

Estimate track momentum using either range or multiple scattering.

Arguments:

trg - Global track whose momentum is to be updated. prop - Propagator. tremom - Track with momentum estimate.

Returns: True if success.

Definition at line 1860 of file KalmanFilterAlg.cxx.

1863 {
1864  mf::LogInfo log("KalmanFilterAlg");
1865  double invp_range = 0.;
1866  double invp_ms = 0.;
1867 
1868  // Get multiple scattering momentum estimate.
1869 
1870  KETrack tremom_ms;
1871  bool ok_ms = false;
1872  if (fFitMomMS) {
1873  ok_ms = fitMomentumMS(trg, prop, tremom_ms);
1874  if (ok_ms) {
1875  KGTrack trg_ms(trg);
1876  ok_ms = updateMomentum(tremom_ms, prop, trg_ms);
1877  if (ok_ms) {
1878  invp_ms = trg_ms.startTrack().getVector()(4);
1879  double var_invp = trg_ms.startTrack().getError()(4, 4);
1880  double p = 0.;
1881  if (invp_ms != 0.) p = 1. / invp_ms;
1882  double err_p = p * p * std::sqrt(var_invp);
1883  log << "Multiple scattering momentum estimate = " << p << "+-" << err_p << "\n";
1884  }
1885  }
1886  }
1887 
1888  // Get range momentum estimate.
1889 
1890  KETrack tremom_range;
1891  bool ok_range = false;
1892  if (fFitMomRange) {
1893  ok_range = fitMomentumRange(trg, tremom_range);
1894  if (ok_range) {
1895  KGTrack trg_range(trg);
1896  ok_range = updateMomentum(tremom_range, prop, trg_range);
1897  if (ok_range) {
1898  invp_range = trg_range.startTrack().getVector()(4);
1899  double var_invp = trg_range.startTrack().getError()(4, 4);
1900  double p = 0.;
1901  if (invp_range != 0.) p = 1. / invp_range;
1902  double err_p = p * p * std::sqrt(var_invp);
1903  log << "Range momentum estimate = " << p << "+-" << err_p << "\n";
1904  }
1905  }
1906  }
1907 
1908  bool result = false;
1909  if (ok_range) {
1910  tremom = tremom_range;
1911  result = ok_range;
1912  }
1913  else if (ok_ms) {
1914  tremom = tremom_ms;
1915  result = ok_ms;
1916  }
1917  return result;
1918 }
bool fFitMomRange
Fit momentum using range.
bool fFitMomMS
Fit momentum using multiple scattering.
pdgs p
Definition: selectors.fcl:22
bool fitMomentumMS(const KGTrack &trg, const Propagator &prop, KETrack &tremom) const
Estimate track momentum using multiple scattering.
bool fitMomentumRange(const KGTrack &trg, KETrack &tremom) const
Estimate track momentum using range.
bool updateMomentum(const KETrack &tremom, const Propagator &prop, KGTrack &trg) const
Set track momentum at each track surface.
bool trkf::KalmanFilterAlg::fitMomentumMS ( const KGTrack trg,
const Propagator prop,
KETrack tremom 
) const

Estimate track momentum using multiple scattering.

Estimate track momentum using multiple scattering.

Arguments:

trg - Global track. prop - Propagator. tremom - Track containing momentum estimate.

Returns: True if success.

This method estimates the momentum of the specified track using multiple scattering. As a result of calling this method, the original global track is not updated, but a KETrack is produced near the starting surface that has the estimated momentum.

The global track passed as argument should have been smoothed prior to calling this method so that all or most fits are optimal. If either the first or last fit is not optimal, return false.

This method assumes that track parameter four is 1/p. This sort of momentum estimation only makes sense if the momentum track parameter is uncorrelated with any other track parameter. The error matrix of the first and last fit is checked for this. If it is found that either error matrix has correlated track parameter four with any other track parameter, this method returns without doing anything (return false).

Definition at line 1665 of file KalmanFilterAlg.cxx.

1668 {
1669  // Get iterators pointing to the first and last tracks.
1670 
1671  const std::multimap<double, KHitTrack>& trackmap = trg.getTrackMap();
1672  if (trackmap.size() < 2) return false;
1673  std::multimap<double, KHitTrack>::const_iterator itend[2];
1674  itend[0] = trackmap.begin();
1675  itend[1] = trackmap.end();
1676  --itend[1];
1677 
1678  // Check the fit status and error matrix of the first and last
1679  // track.
1680 
1681  bool result = true;
1682 
1683  for (int i = 0; result && i < 2; ++i) {
1684  const KHitTrack& trh = itend[i]->second;
1685  KFitTrack::FitStatus stat = trh.getStat();
1686  if (stat != KFitTrack::OPTIMAL) result = false;
1687  const TrackError& err = trh.getError();
1688  for (int j = 0; j < 4; ++j) {
1689  if (err(4, j) != 0.) result = false;
1690  }
1691  }
1692  if (!result) return result;
1693 
1694  // We will periodically sample the track trajectory. At each sample
1695  // point, collect the following information.
1696  //
1697  // 1. The path distance at the sample point.
1698  // 2. The original momentum of the track at the sample point and its error.
1699  // 3. One copy of the track (KETrack) that will be propagated without noise
1700  // (infinite momentum track).
1701  // 3. A second copy of the track (KETrack) that will be propagated with the
1702  // minimum allowed momentum (range out track).
1703  // 4. A third copy of the track (KETrack) that will propagated with noise
1704  // with some intermediate momentum (noise track).
1705  //
1706  // Collect the first sample from the maximum path distance track.
1707 
1708  double s_sample = itend[1]->first;
1709  const KETrack& tre = itend[1]->second;
1710  KETrack tre_inf(tre);
1711  KTrack trk_range(tre);
1712  KETrack tre_noise(tre);
1713  tre_inf.getVector()(4) = 0.;
1714  tre_inf.getError()(4, 4) = 0.;
1715  trk_range.getVector()(4) = 100.;
1716  tre_noise.getError()(4, 4) = 0.;
1717  tre_noise.getVector()(4) = 1.;
1718  tre_noise.getError()(4, 4) = 10.;
1719  double invp0 = tre_noise.getVector()(4);
1720  double var_invp0 = tre_noise.getError()(4, 4);
1721 
1722  // Loop over fits, starting at the high path distance (low momentum)
1723  // end.
1724 
1725  for (auto it = trackmap.rbegin(); it != trackmap.rend(); ++it) {
1726  double s = it->first;
1727  const KHitTrack& trh = it->second;
1728 
1729  // Ignore non-optimal fits.
1730 
1731  KFitTrack::FitStatus stat = trh.getStat();
1732  if (stat != KFitTrack::OPTIMAL) continue;
1733 
1734  // See if this track is far enough from the previous sample to
1735  // make a new sample point.
1736 
1737  if (std::abs(s - s_sample) > fMinSampleDist) {
1738 
1739  // Propagate tracks to the current track surface.
1740 
1741  std::shared_ptr<const Surface> psurf = trh.getSurface();
1742  std::optional<double> dist_inf = prop.err_prop(tre_inf, psurf, Propagator::UNKNOWN, false);
1743  std::optional<double> dist_range =
1744  prop.vec_prop(trk_range, psurf, Propagator::UNKNOWN, false);
1745  std::optional<double> dist_noise =
1746  prop.noise_prop(tre_noise, psurf, Propagator::UNKNOWN, true);
1747 
1748  // All propagations should normally succeed. If they don't,
1749  // ignore this sample for the purpose of updating the momentum.
1750 
1751  bool momentum_updated = false;
1752  if (dist_inf && dist_range && dist_noise) {
1753 
1754  // Extract the momentum at the new sample point.
1755 
1756  double invp1 = tre_noise.getVector()(4);
1757  double var_invp1 = tre_noise.getError()(4, 4);
1758 
1759  // Get the average momentum and error for this pair of
1760  // sample points, and other data.
1761 
1762  double invp = 0.5 * (invp0 + invp1);
1763  double var_invp = 0.5 * (var_invp0 + var_invp1);
1764  double mass = tre_inf.Mass();
1765  double beta = std::sqrt(1. + mass * mass * invp * invp);
1766  double invbp = invp / beta;
1767 
1768  // Extract slope subvectors and sub-error-matrices.
1769  // We have the following variables.
1770  //
1771  // slope0 - Predicted slope vector (from infinite momentum track).
1772  // slope1 - Measured slope vector (from new sample point).
1773  // defl - Deflection (slope residual = difference between measured
1774  // and predicted slope vector).
1775  // err0 - Slope error matrix of prediction.
1776  // err1 - Slope error matrix of measurement
1777  // errc - Slope residual error matrix = err0 + err1.
1778  // errn - Noise slope error matrix divided by (1/beta*momentum).
1779 
1780  KVector<2>::type slope0 = project(tre_inf.getVector(), ublas::range(2, 4));
1781  KVector<2>::type slope1 = project(trh.getVector(), ublas::range(2, 4));
1782  KVector<2>::type defl = slope1 - slope0;
1783  KSymMatrix<2>::type err0 =
1784  project(tre_inf.getError(), ublas::range(2, 4), ublas::range(2, 4));
1785  KSymMatrix<2>::type err1 = project(trh.getError(), ublas::range(2, 4), ublas::range(2, 4));
1786  KSymMatrix<2>::type errc = err0 + err1;
1787  KSymMatrix<2>::type errn =
1788  project(tre_noise.getError(), ublas::range(2, 4), ublas::range(2, 4));
1789  errn -= err0;
1790  errn /= invbp;
1791 
1792  // Calculate updated average momentum and error.
1793 
1794  double new_invp = invp;
1795  double new_var_invp = var_invp;
1796  update_momentum(defl, errc, errn, mass, new_invp, new_var_invp);
1797 
1798  // Calculate updated momentum and error at the second sample
1799  // point.
1800 
1801  double dp = 1. / new_invp - 1. / invp;
1802  invp0 = 1. / (1. / invp1 + dp);
1803  var_invp0 = new_var_invp;
1804  momentum_updated = true;
1805 
1806  // Make sure that updated momentum is not less than minimum
1807  // allowed momentum.
1808 
1809  double invp_range = trk_range.getVector()(4);
1810  if (invp0 > invp_range) invp0 = invp_range;
1811  }
1812 
1813  // Update sample.
1814 
1815  if (momentum_updated) {
1816  s_sample = s;
1817  tre_inf = trh;
1818  tre_inf.getVector()(4) = 0.;
1819  tre_inf.getError()(4, 4) = 0.;
1820  double invp_range = trk_range.getVector()(4);
1821  trk_range = trh;
1822  trk_range.getVector()(4) = invp_range;
1823  tre_noise = trh;
1824  tre_noise.getVector()(4) = invp0;
1825  tre_noise.getError()(4, 4) = var_invp0;
1826  }
1827  }
1828  }
1829 
1830  // Propagate noise track to starting (high momentum) track surface
1831  // to get final starting momentum. This propagation should normally
1832  // always succeed, but if it doesn't, don't update the track.
1833 
1834  const KHitTrack& trh0 = itend[0]->second;
1835  std::shared_ptr<const Surface> psurf = trh0.getSurface();
1836  std::optional<double> dist_noise = prop.noise_prop(tre_noise, psurf, Propagator::UNKNOWN, true);
1837  result = dist_noise.has_value();
1838 
1839  // Update momentum-estimating track.
1840 
1841  mf::LogInfo log("KalmanFilterAlg");
1842  if (result) tremom = tre_noise;
1843 
1844  // Done.
1845 
1846  return result;
1847 }
KSymMatrix< 5 >::type TrackError
Track error matrix, dimension 5x5.
EResult err(const char *call)
FitStatus
Fit status enum.
Definition: KFitTrack.h:41
T abs(T value)
double fMinSampleDist
Minimum sample distance (for momentum measurement).
then echo File list $list not found else cat $list while read file do echo $file sed s
Definition: file_to_url.sh:60
float mass
Definition: dedx.py:47
bool trkf::KalmanFilterAlg::fitMomentumRange ( const KGTrack trg,
KETrack tremom 
) const

Estimate track momentum using range.

Estimate track momentum using range.

Arguments:

trg - Global track. prop - Propagator. tremom - Track with momentum estimate.

Returns: True if success.

This method generates a momentum-estimating track by extracting the last track from a global track, and setting its momentum to some small value.

Definition at line 1611 of file KalmanFilterAlg.cxx.

1612 {
1613  if (!trg.isValid()) return false;
1614 
1615  // Extract track with lowest momentum.
1616 
1617  const KHitTrack& trh = trg.endTrack();
1618  if (trh.getStat() == KFitTrack::INVALID)
1619  throw cet::exception("KalmanFilterAlg")
1620  << "KalmanFilterAlg::fitMomentumRange(): invalid end track\n";
1621  tremom = trh;
1622 
1623  // Set track momentum to a small value.
1624 
1625  tremom.getVector()(4) = 100.;
1626  tremom.getError()(4, 0) = 0.;
1627  tremom.getError()(4, 1) = 0.;
1628  tremom.getError()(4, 2) = 0.;
1629  tremom.getError()(4, 3) = 0.;
1630  tremom.getError()(4, 4) = 10000.;
1631 
1632  // Done.
1633 
1634  return true;
1635 }
int trkf::KalmanFilterAlg::getPlane ( ) const
inline

Preferred view plane.

Definition at line 80 of file KalmanFilterAlg.h.

bool trkf::KalmanFilterAlg::getTrace ( ) const
inline

Trace config parameters.

Definition at line 75 of file KalmanFilterAlg.h.

void trkf::KalmanFilterAlg::setPlane ( int  plane)
inline

Set preferred view plane.

Definition at line 93 of file KalmanFilterAlg.h.

void trkf::KalmanFilterAlg::setTrace ( bool  trace)
inline

Set trace config parameter.

Definition at line 88 of file KalmanFilterAlg.h.

bool trkf::KalmanFilterAlg::smoothTrack ( KGTrack trg,
KGTrack trg1,
const Propagator prop 
) const

Smooth track.

Smooth track.

Arguments:

trg - Track to be smoothed. trg1 - Track to receive result of unidirectional fit. prop - Propagator.

Returns: True if success.

The starting track should be a global track that has been fit in one direction. Fit status should be optimal at (at least) one end. It is an error if the fit status is not optimal at either end. If the fit status is optimal at both ends, do nothing, but return success.

If the second argument is non-null, save the result of the unidirectional track fit produced as a byproduct of the smoothing operation. This track can be smoothed in order to iterate the Kalman fit, etc.

The Kalman smoothing algorithm starts at the optimal end and fits the track in the reverse direction, calculating optimal track parameters at each measurement surface.

All measurements are included in the reverse fit. No incremental chisquare cut is applied.

If any measurement surface can not be reached because of a measurement error, the entire smoothing operation is considered a failure. In that case, false is returned and the track is left in an undefined state.

Definition at line 850 of file KalmanFilterAlg.cxx.

851 {
852  if (not trg.isValid()) {
853  // It is an error if the KGTrack is not valid.
854  return false;
855  }
856 
857  bool result = false;
858 
859  // Examine the track endpoints and figure out which end of the track
860  // to start from. The fit always starts at the optimal end. It is
861  // an error if neither end point is optimal. Do nothing and return
862  // success if both end points are optimal.
863 
864  const KHitTrack& trh0 = trg.startTrack();
865  const KHitTrack& trh1 = trg.endTrack();
866  KFitTrack::FitStatus stat0 = trh0.getStat();
867  KFitTrack::FitStatus stat1 = trh1.getStat();
868  bool dofit = false;
869 
870  // Remember starting direction, track, and distance.
871 
873  const KTrack* trk = 0;
874  double path = 0.;
875 
876  if (stat0 == KFitTrack::OPTIMAL) {
877  if (stat1 == KFitTrack::OPTIMAL) {
878 
879  // Both ends optimal (do nothing, return success).
880 
881  dofit = false;
882  result = true;
883  }
884  else {
885 
886  // Start optimal.
887 
888  dofit = true;
889  dir = Propagator::FORWARD;
890  trk = &trh0;
891  path = 0.;
892  }
893  }
894  else {
895  if (stat1 == KFitTrack::OPTIMAL) {
896 
897  // End optimal.
898 
899  dofit = true;
900  dir = Propagator::BACKWARD;
901  trk = &trh1;
902  path = trh1.getPath();
903  }
904  else {
905 
906  // Neither end optimal (do nothing, return failure).
907 
908  dofit = false;
909  result = false;
910  }
911  }
912  if (dofit) {
913  if (!trk)
914  throw cet::exception("KalmanFilterAlg") << "KalmanFilterAlg::smoothTrack(): no track!\n";
915 
916  // Cumulative chisquare.
917 
918  double tchisq = 0.;
919 
920  // Construct starting KFitTrack with track information and distance
921  // taken from the optimal end, but with "infinite" errors.
922 
923  TrackError err;
924  trk->getSurface()->getStartingError(err);
925  KETrack tre(*trk, err);
926  KFitTrack trf(tre, path, tchisq);
927 
928  // Make initial reference track to be same as initial fit track.
929 
930  KTrack ref(trf);
931 
932  // Loop over KHitTracks contained in KGTrack.
933 
934  std::multimap<double, KHitTrack>::iterator it;
935  std::multimap<double, KHitTrack>::iterator itend;
936  switch (dir) {
937  case Propagator::FORWARD:
938  it = trg.getTrackMap().begin();
939  itend = trg.getTrackMap().end();
940  break;
942  it = trg.getTrackMap().end();
943  itend = trg.getTrackMap().begin();
944  break;
945  default:
946  throw cet::exception("KalmanFilterAlg")
947  << "KalmanFilterAlg::smoothTrack(): invalid direction\n";
948  } // switch
949 
950  mf::LogInfo log("KalmanFilterAlg");
951 
952  // Loop starts here.
953 
954  result = true; // Result success unless we find an error.
955  int step = 0; // Step count.
956  while (dofit && it != itend) {
957  ++step;
958  if (fTrace) {
959  log << "Smooth Step " << step << "\n";
960  log << "Reverse fit track:\n";
961  log << trf;
962  }
963 
964  // For backward fit, decrement iterator at start of loop.
965 
966  if (dir == Propagator::BACKWARD) --it;
967 
968  KHitTrack& trh = (*it).second;
969  if (fTrace) {
970  log << "Forward track:\n";
971  log << trh;
972  }
973 
974  // Extract measurement.
975 
976  const KHitBase& hit = *(trh.getHit());
977 
978  // Propagate KFitTrack to the next track surface.
979 
980  std::shared_ptr<const Surface> psurf = trh.getSurface();
981  std::optional<double> dist = prop.noise_prop(trf, psurf, Propagator::UNKNOWN, true, &ref);
982 
983  // Check if propagation succeeded. If propagation fails, this
984  // measurement will be dropped from the unidirectional fit
985  // track. This measurement will still be in the original
986  // track, but with a status other than optimal.
987 
988  if (dist) {
989 
990  // Propagation succeeded.
991  // Update cumulative path distance and track status.
992 
993  double ds = *dist;
994  path += ds;
995  trf.setPath(path);
996  if (dir == Propagator::FORWARD)
997  trf.setStat(KFitTrack::FORWARD_PREDICTED);
998  else {
999  trf.setStat(KFitTrack::BACKWARD_PREDICTED);
1000  }
1001  if (fTrace) {
1002  log << "Reverse fit track after propagation:\n";
1003  log << " Propagation distance = " << ds << "\n";
1004  log << trf;
1005  }
1006 
1007  // See if we have the proper information to calculate an optimal track
1008  // at this surface (should normally be possible).
1009 
1010  KFitTrack::FitStatus stat = trh.getStat();
1011  KFitTrack::FitStatus newstat = trf.getStat();
1012 
1013  if ((newstat == KFitTrack::FORWARD_PREDICTED && stat == KFitTrack::BACKWARD) ||
1014  (newstat == KFitTrack::BACKWARD_PREDICTED && stat == KFitTrack::FORWARD)) {
1015 
1016  // Update stored KHitTrack to be optimal.
1017 
1018  bool ok = trh.combineFit(trf);
1019 
1020  // Update the stored path distance to be from the currently fitting track.
1021 
1022  trh.setPath(trf.getPath());
1023 
1024  // Update reference track.
1025 
1026  ref = trh;
1027 
1028  // If combination failed, abandon the fit and return failure.
1029 
1030  if (!ok) {
1031  dofit = false;
1032  result = false;
1033  break;
1034  }
1035  if (fTrace) {
1036  log << "Combined track:\n";
1037  log << trh;
1038  }
1039  }
1040 
1041  // Update measurement predction using current track hypothesis.
1042 
1043  if (not hit.predict(trf, prop, &ref)) {
1044 
1045  // Abandon the fit and return failure.
1046 
1047  dofit = false;
1048  result = false;
1049  break;
1050  }
1051 
1052  // Prediction succeeded. Get incremental chisquare. If this
1053  // hit fails incremental chisquare cut, this hit will be
1054  // dropped from the unidirecitonal Kalman fit track, but may
1055  // still be in the smoothed track.
1056 
1057  double chisq = hit.getChisq();
1058  if (chisq < fMaxSmoothIncChisq) {
1059 
1060  // Update the reverse fitting track using the current measurement
1061  // (both track parameters and status).
1062 
1063  KFitTrack trf0(trf);
1064  hit.update(trf);
1065  bool update_ok = trf.isValid();
1066  if (!update_ok)
1067  trf = trf0;
1068  else {
1069  tchisq += chisq;
1070  trf.setChisq(tchisq);
1071 
1072  if (dir == Propagator::FORWARD)
1073  trf.setStat(KFitTrack::FORWARD);
1074  else {
1075  trf.setStat(KFitTrack::BACKWARD);
1076  }
1077  if (fTrace) {
1078  log << "Reverse fit track after update:\n";
1079  log << trf;
1080  }
1081  if (fGTrace && fCanvases.size() > 0) {
1082  auto marker_it = fMarkerMap.find(hit.getID());
1083  if (marker_it != fMarkerMap.end()) {
1084  TMarker* marker = marker_it->second;
1085  marker->SetMarkerColor(kOrange);
1086  }
1087  fCanvases.back()->Update();
1088  }
1089 
1090  // If unidirectional track pointer is not null, make a
1091  // KHitTrack and save it in the unidirectional track.
1092 
1093  if (trg1 != 0) {
1094  KHitTrack trh1(trf, trh.getHit());
1095  trg1->addTrack(trh1);
1096  }
1097  }
1098  }
1099  }
1100 
1101  // For forward fit, increment iterator at end of loop.
1102 
1103  if (dir == Propagator::FORWARD) ++it;
1104 
1105  } // Loop over KHitTracks.
1106 
1107  // If fit was successful and the unidirectional track pointer
1108  // is not null and the track is valid, set the fit status of
1109  // the last added KHitTrack to optimal.
1110 
1111  if (result && trg1 != 0 && trg1->isValid()) {
1112  if (dir == Propagator::FORWARD)
1113  trg1->endTrack().setStat(KFitTrack::OPTIMAL);
1114  else {
1115  trg1->startTrack().setStat(KFitTrack::OPTIMAL);
1116  }
1117  }
1118 
1119  // Recalibrate track map.
1120 
1121  trg.recalibrate();
1122 
1123  } // Do fit.
1124 
1125  // Get the final chisquare.
1126 
1127  double fchisq = 0.5 * (trg.startTrack().getChisq() + trg.endTrack().getChisq());
1128 
1129  // Summary.
1130 
1131  mf::LogInfo log("KalmanFilterAlg");
1132  log << "KalmanFilterAlg smooth track summary.\n"
1133  << "Smooth direction = " << (dir == Propagator::FORWARD ? "FORWARD" : "BACKWARD") << "\n"
1134  << "Track has " << trg.numHits() << " hits.\n"
1135  << "Track length = " << trg.endTrack().getPath() - trg.startTrack().getPath() << "\n"
1136  << "Track chisquare = " << fchisq << "\n";
1137 
1138  return result;
1139 }
bool fGTrace
Graphical trace flag.
std::map< int, TMarker * > fMarkerMap
Markers in current canvas.
KSymMatrix< 5 >::type TrackError
Track error matrix, dimension 5x5.
bool fTrace
Trace flag.
EResult err(const char *call)
FitStatus
Fit status enum.
Definition: KFitTrack.h:41
process_name hit
Definition: cheaterreco.fcl:51
BEGIN_PROLOG triggeremu_data_config_icarus settings PMTADCthresholds sequence::icarus_stage0_multiTPC_TPC physics sequence::icarus_stage0_EastHits_TPC physics sequence::icarus_stage0_WestHits_TPC physics producers purityana0 caloskimCalorimetryCryoE physics caloskimCalorimetryCryoW physics path
std::vector< std::unique_ptr< TCanvas > > fCanvases
Graphical trace canvases.
double fMaxSmoothIncChisq
Maximum incremental chisquare to accept a hit in smooth phase.
tuple dir
Definition: dropbox.py:28
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
PropDirection
Propagation direction enum.
Definition: Propagator.h:94
bool trkf::KalmanFilterAlg::smoothTrackIter ( int  nsmooth,
KGTrack trg,
const Propagator prop 
) const

Iteratively smooth a track.

Iteratively smooth a track.

Arguments:

nsmooth - Number of iterations. trg - Track to be smoothed. prop - Propagator.

Returns: True if success.

The initial track should have been unidirectionally fit.

Definition at line 2078 of file KalmanFilterAlg.cxx.

2079 {
2080  bool ok = true;
2081 
2082  // Call smoothTrack method in a loop.
2083 
2084  for (int ismooth = 0; ok && ismooth < nsmooth - 1; ++ismooth) {
2085  KGTrack trg1{trg.getPrefPlane()};
2086  ok = smoothTrack(trg, &trg1, prop);
2087  if (ok) trg = trg1;
2088  }
2089 
2090  // Do one final smooth without generating a new unidirectional
2091  // track.
2092 
2093  if (ok) ok = smoothTrack(trg, 0, prop);
2094 
2095  // Done.
2096 
2097  return ok;
2098 }
bool smoothTrack(KGTrack &trg, KGTrack *trg1, const Propagator &prop) const
Smooth track.
bool trkf::KalmanFilterAlg::updateMomentum ( const KETrack tremom,
const Propagator prop,
KGTrack trg 
) const

Set track momentum at each track surface.

Set track momentum at each track surface.

Arguments:

tremom - Track containing momentum estimate. prop - Propagator. trg - Global track to be updated.

The track containing the momentum estimate is propagated to the first or last track fit of the global track (whichever is closer), then the momentum estimate is transfered to that track fit. In similar fashion, the momentum estimate is successively transfered from that track fit to each track fit of the global track.

Only momentum track parameters of the global track fits are updated. Other track parameters and their errors are unmodified. Unreachable track fits are deleted from the global track. Overall failure will occur if the momentum-estimating track can't be propagated to the initial track fit, or if the final global track has no valid track fits.

Definition at line 1942 of file KalmanFilterAlg.cxx.

1945 {
1946  // Get modifiable track map.
1947 
1948  std::multimap<double, KHitTrack>& trackmap = trg.getTrackMap();
1949 
1950  // If track map is empty, immediately return failure.
1951 
1952  if (trackmap.size() == 0) return false;
1953 
1954  // Make trial propagations to the first and last track fit to figure
1955  // out which track fit is closer to the momentum estimating track.
1956 
1957  // Find distance to first track fit.
1958 
1959  KETrack tre0(tremom);
1960  std::optional<double> dist0 =
1961  prop.vec_prop(tre0, trackmap.begin()->second.getSurface(), Propagator::UNKNOWN, false, 0, 0);
1962  // Find distance to last track fit.
1963 
1964  KETrack tre1(tremom);
1965  std::optional<double> dist1 =
1966  prop.vec_prop(tre1, trackmap.rbegin()->second.getSurface(), Propagator::UNKNOWN, false, 0, 0);
1967 
1968  // Based on distances, make starting iterator and direction flag.
1969 
1970  bool forward = true;
1971  std::multimap<double, KHitTrack>::iterator it = trackmap.begin();
1972  if (dist0) {
1973 
1974  // Propagation to first track succeeded.
1975 
1976  if (dist1) {
1977 
1978  // Propagation to both ends succeeded. If the last track is
1979  // closer, initialize iterator and direction flag for reverse
1980  // direction.
1981 
1982  if (std::abs(*dist0) > std::abs(*dist1)) {
1983  it = trackmap.end();
1984  --it;
1985  forward = false;
1986  }
1987  }
1988  }
1989  else {
1990 
1991  // Propagation to first track failed. Initialize iterator and
1992  // direction flag for reverse direction, provided that the
1993  // propagation to the last track succeeded.
1994 
1995  if (dist1) {
1996  it = trackmap.end();
1997  --it;
1998  forward = false;
1999  }
2000  else {
2001 
2002  // Propagation to both ends failed. Return failure.
2003 
2004  return false;
2005  }
2006  }
2007 
2008  // Loop over track fits in global track.
2009 
2010  KETrack tre(tremom);
2011  bool done = false;
2012  while (!done) {
2013  KHitTrack& trh = it->second;
2014  if (trh.getStat() == KFitTrack::INVALID)
2015  throw cet::exception("KalmanFilterAlg")
2016  << "KalmanFilterAlg::updateMomentum(): invalid track\n";
2017 
2018  // Propagate momentum-estimating track to current track surface
2019  // and update momentum.
2020 
2021  std::optional<double> dist = prop.noise_prop(tre, trh.getSurface(), Propagator::UNKNOWN, true);
2022 
2023  // Copy momentum to global track.
2024 
2025  std::multimap<double, KHitTrack>::iterator erase_it = trackmap.end();
2026  if (dist) {
2027  trh.getVector()(4) = tre.getVector()(4);
2028  trh.getError()(4, 0) = 0.;
2029  trh.getError()(4, 1) = 0.;
2030  trh.getError()(4, 2) = 0.;
2031  trh.getError()(4, 3) = 0.;
2032  trh.getError()(4, 4) = tre.getError()(4, 4);
2033  }
2034  else {
2035 
2036  // If the propagation failed, remember that we are supposed to
2037  // erase this track from the global track.
2038 
2039  erase_it = it;
2040  }
2041 
2042  // Advance the iterator and set the done flag.
2043 
2044  if (forward) {
2045  ++it;
2046  done = (it == trackmap.end());
2047  }
2048  else {
2049  done = (it == trackmap.begin());
2050  if (!done) --it;
2051  }
2052 
2053  // Update momentum-estimating track from just-updated global track
2054  // fit, or erase global track fit.
2055 
2056  if (erase_it == trackmap.end())
2057  tre = trh;
2058  else
2059  trackmap.erase(erase_it);
2060  }
2061 
2062  return not empty(trackmap);
2063 }
T abs(T value)
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
bool empty(FixedBins< T, C > const &) noexcept
Definition: FixedBins.h:555

Member Data Documentation

std::vector<std::unique_ptr<TCanvas> > trkf::KalmanFilterAlg::fCanvases
mutableprivate

Graphical trace canvases.

Definition at line 179 of file KalmanFilterAlg.h.

bool trkf::KalmanFilterAlg::fFitMomMS
private

Fit momentum using multiple scattering.

Definition at line 167 of file KalmanFilterAlg.h.

bool trkf::KalmanFilterAlg::fFitMomRange
private

Fit momentum using range.

Definition at line 166 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fGapDist
private

Minimum gap distance.

Definition at line 163 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fGoodPErr
private

Pointing error threshold for switching to free propagation.

Definition at line 150 of file KalmanFilterAlg.h.

bool trkf::KalmanFilterAlg::fGTrace
private

Graphical trace flag.

Definition at line 168 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fGTraceWH
private

Window height.

Definition at line 170 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fGTraceWW
private

Window width.

Definition at line 169 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fGTraceXMax
private

Graphical trace maximum x.

Definition at line 172 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fGTraceXMin
private

Graphical trace minimum x.

Definition at line 171 of file KalmanFilterAlg.h.

std::vector<double> trkf::KalmanFilterAlg::fGTraceZMax
private

Graphical trace maximum z for each view.

Definition at line 174 of file KalmanFilterAlg.h.

std::vector<double> trkf::KalmanFilterAlg::fGTraceZMin
private

Graphical trace minimum z for each view.

Definition at line 173 of file KalmanFilterAlg.h.

TVirtualPad* trkf::KalmanFilterAlg::fInfoPad
mutableprivate

Information pad.

Definition at line 181 of file KalmanFilterAlg.h.

std::map<int, TMarker*> trkf::KalmanFilterAlg::fMarkerMap
mutableprivate

Markers in current canvas.

Definition at line 183 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fMaxEndChisq
private

Maximum incremental chisquare for endpoint hit.

Definition at line 154 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fMaxIncChisq
private

Maximum incremental chisquare to accept a hit.

Definition at line 151 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fMaxLDist
private

Maximum distance for linearized propagation.

Definition at line 156 of file KalmanFilterAlg.h.

int trkf::KalmanFilterAlg::fMaxNoiseHits
private

Maximum number of hits in noise cluster.

Definition at line 164 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fMaxPErr
private

Maximum pointing error for free propagation.

Definition at line 149 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fMaxPredDist
private

Maximum prediciton distance to accept a hit.

Definition at line 157 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fMaxPropDist
private

Maximum propagation distance to candidate surface.

Definition at line 159 of file KalmanFilterAlg.h.

int trkf::KalmanFilterAlg::fMaxSamePlane
private

Maximum consecutive hits in same plane.

Definition at line 162 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fMaxSeedIncChisq
private

Maximum incremental chisquare to accept a hit in seed phase.

Definition at line 152 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fMaxSeedPredDist
private

Maximum prediciton distance to accept a hit in seed phase.

Definition at line 158 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fMaxSmoothIncChisq
private

Maximum incremental chisquare to accept a hit in smooth phase.

Definition at line 153 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fMaxSortDist
private

Sort high distance threshold.

Definition at line 161 of file KalmanFilterAlg.h.

TPaveText* trkf::KalmanFilterAlg::fMessages
mutableprivate

Message box.

Definition at line 182 of file KalmanFilterAlg.h.

int trkf::KalmanFilterAlg::fMinLHits
private

Minimum number of hits to turn off linearized propagation.

Definition at line 155 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fMinSampleDist
private

Minimum sample distance (for momentum measurement).

Definition at line 165 of file KalmanFilterAlg.h.

double trkf::KalmanFilterAlg::fMinSortDist
private

Sort low distance threshold.

Definition at line 160 of file KalmanFilterAlg.h.

std::vector<TVirtualPad*> trkf::KalmanFilterAlg::fPads
mutableprivate

View pads in current canvas.

Definition at line 180 of file KalmanFilterAlg.h.

int trkf::KalmanFilterAlg::fPlane
private

Preferred view plane.

Definition at line 178 of file KalmanFilterAlg.h.

bool trkf::KalmanFilterAlg::fTrace
private

Trace flag.

Definition at line 148 of file KalmanFilterAlg.h.


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