LArSoft  v09_90_00
Liquid Argon Software toolkit - https://larsoft.org/
KalmanFilterAlg.cxx
Go to the documentation of this file.
1 
13 #include "boost/numeric/ublas/matrix_proxy.hpp"
14 #include "boost/numeric/ublas/vector_proxy.hpp"
15 #include "cetlib_except/exception.h"
24 
25 #include "Rtypes.h"
26 #include "TCanvas.h"
27 #include "TGaxis.h"
28 #include "TLegend.h"
29 #include "TLegendEntry.h"
30 #include "TMarker.h"
31 #include "TObject.h"
32 #include "TPad.h"
33 #include "TPaveText.h"
34 #include "TText.h"
35 #include "TVirtualPad.h"
36 
37 // Local functions.
38 
39 namespace {
40 
41  void hit_position(const trkf::KHitBase& hit, double& z, double& x)
42  {
43 
44  // Map hit -> (z,x) for plotting.
45 
46  // Default result.
47 
48  z = 0.;
49  x = 0.;
50 
51  // Cast this hit to KHit<1>. Only know how to handle 1D hits.
52 
53  const trkf::KHit<1>* phit1 = dynamic_cast<const trkf::KHit<1>*>(&hit);
54  if (phit1) {
55  const std::shared_ptr<const trkf::Surface>& psurf = hit.getMeasSurface();
56 
57  // Handle SurfYZPlane.
58 
59  if (const trkf::SurfYZPlane* pyz = dynamic_cast<const trkf::SurfYZPlane*>(&*psurf)) {
60 
61  // Now finished doing casts.
62  // We have a kind of hit and measurement surface that we know how to handle.
63 
64  // Get x coordinate from hit and surface.
65 
66  x = pyz->x0() + phit1->getMeasVector()(0);
67 
68  // Get z position from surface parameters.
69  // The "z" position is actually calculated as the perpendicular distance
70  // from a corner, which is a proxy for wire number.
71 
72  double z0 = pyz->z0();
73  double y0 = pyz->y0();
74  double phi = pyz->phi();
76  double ymax = geom->DetHalfWidth();
77  if (phi > 0.)
78  z = z0 * std::cos(phi) + (ymax - y0) * std::sin(phi);
79  else
80  z = z0 * std::cos(phi) - (ymax + y0) * std::sin(phi);
81  }
82  else if (const trkf::SurfYZLine* pyz = dynamic_cast<const trkf::SurfYZLine*>(&*psurf)) {
83 
84  // Now finished doing casts.
85  // We have a kind of hit and measurement surface that we know how to handle.
86 
87  // Get x coordinate from surface.
88 
89  x = pyz->x0();
90 
91  // Get z position from surface parameters.
92  // The "z" position is actually calculated as the perpendicular distance
93  // from a corner, which is a proxy for wire number.
94 
95  double z0 = pyz->z0();
96  double y0 = pyz->y0();
97  double phi = pyz->phi();
99  double ymax = geom->DetHalfWidth();
100  if (phi > 0.)
101  z = z0 * std::cos(phi) + (ymax - y0) * std::sin(phi);
102  else
103  z = z0 * std::cos(phi) - (ymax + y0) * std::sin(phi);
104  }
105  }
106  }
107 
108  void update_momentum(const trkf::KVector<2>::type& defl,
109  const trkf::KSymMatrix<2>::type& errc,
110  const trkf::KSymMatrix<2>::type& errn,
111  double mass,
112  double& invp,
113  double& var_invp)
114  // Momentum updater.
115  //
116  // Arguments: defl - Deflection (2D slope residual between two
117  // sampled points on track).
118  // errc - Error matrix of deflection residual
119  // exclusive of multiple scattering.
120  // errn - Part of deflection noise error matrix
121  // proportional to 1/(beta*momentum).
122  // mass - Mass of particle.
123  // invp - Modified 1/momentum track parameter.
124  // var_invp - Modified 1/momentum variance.
125  //
126  // Returns: True if success.
127  //
128  // The total deflection residual error matrix is
129  //
130  // err = errc + [1/(beta*momentum)] * errn.
131  //
132  // The inverse momentum and error are updated using using log-likelihood
133  //
134  // log(L) = -0.5 * log(det(err)) - 0.5 * defl^T * err^(-1) * defl.
135  // log(L) = -0.5 * tr(log(err)) - 0.5 * defl^T * err^(-1) * defl.
136  //
137  {
138  // Calculate original k = 1./(beta*p), and original variance of k.
139 
140  double invp2 = invp * invp;
141  double invp3 = invp * invp2;
142  double invp4 = invp2 * invp2;
143  double mass2 = mass * mass;
144  double k = std::sqrt(invp2 + mass2 * invp4);
145  double dkdinvp = (invp + 2. * mass2 * invp3) / k;
146  double vark = var_invp * dkdinvp * dkdinvp;
147 
148  // First, find current inverse error matrix using momentum hypothesis.
149 
150  trkf::KSymMatrix<2>::type inverr = errc + k * errn;
151  trkf::syminvert(inverr);
152 
153  // Find the first and second derivatives of the log likelihood
154  // with respact to k.
155 
156  trkf::KMatrix<2, 2>::type temp1 = prod(inverr, errn);
157  trkf::KMatrix<2, 2>::type temp2 = prod(temp1, temp1);
158 
159  // KJK-FIXME: Ideally, the explicit type specification below
160  // ('trk::KVector<T>::type') should be replaced with the 'auto'
161  // keyword. However, Boost 1.66.0a runs into a static-assert
162  // failure with the 'complexity' of vtemp2 not being 0. I am not
163  // sure what the complexity means.
164  //
165  // If 'trkf::KVector<2>::type vtemp2 = prod(temp1, vtemp1)' is
166  // used below, then the GCC 7.3 compiler thinks there may be an
167  // uninitialized variable in vtemp2. A workaround is to use the
168  // Boost interface that fills a default-constructed vector. Note
169  // that this workaround may hide an actual problem. It is the
170  // responsibility of the maintainer of this code to determine if
171  // an actual problem exists.
172 
173  trkf::KVector<2>::type vtemp1 = prod(inverr, defl);
174  trkf::KVector<2>::type vtemp2 = prod(temp1, vtemp1);
175  trkf::KVector<2>::type vtemp3 = prod(temp1, vtemp2);
176  double derivk1 = -0.5 * trkf::trace(temp1) + 0.5 * inner_prod(defl, vtemp2);
177  double derivk2 = 0.5 * trkf::trace(temp2) - inner_prod(defl, vtemp3);
178 
179  // We expect the log-likelihood to be most nearly Gaussian
180  // with respect to variable q = k^(-1/2) = std::sqrt(beta*p).
181  // Therefore, transform the original variables and log-likelihood
182  // derivatives to q-space.
183 
184  double q = 1. / std::sqrt(k);
185  double varq = vark / (4. * k * k * k);
186  double derivq1 = (-2. * k / q) * derivk1;
187  double derivq2 = 6. * k * k * derivk1 + 4. * k * k * k * derivk2;
188 
189  if (derivq2 < 0.) {
190 
191  // Estimate the measurement in q-space.
192 
193  double q1 = q - derivq1 / derivq2;
194  double varq1 = -1. / derivq2;
195 
196  // Get updated estimated q, combining original estimate and
197  // update from the current measurement.
198 
199  double newvarq = 1. / (1. / varq + 1. / varq1);
200  double newq = newvarq * (q / varq + q1 / varq1);
201  q = newq;
202  varq = newvarq;
203 
204  // Calculate updated c = 1./p and variance.
205 
206  double q2 = q * q;
207  double q4 = q2 * q2;
208  double c2 = 2. / (q4 + q2 * std::sqrt(q4 + 4. * mass2));
209  double c = std::sqrt(c2);
210  double dcdq = -2. * (c / q) * (1. + mass2 * c2) / (1. + 2. * mass2 * c2);
211  double varc = varq * dcdq * dcdq;
212 
213  // Update result.
214 
215  invp = c;
216  var_invp = varc;
217  }
218  }
219 }
220 
222 
224  : fTrace(false)
225  , fMaxPErr(0.)
226  , fGoodPErr(0.)
227  , fMaxIncChisq(0.)
228  , fMaxSeedIncChisq(0.)
229  , fMaxSmoothIncChisq(0.)
230  , fMaxEndChisq(0.)
231  , fMinLHits(0)
232  , fMaxLDist(0.)
233  , fMaxPredDist(0.)
234  , fMaxSeedPredDist(0.)
235  , fMaxPropDist(0.)
236  , fMinSortDist(0.)
237  , fMaxSortDist(0.)
238  , fMaxSamePlane(0)
239  , fGapDist(0.)
240  , fMaxNoiseHits(0)
241  , fMinSampleDist(0.)
242  , fFitMomRange(false)
243  , fFitMomMS(false)
244  , fGTrace(false)
245  , fGTraceWW(0)
246  , fGTraceWH(0)
247  , fPlane(-1)
248  , fInfoPad(0)
249  , fMessages(0)
250 {
251  mf::LogInfo("KalmanFilterAlg") << "KalmanFilterAlg instantiated.";
252 
253  fTrace = pset.get<bool>("Trace");
254  fMaxPErr = pset.get<double>("MaxPErr");
255  fGoodPErr = pset.get<double>("GoodPErr");
256  fMaxIncChisq = pset.get<double>("MaxIncChisq");
257  fMaxSeedIncChisq = pset.get<double>("MaxSeedIncChisq");
258  fMaxSmoothIncChisq = pset.get<double>("MaxSmoothIncChisq");
259  fMaxEndChisq = pset.get<double>("MaxEndChisq");
260  fMinLHits = pset.get<int>("MinLHits");
261  fMaxLDist = pset.get<double>("MaxLDist");
262  fMaxPredDist = pset.get<double>("MaxPredDist");
263  fMaxSeedPredDist = pset.get<double>("MaxSeedPredDist");
264  fMaxPropDist = pset.get<double>("MaxPropDist");
265  fMinSortDist = pset.get<double>("MinSortDist");
266  fMaxSortDist = pset.get<double>("MaxSortDist");
267  fMaxSamePlane = pset.get<int>("MaxSamePlane");
268  fGapDist = pset.get<double>("GapDist");
269  fMaxNoiseHits = pset.get<double>("MaxNoiseHits");
270  fMinSampleDist = pset.get<double>("MinSampleDist");
271  fFitMomRange = pset.get<bool>("FitMomRange");
272  fFitMomMS = pset.get<bool>("FitMomMS");
273  fGTrace = pset.get<bool>("GTrace");
274  if (fGTrace) {
275  fGTraceWW = pset.get<double>("GTraceWW");
276  fGTraceWH = pset.get<double>("GTraceWH");
277  fGTraceXMin = pset.get<double>("GTraceXMin");
278  fGTraceXMax = pset.get<double>("GTraceXMax");
279  fGTraceZMin = pset.get<std::vector<double>>("GTraceZMin");
280  fGTraceZMax = pset.get<std::vector<double>>("GTraceZMax");
281  }
282 }
283 
307  KGTrack& trg,
308  const Propagator& prop,
311  bool linear) const
312 {
313  // Direction must be forward or backward (unknown is not allowed).
314 
315  if (dir != Propagator::FORWARD && dir != Propagator::BACKWARD)
316  throw cet::exception("KalmanFilterAlg") << "No direction for Kalman fit.\n";
317 
318  // Set up canvas for graphical trace.
319 
320  if (fGTrace) {
321 
322  // Make a new canvas with a unique name.
323 
324  static int cnum = 0;
325  ++cnum;
326  std::ostringstream ostr;
327  ostr << "khit" << cnum;
328  fCanvases.emplace_back(
329  new TCanvas(ostr.str().c_str(), ostr.str().c_str(), fGTraceWW, fGTraceWH));
330  fPads.clear();
331  fMarkerMap.clear();
332  int nview = 3;
333  fCanvases.back()->Divide(2, nview / 2 + 1);
334 
335  // Make subpad for each view.
336 
337  for (int iview = 0; iview < nview; ++iview) {
338 
339  std::ostringstream ostr;
340  ostr << "Plane " << iview;
341 
342  fCanvases.back()->cd(iview + 1);
343  double zmin = 0.06;
344  double zmax = 0.94;
345  double xmin = 0.04;
346  double xmax = 0.95;
347  TPad* p = new TPad(ostr.str().c_str(), ostr.str().c_str(), zmin, xmin, zmax, xmax);
348  p->SetBit(kCanDelete); // Give away ownership.
349  p->Range(fGTraceZMin[iview], fGTraceXMin, fGTraceZMax[iview], fGTraceXMax);
350  p->SetFillStyle(4000); // Transparent.
351  p->Draw();
352  fPads.push_back(p);
353 
354  // Draw label.
355 
356  TText* t = new TText(zmax - 0.02, xmax - 0.03, ostr.str().c_str());
357  t->SetBit(kCanDelete); // Give away ownership.
358  t->SetTextAlign(33);
359  t->SetTextFont(42);
360  t->SetTextSize(0.04);
361  t->Draw();
362 
363  // Draw axes.
364 
365  TGaxis* pz1 =
366  new TGaxis(zmin, xmin, zmax, xmin, fGTraceZMin[iview], fGTraceZMax[iview], 510, "");
367  pz1->SetBit(kCanDelete); // Give away ownership.
368  pz1->Draw();
369 
370  TGaxis* px1 = new TGaxis(zmin, xmin, zmin, xmax, fGTraceXMin, fGTraceXMax, 510, "");
371  px1->SetBit(kCanDelete); // Give away ownership.
372  px1->Draw();
373 
374  TGaxis* pz2 =
375  new TGaxis(zmin, xmax, zmax, xmax, fGTraceZMin[iview], fGTraceZMax[iview], 510, "-");
376  pz2->SetBit(kCanDelete); // Give away ownership.
377  pz2->Draw();
378 
379  TGaxis* px2 = new TGaxis(zmax, xmin, zmax, xmax, fGTraceXMin, fGTraceXMax, 510, "L+");
380  px2->SetBit(kCanDelete); // Give away ownership.
381  px2->Draw();
382  }
383 
384  // Draw legend.
385 
386  fCanvases.back()->cd(nview + 1);
387  fInfoPad = gPad;
388  TLegend* leg = new TLegend(0.6, 0.5, 0.99, 0.99);
389  leg->SetBit(kCanDelete); // Give away ownership.
390 
391  TLegendEntry* entry = 0;
392  TMarker* m = 0;
393 
394  m = new TMarker(0., 0., 20);
395  m->SetBit(kCanDelete);
396  m->SetMarkerSize(1.2);
397  m->SetMarkerColor(kRed);
398  entry = leg->AddEntry(m, "Hits on Track", "P");
399  entry->SetBit(kCanDelete);
400 
401  m = new TMarker(0., 0., 20);
402  m->SetBit(kCanDelete);
403  m->SetMarkerSize(1.2);
404  m->SetMarkerColor(kOrange);
405  entry = leg->AddEntry(m, "Smoothed Hits on Track", "P");
406  entry->SetBit(kCanDelete);
407 
408  m = new TMarker(0., 0., 20);
409  m->SetBit(kCanDelete);
410  m->SetMarkerSize(1.2);
411  m->SetMarkerColor(kBlack);
412  entry = leg->AddEntry(m, "Available Hits", "P");
413  entry->SetBit(kCanDelete);
414 
415  m = new TMarker(0., 0., 20);
416  m->SetBit(kCanDelete);
417  m->SetMarkerSize(1.2);
418  m->SetMarkerColor(kBlue);
419  entry = leg->AddEntry(m, "Rejected Hits", "P");
420  entry->SetBit(kCanDelete);
421 
422  m = new TMarker(0., 0., 20);
423  m->SetBit(kCanDelete);
424  m->SetMarkerSize(1.2);
425  m->SetMarkerColor(kGreen);
426  entry = leg->AddEntry(m, "Removed Hits", "P");
427  entry->SetBit(kCanDelete);
428 
429  m = new TMarker(0., 0., 20);
430  m->SetBit(kCanDelete);
431  m->SetMarkerSize(1.2);
432  m->SetMarkerColor(kMagenta);
433  entry = leg->AddEntry(m, "Seed Hits", "P");
434  entry->SetBit(kCanDelete);
435 
436  m = new TMarker(0., 0., 20);
437  m->SetBit(kCanDelete);
438  m->SetMarkerSize(1.2);
439  m->SetMarkerColor(kCyan);
440  entry = leg->AddEntry(m, "Unreachable Seed Hits", "P");
441  entry->SetBit(kCanDelete);
442 
443  leg->Draw();
444 
445  // Draw message box.
446 
447  fMessages = new TPaveText(0.01, 0.01, 0.55, 0.99, "");
448  fMessages->SetBit(kCanDelete);
449  fMessages->SetTextFont(42);
450  fMessages->SetTextSize(0.04);
451  fMessages->SetTextAlign(12);
452  fMessages->Draw();
453  TText* t = fMessages->AddText("Enter buildTrack");
454  t->SetBit(kCanDelete);
455 
456  fCanvases.back()->Update();
457  }
458 
459  // Sort container using this seed track.
460 
461  hits.sort(trk, true, prop, Propagator::UNKNOWN);
462 
463  // Draw hits and populate hit->marker map.
464 
465  if (fGTrace && fCanvases.size() > 0) {
466 
467  // Loop over sorted KHitGroups.
468  // Paint sorted seed hits magenta.
469 
470  const std::list<KHitGroup>& groups = hits.getSorted();
471  for (auto const& gr : groups) {
472 
473  // Loop over hits in this group.
474 
475  const std::vector<std::shared_ptr<const KHitBase>>& phits = gr.getHits();
476  for (auto const& phit : phits) {
477  const KHitBase& hit = *phit;
478  int pl = hit.getMeasPlane();
479  if (pl >= 0 && pl < int(fPads.size())) {
480  double z = 0.;
481  double x = 0.;
482  hit_position(hit, z, x);
483  TMarker* marker = new TMarker(z, x, 20);
484  fMarkerMap[hit.getID()] = marker;
485  fPads[pl]->cd();
486  marker->SetBit(kCanDelete); // Give away ownership.
487  marker->SetMarkerSize(0.5);
488  marker->SetMarkerColor(kMagenta);
489  marker->Draw();
490  }
491  }
492  }
493 
494  // Loop over unsorted KHitGroups.
495  // Paint unsorted seed hits cyan.
496  // There should be few, if any, unsorted seed hits.
497 
498  const std::list<KHitGroup>& ugroups = hits.getUnsorted();
499  for (auto const& gr : ugroups) {
500 
501  // Loop over hits in this group.
502 
503  const std::vector<std::shared_ptr<const KHitBase>>& phits = gr.getHits();
504  for (auto const& phit : phits) {
505  const KHitBase& hit = *phit;
506  int pl = hit.getMeasPlane();
507  if (pl >= 0 && pl < int(fPads.size())) {
508  double z = 0.;
509  double x = 0.;
510  hit_position(hit, z, x);
511  TMarker* marker = new TMarker(z, x, 20);
512  fMarkerMap[hit.getID()] = marker;
513  fPads[pl]->cd();
514  marker->SetBit(kCanDelete); // Give away ownership.
515  marker->SetMarkerSize(0.5);
516  marker->SetMarkerColor(kCyan);
517  marker->Draw();
518  }
519  }
520  }
521  fCanvases.back()->Update();
522  }
523 
524  // Loop over measurements (KHitGroup) from sorted list.
525 
526  double tchisq = 0.; // Cumulative chisquare.
527  double path = 0.; // Cumulative path distance.
528  int step = 0; // Step count.
529  int nsame = 0; // Number of consecutive measurements in same plane.
530  int last_plane = -1; // Last plane.
531  bool has_pref_plane = false; // Set to true when first preferred plane hit is added to track.
532 
533  // Make a copy of the starting track, in the form of a KFitTrack,
534  // which we will update as we go.
535 
536  TrackError err;
537  trk.getSurface()->getStartingError(err);
538  KETrack tre(trk, err);
539  KFitTrack trf(tre, path, tchisq);
540 
541  // Also use the starting track as the reference track for linearized
542  // propagation, until the track is established with reasonably small
543  // errors.
544 
545  KTrack ref(trk);
546  KTrack* pref = &ref;
547 
548  mf::LogInfo log("KalmanFilterAlg");
549 
550  // Loop over measurement groups (KHitGroups).
551 
552  while (hits.getSorted().size() > 0) {
553  ++step;
554  if (fTrace) {
555  log << "Build Step " << step << "\n";
556  log << "KGTrack has " << trg.numHits() << " hits.\n";
557  log << trf;
558  }
559 
560  // Get an iterator for the next KHitGroup.
561 
563  if (dir == Propagator::FORWARD)
564  it = hits.getSorted().begin();
565  else {
566  it = hits.getSorted().end();
567  --it;
568  }
569  const KHitGroup& gr = *it;
570 
571  if (fTrace) {
572  double path_est = gr.getPath();
573  log << "Next surface: " << *(gr.getSurface()) << "\n";
574  log << " Estimated path length of hit group = " << path_est << "\n";
575  }
576 
577  // Get the next prediction surface. If this KHitGroup is on the
578  // preferred plane, use that as the prediction surface.
579  // Otherwise, use the current track surface as the prediction
580  // surface.
581 
582  std::shared_ptr<const Surface> psurf = trf.getSurface();
583  if (gr.getPlane() < 0) throw cet::exception("KalmanFilterAlg") << "negative plane?\n";
584  if (fPlane < 0 || gr.getPlane() < 0 || fPlane == gr.getPlane()) psurf = gr.getSurface();
585 
586  // Propagate track to the prediction surface.
587 
588  std::optional<double> dist = prop.noise_prop(trf, psurf, Propagator::UNKNOWN, true, pref);
589  if (dist && std::abs(*dist) > fMaxPropDist) dist = std::nullopt;
590  double ds = 0.;
591 
592  if (dist) {
593 
594  // Propagation succeeded.
595  // Update cumulative path distance and track status.
596 
597  ds = *dist;
598  path += ds;
599  trf.setPath(path);
600  if (dir == Propagator::FORWARD)
602  else {
604  }
605  if (fTrace) {
606  log << "After propagation\n";
607  log << " Incremental propagation distance = " << ds << "\n";
608  log << " Path length of prediction surface = " << path << "\n";
609  log << "KGTrack has " << trg.numHits() << " hits.\n";
610  log << trf;
611  }
612 
613  // Loop over measurements in this group.
614 
615  const std::vector<std::shared_ptr<const KHitBase>>& hits = gr.getHits();
616  double best_chisq = 0.;
617  std::shared_ptr<const KHitBase> best_hit;
618  for (std::vector<std::shared_ptr<const KHitBase>>::const_iterator ihit = hits.begin();
619  ihit != hits.end();
620  ++ihit) {
621  const KHitBase& hit = **ihit;
622 
623  // Turn this hit blue.
624 
625  if (fGTrace && fCanvases.size() > 0) {
626  auto marker_it = fMarkerMap.find(hit.getID());
627  if (marker_it != fMarkerMap.end()) {
628  TMarker* marker = marker_it->second;
629  marker->SetMarkerColor(kBlue);
630  }
631  //fCanvases.back()->Update();
632  }
633 
634  // Update predction using current track hypothesis and get
635  // incremental chisquare.
636 
637  if (hit.predict(trf, prop)) {
638  double chisq = hit.getChisq();
639  double preddist = hit.getPredDistance();
640  if (fTrace) {
641  log << "Trying Hit.\n"
642  << hit << "\nchisq = " << chisq << "\n"
643  << "prediction distance = " << preddist << "\n";
644  }
645  if ((!has_pref_plane || abs(preddist) < fMaxSeedPredDist) &&
646  (best_hit.get() == 0 || chisq < best_chisq)) {
647  best_chisq = chisq;
648  if (chisq < fMaxSeedIncChisq) best_hit = *ihit;
649  }
650  }
651  }
652  if (fTrace) {
653  log << "Best hit incremental chisquare = " << best_chisq << "\n";
654  if (best_hit.get() != 0) {
655  log << "Hit after prediction\n";
656  log << *best_hit;
657  }
658  else
659  log << "No hit passed chisquare cut.\n";
660  }
661  if (fGTrace && fCanvases.size() > 0) fCanvases.back()->Update();
662 
663  // If we found a best measurement, and if the incremental
664  // chisquare passes the cut, add it to the track and update
665  // fit information.
666 
667  bool update_ok = false;
668  if (best_hit.get() != 0) {
669  KFitTrack trf0(trf);
670  best_hit->update(trf);
671  update_ok = trf.isValid();
672  if (!update_ok) trf = trf0;
673  }
674  if (update_ok) {
675  ds += best_hit->getPredDistance();
676  tchisq += best_chisq;
677  trf.setChisq(tchisq);
678  if (dir == Propagator::FORWARD)
680  else {
682  }
683 
684  // If the pointing error got too large, and there is no
685  // reference track, quit.
686 
687  if (pref == 0 && trf.PointingError() > fMaxPErr) {
688  if (fTrace) log << "Quitting because pointing error got too large.\n";
689  break;
690  }
691 
692  // Test number of consecutive measurements in the same plane.
693 
694  if (gr.getPlane() >= 0) {
695  if (gr.getPlane() == last_plane)
696  ++nsame;
697  else {
698  nsame = 1;
699  last_plane = gr.getPlane();
700  }
701  }
702  else {
703  nsame = 0;
704  last_plane = -1;
705  }
706  if (nsame <= fMaxSamePlane) {
707 
708  // Turn best hit red.
709 
710  if (fGTrace && fCanvases.size() > 0) {
711  int pl = best_hit->getMeasPlane();
712  if (pl >= 0 && pl < int(fPads.size())) {
713  auto marker_it = fMarkerMap.find(best_hit->getID());
714  if (marker_it != fMarkerMap.end()) {
715  TMarker* marker = marker_it->second;
716  marker->SetMarkerColor(kRed);
717 
718  // Redraw marker so that it will be on top.
719 
720  fPads[pl]->cd();
721  marker->Draw();
722  }
723  }
724  fCanvases.back()->Update();
725  }
726 
727  // Make a KHitTrack and add it to the KGTrack.
728 
729  KHitTrack trh(trf, best_hit);
730  trg.addTrack(trh);
731  if (fPlane == gr.getPlane()) has_pref_plane = true;
732 
733  // Decide if we want to kill the reference track.
734 
735  if (!linear && pref != 0 && int(trg.numHits()) >= fMinLHits &&
736  (trf.PointingError() < fGoodPErr || path > fMaxLDist)) {
737  pref = 0;
738  if (fTrace) log << "Killing reference track.\n";
739  }
740 
741  if (fTrace) {
742  log << "After update\n";
743  log << "KGTrack has " << trg.numHits() << " hits.\n";
744  log << trf;
745  }
746  }
747  }
748  }
749 
750  // The current KHitGroup is now resolved.
751  // Move it to unused list.
752 
753  hits.getUnused().splice(hits.getUnused().end(), hits.getSorted(), it);
754 
755  // If the propagation distance was the wrong direction, resort the measurements.
756 
757  if (pref == 0 && dist &&
758  ((dir == Propagator::FORWARD && (ds < fMinSortDist || ds > fMaxSortDist)) ||
759  (dir == Propagator::BACKWARD && (-ds < fMinSortDist || -ds > fMaxSortDist)))) {
760  if (fTrace) log << "Resorting measurements.\n";
761  hits.sort(trf, true, prop, dir);
762  }
763  }
764 
765  // Clean track.
766 
767  cleanTrack(trg);
768 
769  // Set the fit status of the last added KHitTrack to optimal and get
770  // the final chisquare and path length.
771 
772  double fchisq = 0.;
773  path = 0.;
774  if (trg.isValid()) {
775  path = trg.endTrack().getPath() - trg.startTrack().getPath();
776  if (dir == Propagator::FORWARD) {
778  fchisq = trg.endTrack().getChisq();
779  }
780  else {
782  fchisq = trg.startTrack().getChisq();
783  }
784  }
785 
786  // Summary.
787 
788  log << "KalmanFilterAlg build track summary.\n"
789  << "Build direction = " << (dir == Propagator::FORWARD ? "FORWARD" : "BACKWARD") << "\n"
790  << "Track has " << trg.numHits() << " hits.\n"
791  << "Track length = " << path << "\n"
792  << "Track chisquare = " << fchisq << "\n";
793 
794  // Done. Return success if we added at least one measurement.
795 
796  bool ok = trg.numHits() > 0;
797 
798  if (fGTrace && fCanvases.size() > 0) {
799  TText* t = 0;
800  if (ok)
801  t = fMessages->AddText("Exit buildTrack, status success");
802  else
803  t = fMessages->AddText("Exit buildTrack, status fail");
804  t->SetBit(kCanDelete);
805  fInfoPad->Modified();
806  fCanvases.back()->Update();
807  }
808 
809  return ok;
810 }
811 
845 bool trkf::KalmanFilterAlg::smoothTrack(KGTrack& trg, KGTrack* trg1, const Propagator& prop) const
846 {
847  if (not trg.isValid()) {
848  // It is an error if the KGTrack is not valid.
849  return false;
850  }
851 
852  bool result = false;
853 
854  // Examine the track endpoints and figure out which end of the track
855  // to start from. The fit always starts at the optimal end. It is
856  // an error if neither end point is optimal. Do nothing and return
857  // success if both end points are optimal.
858 
859  const KHitTrack& trh0 = trg.startTrack();
860  const KHitTrack& trh1 = trg.endTrack();
861  KFitTrack::FitStatus stat0 = trh0.getStat();
862  KFitTrack::FitStatus stat1 = trh1.getStat();
863  bool dofit = false;
864 
865  // Remember starting direction, track, and distance.
866 
868  const KTrack* trk = 0;
869  double path = 0.;
870 
871  if (stat0 == KFitTrack::OPTIMAL) {
872  if (stat1 == KFitTrack::OPTIMAL) {
873 
874  // Both ends optimal (do nothing, return success).
875 
876  dofit = false;
877  result = true;
878  }
879  else {
880 
881  // Start optimal.
882 
883  dofit = true;
884  dir = Propagator::FORWARD;
885  trk = &trh0;
886  path = 0.;
887  }
888  }
889  else {
890  if (stat1 == KFitTrack::OPTIMAL) {
891 
892  // End optimal.
893 
894  dofit = true;
895  dir = Propagator::BACKWARD;
896  trk = &trh1;
897  path = trh1.getPath();
898  }
899  else {
900 
901  // Neither end optimal (do nothing, return failure).
902 
903  dofit = false;
904  result = false;
905  }
906  }
907  if (dofit) {
908  if (!trk)
909  throw cet::exception("KalmanFilterAlg") << "KalmanFilterAlg::smoothTrack(): no track!\n";
910 
911  // Cumulative chisquare.
912 
913  double tchisq = 0.;
914 
915  // Construct starting KFitTrack with track information and distance
916  // taken from the optimal end, but with "infinite" errors.
917 
918  TrackError err;
919  trk->getSurface()->getStartingError(err);
920  KETrack tre(*trk, err);
921  KFitTrack trf(tre, path, tchisq);
922 
923  // Make initial reference track to be same as initial fit track.
924 
925  KTrack ref(trf);
926 
927  // Loop over KHitTracks contained in KGTrack.
928 
931  switch (dir) {
932  case Propagator::FORWARD:
933  it = trg.getTrackMap().begin();
934  itend = trg.getTrackMap().end();
935  break;
937  it = trg.getTrackMap().end();
938  itend = trg.getTrackMap().begin();
939  break;
940  default:
941  throw cet::exception("KalmanFilterAlg")
942  << "KalmanFilterAlg::smoothTrack(): invalid direction\n";
943  } // switch
944 
945  mf::LogInfo log("KalmanFilterAlg");
946 
947  // Loop starts here.
948 
949  result = true; // Result success unless we find an error.
950  int step = 0; // Step count.
951  while (dofit && it != itend) {
952  ++step;
953  if (fTrace) {
954  log << "Smooth Step " << step << "\n";
955  log << "Reverse fit track:\n";
956  log << trf;
957  }
958 
959  // For backward fit, decrement iterator at start of loop.
960 
961  if (dir == Propagator::BACKWARD) --it;
962 
963  KHitTrack& trh = (*it).second;
964  if (fTrace) {
965  log << "Forward track:\n";
966  log << trh;
967  }
968 
969  // Extract measurement.
970 
971  const KHitBase& hit = *(trh.getHit());
972 
973  // Propagate KFitTrack to the next track surface.
974 
975  std::shared_ptr<const Surface> psurf = trh.getSurface();
976  std::optional<double> dist = prop.noise_prop(trf, psurf, Propagator::UNKNOWN, true, &ref);
977 
978  // Check if propagation succeeded. If propagation fails, this
979  // measurement will be dropped from the unidirectional fit
980  // track. This measurement will still be in the original
981  // track, but with a status other than optimal.
982 
983  if (dist) {
984 
985  // Propagation succeeded.
986  // Update cumulative path distance and track status.
987 
988  double ds = *dist;
989  path += ds;
990  trf.setPath(path);
991  if (dir == Propagator::FORWARD)
993  else {
995  }
996  if (fTrace) {
997  log << "Reverse fit track after propagation:\n";
998  log << " Propagation distance = " << ds << "\n";
999  log << trf;
1000  }
1001 
1002  // See if we have the proper information to calculate an optimal track
1003  // at this surface (should normally be possible).
1004 
1005  KFitTrack::FitStatus stat = trh.getStat();
1006  KFitTrack::FitStatus newstat = trf.getStat();
1007 
1008  if ((newstat == KFitTrack::FORWARD_PREDICTED && stat == KFitTrack::BACKWARD) ||
1009  (newstat == KFitTrack::BACKWARD_PREDICTED && stat == KFitTrack::FORWARD)) {
1010 
1011  // Update stored KHitTrack to be optimal.
1012 
1013  bool ok = trh.combineFit(trf);
1014 
1015  // Update the stored path distance to be from the currently fitting track.
1016 
1017  trh.setPath(trf.getPath());
1018 
1019  // Update reference track.
1020 
1021  ref = trh;
1022 
1023  // If combination failed, abandon the fit and return failure.
1024 
1025  if (!ok) {
1026  dofit = false;
1027  result = false;
1028  break;
1029  }
1030  if (fTrace) {
1031  log << "Combined track:\n";
1032  log << trh;
1033  }
1034  }
1035 
1036  // Update measurement predction using current track hypothesis.
1037 
1038  if (not hit.predict(trf, prop, &ref)) {
1039 
1040  // Abandon the fit and return failure.
1041 
1042  dofit = false;
1043  result = false;
1044  break;
1045  }
1046 
1047  // Prediction succeeded. Get incremental chisquare. If this
1048  // hit fails incremental chisquare cut, this hit will be
1049  // dropped from the unidirecitonal Kalman fit track, but may
1050  // still be in the smoothed track.
1051 
1052  double chisq = hit.getChisq();
1053  if (chisq < fMaxSmoothIncChisq) {
1054 
1055  // Update the reverse fitting track using the current measurement
1056  // (both track parameters and status).
1057 
1058  KFitTrack trf0(trf);
1059  hit.update(trf);
1060  bool update_ok = trf.isValid();
1061  if (!update_ok)
1062  trf = trf0;
1063  else {
1064  tchisq += chisq;
1065  trf.setChisq(tchisq);
1066 
1067  if (dir == Propagator::FORWARD)
1069  else {
1071  }
1072  if (fTrace) {
1073  log << "Reverse fit track after update:\n";
1074  log << trf;
1075  }
1076  if (fGTrace && fCanvases.size() > 0) {
1077  auto marker_it = fMarkerMap.find(hit.getID());
1078  if (marker_it != fMarkerMap.end()) {
1079  TMarker* marker = marker_it->second;
1080  marker->SetMarkerColor(kOrange);
1081  }
1082  fCanvases.back()->Update();
1083  }
1084 
1085  // If unidirectional track pointer is not null, make a
1086  // KHitTrack and save it in the unidirectional track.
1087 
1088  if (trg1 != 0) {
1089  KHitTrack trh1(trf, trh.getHit());
1090  trg1->addTrack(trh1);
1091  }
1092  }
1093  }
1094  }
1095 
1096  // For forward fit, increment iterator at end of loop.
1097 
1098  if (dir == Propagator::FORWARD) ++it;
1099 
1100  } // Loop over KHitTracks.
1101 
1102  // If fit was successful and the unidirectional track pointer
1103  // is not null and the track is valid, set the fit status of
1104  // the last added KHitTrack to optimal.
1105 
1106  if (result && trg1 != 0 && trg1->isValid()) {
1107  if (dir == Propagator::FORWARD)
1109  else {
1111  }
1112  }
1113 
1114  // Recalibrate track map.
1115 
1116  trg.recalibrate();
1117 
1118  } // Do fit.
1119 
1120  // Get the final chisquare.
1121 
1122  double fchisq = 0.5 * (trg.startTrack().getChisq() + trg.endTrack().getChisq());
1123 
1124  // Summary.
1125 
1126  mf::LogInfo log("KalmanFilterAlg");
1127  log << "KalmanFilterAlg smooth track summary.\n"
1128  << "Smooth direction = " << (dir == Propagator::FORWARD ? "FORWARD" : "BACKWARD") << "\n"
1129  << "Track has " << trg.numHits() << " hits.\n"
1130  << "Track length = " << trg.endTrack().getPath() - trg.startTrack().getPath() << "\n"
1131  << "Track chisquare = " << fchisq << "\n";
1132 
1133  return result;
1134 }
1135 
1153  const Propagator& prop,
1154  KHitContainer& hits) const
1155 {
1156  // Default result failure.
1157 
1158  bool result = false;
1159 
1160  if (fGTrace && fCanvases.size() > 0) {
1161  TText* t = fMessages->AddText("Enter extendTrack");
1162  t->SetBit(kCanDelete);
1163  fInfoPad->Modified();
1164  fCanvases.back()->Update();
1165  }
1166 
1167  // Remember the original number of measurement.
1168 
1169  unsigned int nhits0 = trg.numHits();
1170 
1171  // It is an error if the KGTrack is not valid.
1172 
1173  if (trg.isValid()) {
1174  mf::LogInfo log("KalmanFilterAlg");
1175 
1176  // Examine the track endpoints and figure out which end of the
1177  // track to extend. The track is always extended from the optimal
1178  // end. It is an error if neither end point is optimal, or both
1179  // endoints are optimal. Reset the status of the optimal, and
1180  // make a copy of the starting track fit. Also get starting path
1181  // and chisquare.
1182 
1183  KHitTrack& trh0 = trg.startTrack();
1184  KHitTrack& trh1 = trg.endTrack();
1185  KFitTrack::FitStatus stat0 = trh0.getStat();
1186  KFitTrack::FitStatus stat1 = trh1.getStat();
1187  bool dofit = false;
1189  KFitTrack trf;
1190  double path = 0.;
1191  double tchisq = 0.;
1192 
1193  if (stat0 == KFitTrack::OPTIMAL) {
1194  if (stat1 == KFitTrack::OPTIMAL) {
1195 
1196  // Both ends optimal (do nothing, return failure).
1197 
1198  dofit = false;
1199  result = false;
1200  return result;
1201  }
1202  else {
1203 
1204  // Start optimal. Extend backward.
1205 
1206  dofit = true;
1207  dir = Propagator::BACKWARD;
1209  trf = trh0;
1210  path = trh0.getPath();
1211  tchisq = trh0.getChisq();
1212  }
1213  }
1214  else {
1215  if (stat1 == KFitTrack::OPTIMAL) {
1216 
1217  // End optimal. Extend forward.
1218 
1219  dofit = true;
1220  dir = Propagator::FORWARD;
1222  trf = trh1;
1223  path = trh1.getPath();
1224  tchisq = trh1.getChisq();
1225 
1226  // Make sure forward extend track momentum is over some
1227  // minimum value.
1228 
1229  if (trf.getVector()(4) > 5.) {
1230  trf.getVector()(4) = 5.;
1231  trf.getError()(4, 4) = 5.;
1232  }
1233  }
1234  else {
1235 
1236  // Neither end optimal (do nothing, return failure).
1237 
1238  dofit = false;
1239  result = false;
1240  return result;
1241  }
1242  }
1243  if (dofit) {
1244 
1245  // Sort hit container using starting track.
1246 
1247  hits.sort(trf, true, prop, dir);
1248 
1249  // Draw and add hits in hit->marker map that are not already there.
1250 
1251  if (fGTrace && fCanvases.size() > 0) {
1252 
1253  // Loop over sorted KHitGroups.
1254  // Paint sorted hits black.
1255 
1256  const std::list<KHitGroup>& groups = hits.getSorted();
1257  for (auto const& gr : groups) {
1258 
1259  // Loop over hits in this group.
1260 
1261  const std::vector<std::shared_ptr<const KHitBase>>& phits = gr.getHits();
1262  for (auto const& phit : phits) {
1263  const KHitBase& hit = *phit;
1264  int pl = hit.getMeasPlane();
1265  if (pl >= 0 && pl < int(fPads.size())) {
1266 
1267  // Is this hit already in map?
1268 
1269  TMarker* marker = 0;
1270  auto marker_it = fMarkerMap.find(hit.getID());
1271  if (marker_it == fMarkerMap.end()) {
1272  double z = 0.;
1273  double x = 0.;
1274  hit_position(hit, z, x);
1275  marker = new TMarker(z, x, 20);
1276  fMarkerMap[hit.getID()] = marker;
1277  fPads[pl]->cd();
1278  marker->SetBit(kCanDelete); // Give away ownership.
1279  marker->SetMarkerSize(0.5);
1280  marker->Draw();
1281  }
1282  else
1283  marker = marker_it->second;
1284  marker->SetMarkerColor(kBlack);
1285  }
1286  }
1287  }
1288 
1289  // Loop over unsorted KHitGroups.
1290  // Paint unsorted hits blue.
1291 
1292  const std::list<KHitGroup>& ugroups = hits.getUnsorted();
1293  for (auto const& gr : ugroups) {
1294 
1295  // Loop over hits in this group.
1296 
1297  const std::vector<std::shared_ptr<const KHitBase>>& phits = gr.getHits();
1298  for (auto const& phit : phits) {
1299  const KHitBase& hit = *phit;
1300  int pl = hit.getMeasPlane();
1301  if (pl >= 0 && pl < int(fPads.size())) {
1302 
1303  // Is this hit already in map?
1304 
1305  TMarker* marker = 0;
1306  auto marker_it = fMarkerMap.find(hit.getID());
1307  if (marker_it == fMarkerMap.end()) {
1308  double z = 0.;
1309  double x = 0.;
1310  hit_position(hit, z, x);
1311  marker = new TMarker(z, x, 20);
1312  fMarkerMap[hit.getID()] = marker;
1313  fPads[pl]->cd();
1314  marker->SetBit(kCanDelete); // Give away ownership.
1315  marker->SetMarkerSize(0.5);
1316  marker->Draw();
1317  }
1318  else
1319  marker = marker_it->second;
1320  marker->SetMarkerColor(kBlue);
1321  }
1322  }
1323  }
1324  fCanvases.back()->Update();
1325  }
1326 
1327  //std::cout << "extendTrack: marker map has " << fMarkerMap.size() << " entries." << std::endl;
1328 
1329  // Extend loop starts here.
1330 
1331  int step = 0;
1332  int nsame = 0;
1333  int last_plane = -1;
1334  while (hits.getSorted().size() > 0) {
1335  ++step;
1336  if (fTrace) {
1337  log << "Extend Step " << step << "\n";
1338  log << "KGTrack has " << trg.numHits() << " hits.\n";
1339  log << trf;
1340  }
1341 
1342  // Get an iterator for the next KHitGroup.
1343 
1345  switch (dir) {
1346  case Propagator::FORWARD: it = hits.getSorted().begin(); break;
1347  case Propagator::BACKWARD:
1348  it = hits.getSorted().end();
1349  --it;
1350  break;
1351  default:
1352  throw cet::exception("KalmanFilterAlg")
1353  << "KalmanFilterAlg::extendTrack(): invalid direction\n";
1354  } // switch
1355  const KHitGroup& gr = *it;
1356 
1357  if (fTrace) {
1358  double path_est = gr.getPath();
1359  log << "Next surface: " << *(gr.getSurface()) << "\n";
1360  log << " Estimated path length of hit group = " << path_est << "\n";
1361  }
1362 
1363  // Get the next prediction surface. If this KHitGroup is on the
1364  // preferred plane, use that as the prediction surface.
1365  // Otherwise, use the current track surface as the prediction
1366  // surface.
1367 
1368  std::shared_ptr<const Surface> psurf = trf.getSurface();
1369  if (gr.getPlane() < 0)
1370  throw cet::exception("KalmanFilterAlg")
1371  << "KalmanFilterAlg::extendTrack(): negative plane?\n";
1372  if (fPlane < 0 || gr.getPlane() < 0 || fPlane == gr.getPlane()) psurf = gr.getSurface();
1373 
1374  // Propagate track to the prediction surface.
1375 
1376  std::optional<double> dist = prop.noise_prop(trf, psurf, Propagator::UNKNOWN, true);
1377  if (dist && std::abs(*dist) > fMaxPropDist) dist = std::nullopt;
1378  double ds = 0.;
1379 
1380  if (dist) {
1381 
1382  // Propagation succeeded.
1383  // Update cumulative path distance and track status.
1384 
1385  ds = *dist;
1386  path += ds;
1387  trf.setPath(path);
1388  if (dir == Propagator::FORWARD)
1390  else {
1392  }
1393  if (fTrace) {
1394  log << "After propagation\n";
1395  log << " Incremental distance = " << ds << "\n";
1396  log << " Track path length = " << path << "\n";
1397  log << "KGTrack has " << trg.numHits() << " hits.\n";
1398  log << trf;
1399  }
1400 
1401  // Loop over measurements in this group.
1402 
1403  const std::vector<std::shared_ptr<const KHitBase>>& hits = gr.getHits();
1404  double best_chisq = 0.;
1405  std::shared_ptr<const KHitBase> best_hit;
1406  for (std::vector<std::shared_ptr<const KHitBase>>::const_iterator ihit = hits.begin();
1407  ihit != hits.end();
1408  ++ihit) {
1409  const KHitBase& hit = **ihit;
1410 
1411  // Turn this hit blue.
1412 
1413  if (fGTrace && fCanvases.size() > 0) {
1414  auto marker_it = fMarkerMap.find(hit.getID());
1415  if (marker_it != fMarkerMap.end()) {
1416  TMarker* marker = marker_it->second;
1417  marker->SetMarkerColor(kBlue);
1418  }
1419  //fCanvases.back()->Update();
1420  }
1421 
1422  // Update predction using current track hypothesis and get
1423  // incremental chisquare.
1424 
1425  bool ok = hit.predict(trf, prop);
1426  if (ok) {
1427  double chisq = hit.getChisq();
1428  double preddist = hit.getPredDistance();
1429  if (abs(preddist) < fMaxPredDist && (best_hit.get() == 0 || chisq < best_chisq)) {
1430  best_chisq = chisq;
1431  if (chisq < fMaxIncChisq) best_hit = *ihit;
1432  }
1433  }
1434  }
1435  if (fTrace) {
1436  log << "Best hit incremental chisquare = " << best_chisq << "\n";
1437  if (best_hit.get() != 0) {
1438  log << "Hit after prediction\n";
1439  log << *best_hit;
1440  }
1441  else
1442  log << "No hit passed chisquare cut.\n";
1443  }
1444  if (fGTrace && fCanvases.size() > 0) fCanvases.back()->Update();
1445 
1446  // If we found a best measurement, and if the incremental
1447  // chisquare passes the cut, add it to the track and update
1448  // fit information.
1449 
1450  bool update_ok = false;
1451  if (best_hit.get() != 0) {
1452  KFitTrack trf0(trf);
1453  best_hit->update(trf);
1454  update_ok = trf.isValid();
1455  if (!update_ok) trf = trf0;
1456  }
1457  if (update_ok) {
1458  ds += best_hit->getPredDistance();
1459  tchisq += best_chisq;
1460  trf.setChisq(tchisq);
1461  if (dir == Propagator::FORWARD)
1463  else {
1465  }
1466 
1467  // If the pointing error got too large, quit.
1468 
1469  if (trf.PointingError() > fMaxPErr) {
1470  if (fTrace) log << "Quitting because pointing error got too large.\n";
1471  break;
1472  }
1473 
1474  // Test number of consecutive measurements in the same plane.
1475 
1476  if (gr.getPlane() >= 0) {
1477  if (gr.getPlane() == last_plane)
1478  ++nsame;
1479  else {
1480  nsame = 1;
1481  last_plane = gr.getPlane();
1482  }
1483  }
1484  else {
1485  nsame = 0;
1486  last_plane = -1;
1487  }
1488  if (nsame <= fMaxSamePlane) {
1489 
1490  // Turn best hit red.
1491 
1492  if (fGTrace && fCanvases.size() > 0) {
1493  int pl = best_hit->getMeasPlane();
1494  if (pl >= 0 && pl < int(fPads.size())) {
1495  auto marker_it = fMarkerMap.find(best_hit->getID());
1496  if (marker_it != fMarkerMap.end()) {
1497  TMarker* marker = marker_it->second;
1498  marker->SetMarkerColor(kRed);
1499 
1500  // Redraw marker so that it will be on top.
1501 
1502  fPads[pl]->cd();
1503  marker->Draw();
1504  }
1505  }
1506  fCanvases.back()->Update();
1507  }
1508 
1509  // Make a KHitTrack and add it to the KGTrack.
1510 
1511  KHitTrack trh(trf, best_hit);
1512  trg.addTrack(trh);
1513 
1514  if (fTrace) {
1515  log << "After update\n";
1516  log << "KGTrack has " << trg.numHits() << " hits.\n";
1517  log << trf;
1518  }
1519  }
1520  }
1521  }
1522 
1523  // The current KHitGroup is now resolved.
1524  // Move it to unused list.
1525 
1526  hits.getUnused().splice(hits.getUnused().end(), hits.getSorted(), it);
1527 
1528  // If the propagation distance was the wrong direction, resort the measurements.
1529 
1530  if (dist && ((dir == Propagator::FORWARD && (ds < fMinSortDist || ds > fMaxSortDist)) ||
1531  (dir == Propagator::BACKWARD && (-ds < fMinSortDist || -ds > fMaxSortDist)))) {
1532  if (fTrace) log << "Resorting measurements.\n";
1533  hits.sort(trf, true, prop, dir);
1534  }
1535  }
1536  }
1537 
1538  // Clean track.
1539 
1540  cleanTrack(trg);
1541 
1542  // Set the fit status of the last added KHitTrack to optimal and
1543  // get the final chisquare and path length.
1544 
1545  double fchisq = 0.;
1546  path = 0.;
1547  if (trg.isValid()) {
1548  path = trg.endTrack().getPath() - trg.startTrack().getPath();
1549  switch (dir) {
1550  case Propagator::FORWARD:
1552  fchisq = trg.endTrack().getChisq();
1553  break;
1554  case Propagator::BACKWARD:
1556  fchisq = trg.startTrack().getChisq();
1557  break;
1558  default:
1559  throw cet::exception("KalmanFilterAlg")
1560  << "KalmanFilterAlg::extendTrack(): invalid direction [II]\n";
1561  } // switch
1562  }
1563 
1564  // Summary.
1565 
1566  log << "KalmanFilterAlg extend track summary.\n"
1567  << "Extend direction = " << (dir == Propagator::FORWARD ? "FORWARD" : "BACKWARD") << "\n"
1568  << "Track has " << trg.numHits() << " hits.\n"
1569  << "Track length = " << path << "\n"
1570  << "Track chisquare = " << fchisq << "\n";
1571  // log << trg << "\n";
1572  }
1573 
1574  // Done.
1575 
1576  result = (trg.numHits() > nhits0);
1577 
1578  if (fGTrace && fCanvases.size() > 0) {
1579  TText* t = 0;
1580  if (result)
1581  t = fMessages->AddText("Exit extendTrack, status success");
1582  else
1583  t = fMessages->AddText("Exit extendTrack, status fail");
1584  t->SetBit(kCanDelete);
1585  fInfoPad->Modified();
1586  fCanvases.back()->Update();
1587  }
1588 
1589  return result;
1590 }
1591 
1607 {
1608  if (!trg.isValid()) return false;
1609 
1610  // Extract track with lowest momentum.
1611 
1612  const KHitTrack& trh = trg.endTrack();
1613  if (trh.getStat() == KFitTrack::INVALID)
1614  throw cet::exception("KalmanFilterAlg")
1615  << "KalmanFilterAlg::fitMomentumRange(): invalid end track\n";
1616  tremom = trh;
1617 
1618  // Set track momentum to a small value.
1619 
1620  tremom.getVector()(4) = 100.;
1621  tremom.getError()(4, 0) = 0.;
1622  tremom.getError()(4, 1) = 0.;
1623  tremom.getError()(4, 2) = 0.;
1624  tremom.getError()(4, 3) = 0.;
1625  tremom.getError()(4, 4) = 10000.;
1626 
1627  // Done.
1628 
1629  return true;
1630 }
1631 
1660  const Propagator& prop,
1661  KETrack& tremom) const
1662 {
1663  // Get iterators pointing to the first and last tracks.
1664 
1665  const std::multimap<double, KHitTrack>& trackmap = trg.getTrackMap();
1666  if (trackmap.size() < 2) return false;
1668  itend[0] = trackmap.begin();
1669  itend[1] = trackmap.end();
1670  --itend[1];
1671 
1672  // Check the fit status and error matrix of the first and last
1673  // track.
1674 
1675  bool result = true;
1676 
1677  for (int i = 0; result && i < 2; ++i) {
1678  const KHitTrack& trh = itend[i]->second;
1679  KFitTrack::FitStatus stat = trh.getStat();
1680  if (stat != KFitTrack::OPTIMAL) result = false;
1681  const TrackError& err = trh.getError();
1682  for (int j = 0; j < 4; ++j) {
1683  if (err(4, j) != 0.) result = false;
1684  }
1685  }
1686  if (!result) return result;
1687 
1688  // We will periodically sample the track trajectory. At each sample
1689  // point, collect the following information.
1690  //
1691  // 1. The path distance at the sample point.
1692  // 2. The original momentum of the track at the sample point and its error.
1693  // 3. One copy of the track (KETrack) that will be propagated without noise
1694  // (infinite momentum track).
1695  // 3. A second copy of the track (KETrack) that will be propagated with the
1696  // minimum allowed momentum (range out track).
1697  // 4. A third copy of the track (KETrack) that will propagated with noise
1698  // with some intermediate momentum (noise track).
1699  //
1700  // Collect the first sample from the maximum path distance track.
1701 
1702  double s_sample = itend[1]->first;
1703  const KETrack& tre = itend[1]->second;
1704  KETrack tre_inf(tre);
1705  KTrack trk_range(tre);
1706  KETrack tre_noise(tre);
1707  tre_inf.getVector()(4) = 0.;
1708  tre_inf.getError()(4, 4) = 0.;
1709  trk_range.getVector()(4) = 100.;
1710  tre_noise.getError()(4, 4) = 0.;
1711  tre_noise.getVector()(4) = 1.;
1712  tre_noise.getError()(4, 4) = 10.;
1713  double invp0 = tre_noise.getVector()(4);
1714  double var_invp0 = tre_noise.getError()(4, 4);
1715 
1716  // Loop over fits, starting at the high path distance (low momentum)
1717  // end.
1718 
1719  for (auto it = trackmap.rbegin(); it != trackmap.rend(); ++it) {
1720  double s = it->first;
1721  const KHitTrack& trh = it->second;
1722 
1723  // Ignore non-optimal fits.
1724 
1725  KFitTrack::FitStatus stat = trh.getStat();
1726  if (stat != KFitTrack::OPTIMAL) continue;
1727 
1728  // See if this track is far enough from the previous sample to
1729  // make a new sample point.
1730 
1731  if (std::abs(s - s_sample) > fMinSampleDist) {
1732 
1733  // Propagate tracks to the current track surface.
1734 
1735  std::shared_ptr<const Surface> psurf = trh.getSurface();
1736  std::optional<double> dist_inf = prop.err_prop(tre_inf, psurf, Propagator::UNKNOWN, false);
1737  std::optional<double> dist_range =
1738  prop.vec_prop(trk_range, psurf, Propagator::UNKNOWN, false);
1739  std::optional<double> dist_noise =
1740  prop.noise_prop(tre_noise, psurf, Propagator::UNKNOWN, true);
1741 
1742  // All propagations should normally succeed. If they don't,
1743  // ignore this sample for the purpose of updating the momentum.
1744 
1745  bool momentum_updated = false;
1746  if (dist_inf && dist_range && dist_noise) {
1747 
1748  // Extract the momentum at the new sample point.
1749 
1750  double invp1 = tre_noise.getVector()(4);
1751  double var_invp1 = tre_noise.getError()(4, 4);
1752 
1753  // Get the average momentum and error for this pair of
1754  // sample points, and other data.
1755 
1756  double invp = 0.5 * (invp0 + invp1);
1757  double var_invp = 0.5 * (var_invp0 + var_invp1);
1758  double mass = tre_inf.Mass();
1759  double beta = std::sqrt(1. + mass * mass * invp * invp);
1760  double invbp = invp / beta;
1761 
1762  // Extract slope subvectors and sub-error-matrices.
1763  // We have the following variables.
1764  //
1765  // slope0 - Predicted slope vector (from infinite momentum track).
1766  // slope1 - Measured slope vector (from new sample point).
1767  // defl - Deflection (slope residual = difference between measured
1768  // and predicted slope vector).
1769  // err0 - Slope error matrix of prediction.
1770  // err1 - Slope error matrix of measurement
1771  // errc - Slope residual error matrix = err0 + err1.
1772  // errn - Noise slope error matrix divided by (1/beta*momentum).
1773 
1774  KVector<2>::type slope0 = project(tre_inf.getVector(), ublas::range(2, 4));
1775  KVector<2>::type slope1 = project(trh.getVector(), ublas::range(2, 4));
1776  KVector<2>::type defl = slope1 - slope0;
1777  KSymMatrix<2>::type err0 =
1778  project(tre_inf.getError(), ublas::range(2, 4), ublas::range(2, 4));
1779  KSymMatrix<2>::type err1 = project(trh.getError(), ublas::range(2, 4), ublas::range(2, 4));
1780  KSymMatrix<2>::type errc = err0 + err1;
1781  KSymMatrix<2>::type errn =
1782  project(tre_noise.getError(), ublas::range(2, 4), ublas::range(2, 4));
1783  errn -= err0;
1784  errn /= invbp;
1785 
1786  // Calculate updated average momentum and error.
1787 
1788  double new_invp = invp;
1789  double new_var_invp = var_invp;
1790  update_momentum(defl, errc, errn, mass, new_invp, new_var_invp);
1791 
1792  // Calculate updated momentum and error at the second sample
1793  // point.
1794 
1795  double dp = 1. / new_invp - 1. / invp;
1796  invp0 = 1. / (1. / invp1 + dp);
1797  var_invp0 = new_var_invp;
1798  momentum_updated = true;
1799 
1800  // Make sure that updated momentum is not less than minimum
1801  // allowed momentum.
1802 
1803  double invp_range = trk_range.getVector()(4);
1804  if (invp0 > invp_range) invp0 = invp_range;
1805  }
1806 
1807  // Update sample.
1808 
1809  if (momentum_updated) {
1810  s_sample = s;
1811  tre_inf = trh;
1812  tre_inf.getVector()(4) = 0.;
1813  tre_inf.getError()(4, 4) = 0.;
1814  double invp_range = trk_range.getVector()(4);
1815  trk_range = trh;
1816  trk_range.getVector()(4) = invp_range;
1817  tre_noise = trh;
1818  tre_noise.getVector()(4) = invp0;
1819  tre_noise.getError()(4, 4) = var_invp0;
1820  }
1821  }
1822  }
1823 
1824  // Propagate noise track to starting (high momentum) track surface
1825  // to get final starting momentum. This propagation should normally
1826  // always succeed, but if it doesn't, don't update the track.
1827 
1828  const KHitTrack& trh0 = itend[0]->second;
1829  std::shared_ptr<const Surface> psurf = trh0.getSurface();
1830  std::optional<double> dist_noise = prop.noise_prop(tre_noise, psurf, Propagator::UNKNOWN, true);
1831  result = dist_noise.has_value();
1832 
1833  // Update momentum-estimating track.
1834 
1835  mf::LogInfo log("KalmanFilterAlg");
1836  if (result) tremom = tre_noise;
1837 
1838  // Done.
1839 
1840  return result;
1841 }
1842 
1854  const Propagator& prop,
1855  KETrack& tremom) const
1856 {
1857  mf::LogInfo log("KalmanFilterAlg");
1858  double invp_range = 0.;
1859  double invp_ms = 0.;
1860 
1861  // Get multiple scattering momentum estimate.
1862 
1863  KETrack tremom_ms;
1864  bool ok_ms = false;
1865  if (fFitMomMS) {
1866  ok_ms = fitMomentumMS(trg, prop, tremom_ms);
1867  if (ok_ms) {
1868  KGTrack trg_ms(trg);
1869  ok_ms = updateMomentum(tremom_ms, prop, trg_ms);
1870  if (ok_ms) {
1871  invp_ms = trg_ms.startTrack().getVector()(4);
1872  double var_invp = trg_ms.startTrack().getError()(4, 4);
1873  double p = 0.;
1874  if (invp_ms != 0.) p = 1. / invp_ms;
1875  double err_p = p * p * std::sqrt(var_invp);
1876  log << "Multiple scattering momentum estimate = " << p << "+-" << err_p << "\n";
1877  }
1878  }
1879  }
1880 
1881  // Get range momentum estimate.
1882 
1883  KETrack tremom_range;
1884  bool ok_range = false;
1885  if (fFitMomRange) {
1886  ok_range = fitMomentumRange(trg, tremom_range);
1887  if (ok_range) {
1888  KGTrack trg_range(trg);
1889  ok_range = updateMomentum(tremom_range, prop, trg_range);
1890  if (ok_range) {
1891  invp_range = trg_range.startTrack().getVector()(4);
1892  double var_invp = trg_range.startTrack().getError()(4, 4);
1893  double p = 0.;
1894  if (invp_range != 0.) p = 1. / invp_range;
1895  double err_p = p * p * std::sqrt(var_invp);
1896  log << "Range momentum estimate = " << p << "+-" << err_p << "\n";
1897  }
1898  }
1899  }
1900 
1901  bool result = false;
1902  if (ok_range) {
1903  tremom = tremom_range;
1904  result = ok_range;
1905  }
1906  else if (ok_ms) {
1907  tremom = tremom_ms;
1908  result = ok_ms;
1909  }
1910  return result;
1911 }
1912 
1935  const Propagator& prop,
1936  KGTrack& trg) const
1937 {
1938  // Get modifiable track map.
1939 
1940  std::multimap<double, KHitTrack>& trackmap = trg.getTrackMap();
1941 
1942  // If track map is empty, immediately return failure.
1943 
1944  if (trackmap.size() == 0) return false;
1945 
1946  // Make trial propagations to the first and last track fit to figure
1947  // out which track fit is closer to the momentum estimating track.
1948 
1949  // Find distance to first track fit.
1950 
1951  KETrack tre0(tremom);
1952  std::optional<double> dist0 =
1953  prop.vec_prop(tre0, trackmap.begin()->second.getSurface(), Propagator::UNKNOWN, false, 0, 0);
1954  // Find distance to last track fit.
1955 
1956  KETrack tre1(tremom);
1957  std::optional<double> dist1 =
1958  prop.vec_prop(tre1, trackmap.rbegin()->second.getSurface(), Propagator::UNKNOWN, false, 0, 0);
1959 
1960  // Based on distances, make starting iterator and direction flag.
1961 
1962  bool forward = true;
1963  std::multimap<double, KHitTrack>::iterator it = trackmap.begin();
1964  if (dist0) {
1965 
1966  // Propagation to first track succeeded.
1967 
1968  if (dist1) {
1969 
1970  // Propagation to both ends succeeded. If the last track is
1971  // closer, initialize iterator and direction flag for reverse
1972  // direction.
1973 
1974  if (std::abs(*dist0) > std::abs(*dist1)) {
1975  it = trackmap.end();
1976  --it;
1977  forward = false;
1978  }
1979  }
1980  }
1981  else {
1982 
1983  // Propagation to first track failed. Initialize iterator and
1984  // direction flag for reverse direction, provided that the
1985  // propagation to the last track succeeded.
1986 
1987  if (dist1) {
1988  it = trackmap.end();
1989  --it;
1990  forward = false;
1991  }
1992  else {
1993 
1994  // Propagation to both ends failed. Return failure.
1995 
1996  return false;
1997  }
1998  }
1999 
2000  // Loop over track fits in global track.
2001 
2002  KETrack tre(tremom);
2003  bool done = false;
2004  while (!done) {
2005  KHitTrack& trh = it->second;
2006  if (trh.getStat() == KFitTrack::INVALID)
2007  throw cet::exception("KalmanFilterAlg")
2008  << "KalmanFilterAlg::updateMomentum(): invalid track\n";
2009 
2010  // Propagate momentum-estimating track to current track surface
2011  // and update momentum.
2012 
2013  std::optional<double> dist = prop.noise_prop(tre, trh.getSurface(), Propagator::UNKNOWN, true);
2014 
2015  // Copy momentum to global track.
2016 
2017  std::multimap<double, KHitTrack>::iterator erase_it = trackmap.end();
2018  if (dist) {
2019  trh.getVector()(4) = tre.getVector()(4);
2020  trh.getError()(4, 0) = 0.;
2021  trh.getError()(4, 1) = 0.;
2022  trh.getError()(4, 2) = 0.;
2023  trh.getError()(4, 3) = 0.;
2024  trh.getError()(4, 4) = tre.getError()(4, 4);
2025  }
2026  else {
2027 
2028  // If the propagation failed, remember that we are supposed to
2029  // erase this track from the global track.
2030 
2031  erase_it = it;
2032  }
2033 
2034  // Advance the iterator and set the done flag.
2035 
2036  if (forward) {
2037  ++it;
2038  done = (it == trackmap.end());
2039  }
2040  else {
2041  done = (it == trackmap.begin());
2042  if (!done) --it;
2043  }
2044 
2045  // Update momentum-estimating track from just-updated global track
2046  // fit, or erase global track fit.
2047 
2048  if (erase_it == trackmap.end())
2049  tre = trh;
2050  else
2051  trackmap.erase(erase_it);
2052  }
2053 
2054  return not empty(trackmap);
2055 }
2056 
2069 bool trkf::KalmanFilterAlg::smoothTrackIter(int nsmooth, KGTrack& trg, const Propagator& prop) const
2070 {
2071  bool ok = true;
2072 
2073  // Call smoothTrack method in a loop.
2074 
2075  for (int ismooth = 0; ok && ismooth < nsmooth - 1; ++ismooth) {
2076  KGTrack trg1{trg.getPrefPlane()};
2077  ok = smoothTrack(trg, &trg1, prop);
2078  if (ok) trg = trg1;
2079  }
2080 
2081  // Do one final smooth without generating a new unidirectional
2082  // track.
2083 
2084  if (ok) ok = smoothTrack(trg, 0, prop);
2085 
2086  // Done.
2087 
2088  return ok;
2089 }
2090 
2098 {
2099  // Get hold of a modifiable track map from trg.
2100 
2101  std::multimap<double, KHitTrack>& trackmap = trg.getTrackMap();
2102 
2103  // Do an indefinite loop until we no longer find any dirt.
2104 
2105  bool done = false;
2106  while (!done) {
2107 
2108  // If track map has fewer than fMaxNoiseHits tracks, then this is a
2109  // noise track. Clear the map, making the whole track invalid.
2110 
2111  int ntrack = trackmap.size();
2112  if (ntrack <= fMaxNoiseHits) {
2113  for (auto const& trackmap_entry : trackmap) {
2114  const KHitTrack& trh = trackmap_entry.second;
2115  const KHitBase& hit = *(trh.getHit());
2116  auto marker_it = fMarkerMap.find(hit.getID());
2117  if (marker_it != fMarkerMap.end()) {
2118  TMarker* marker = marker_it->second;
2119  marker->SetMarkerColor(kGreen);
2120  }
2121  }
2122  trackmap.clear();
2123  done = true;
2124  break;
2125  }
2126 
2127  // Make sure the first two and last two tracks belong to different
2128  // views. If not, remove the first or last track.
2129 
2130  if (ntrack >= 2) {
2131 
2132  // Check start.
2133 
2134  std::multimap<double, KHitTrack>::iterator it = trackmap.begin();
2135  const KHitTrack& trh1a = (*it).second;
2136  const KHitBase& hit1a = *(trh1a.getHit());
2137  int plane1 = hit1a.getMeasPlane();
2138  double chisq1 = hit1a.getChisq();
2139  ++it;
2140  const KHitTrack& trh2a = (*it).second;
2141  const KHitBase& hit2a = *(trh2a.getHit());
2142  int plane2 = hit2a.getMeasPlane();
2143  double chisq2 = hit2a.getChisq();
2144  if ((plane1 >= 0 && plane2 >= 0 && plane1 == plane2) || chisq1 > fMaxEndChisq ||
2145  chisq2 > fMaxEndChisq) {
2146  auto marker_it = fMarkerMap.find(hit1a.getID());
2147  if (marker_it != fMarkerMap.end()) {
2148  TMarker* marker = marker_it->second;
2149  marker->SetMarkerColor(kGreen);
2150  }
2151  trackmap.erase(trackmap.begin(), it);
2152  done = false;
2153  continue;
2154  }
2155 
2156  // Check end.
2157 
2158  it = trackmap.end();
2159  --it;
2160  const KHitTrack& trh1b = (*it).second;
2161  const KHitBase& hit1b = *(trh1b.getHit());
2162  plane1 = hit1b.getMeasPlane();
2163  chisq1 = hit1b.getChisq();
2164  --it;
2165  const KHitTrack& trh2b = (*it).second;
2166  const KHitBase& hit2b = *(trh2b.getHit());
2167  plane2 = hit2b.getMeasPlane();
2168  chisq2 = hit2b.getChisq();
2169  if ((plane1 >= 0 && plane2 >= 0 && plane1 == plane2) || chisq1 > fMaxEndChisq ||
2170  chisq2 > fMaxEndChisq) {
2171  ++it;
2172  auto marker_it = fMarkerMap.find(hit1b.getID());
2173  if (marker_it != fMarkerMap.end()) {
2174  TMarker* marker = marker_it->second;
2175  marker->SetMarkerColor(kGreen);
2176  }
2177  trackmap.erase(it, trackmap.end());
2178  done = false;
2179  continue;
2180  }
2181  }
2182 
2183  // Loop over successive pairs of elements of track map. Look for
2184  // adjacent elements with distance separation greater than
2185  // fGapDist.
2186 
2187  std::multimap<double, KHitTrack>::iterator it = trackmap.begin();
2188  std::multimap<double, KHitTrack>::iterator jt = trackmap.end();
2189  int nb = 0; // Number of elements from begin to jt.
2190  int ne = ntrack; // Number of elements it to end.
2191  bool found_noise = false;
2192  for (; it != trackmap.end(); ++it, ++nb, --ne) {
2193  if (jt == trackmap.end())
2194  jt = trackmap.begin();
2195  else {
2196  if (nb < 1)
2197  throw cet::exception("KalmanFilterAlg")
2198  << "KalmanFilterAlg::cleanTrack(): nb not positive\n";
2199  if (ne < 1)
2200  throw cet::exception("KalmanFilterAlg")
2201  << "KalmanFilterAlg::cleanTrack(): ne not positive\n";
2202 
2203  double disti = (*it).first;
2204  double distj = (*jt).first;
2205  double sep = disti - distj;
2206  if (sep < 0.)
2207  throw cet::exception("KalmanFilterAlg")
2208  << "KalmanFilterAlg::fitMomentumRange(): negative separation\n";
2209 
2210  if (sep > fGapDist) {
2211 
2212  // Found a gap. See if we want to trim track.
2213 
2214  if (nb <= fMaxNoiseHits) {
2215 
2216  // Trim front.
2217 
2218  found_noise = true;
2219  for (auto jt = trackmap.begin(); jt != it; ++jt) {
2220  const KHitTrack& trh = jt->second;
2221  const KHitBase& hit = *(trh.getHit());
2222  auto marker_it = fMarkerMap.find(hit.getID());
2223  if (marker_it != fMarkerMap.end()) {
2224  TMarker* marker = marker_it->second;
2225  marker->SetMarkerColor(kGreen);
2226  }
2227  }
2228  trackmap.erase(trackmap.begin(), it);
2229  break;
2230  }
2231  if (ne <= fMaxNoiseHits) {
2232 
2233  // Trim back.
2234 
2235  found_noise = true;
2236  for (auto jt = it; jt != trackmap.end(); ++jt) {
2237  const KHitTrack& trh = jt->second;
2238  const KHitBase& hit = *(trh.getHit());
2239  auto marker_it = fMarkerMap.find(hit.getID());
2240  if (marker_it != fMarkerMap.end()) {
2241  TMarker* marker = marker_it->second;
2242  marker->SetMarkerColor(kGreen);
2243  }
2244  }
2245  trackmap.erase(it, trackmap.end());
2246  break;
2247  }
2248  }
2249  ++jt;
2250  }
2251  }
2252 
2253  // Decide if we are done.
2254 
2255  done = !found_noise;
2256  }
2257  if (fGTrace && fCanvases.size() > 0) fCanvases.back()->Update();
2258 }
Float_t x
Definition: compare.C:6
bool fGTrace
Graphical trace flag.
intermediate_table::iterator iterator
int getPrefPlane() const
Definition: KGTrack.h:51
const TrackError & getError() const
Track error matrix.
Definition: KETrack.h:52
const std::shared_ptr< const KHitBase > & getHit() const
Measurement.
Definition: KHitTrack.h:51
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 Mass() const
Based on pdg code.
Definition: KTrack.cxx:116
const std::shared_ptr< const Surface > & getSurface() const
Surface.
Definition: KTrack.h:53
double PointingError() const
Pointing error (radians).
Definition: KETrack.cxx:66
double getPredDistance() const
Prediction distance.
Definition: KHitBase.h:79
void recalibrate()
Recalibrate track map.
Definition: KGTrack.cxx:88
std::map< int, TMarker * > fMarkerMap
Markers in current canvas.
int getMeasPlane() const
Measurement plane index.
Definition: KHitBase.h:85
const std::list< KHitGroup > & getUnused() const
Definition: KHitContainer.h:71
bool smoothTrackIter(int niter, KGTrack &trg, const Propagator &prop) const
Iteratively smooth a track.
const KHitTrack & endTrack() const
Track at end point.
Definition: KGTrack.cxx:40
KSymMatrix< 5 >::type TrackError
Track error matrix, dimension 5x5.
const std::vector< std::shared_ptr< const KHitBase > > & getHits() const
Measurement collection accessor.
Definition: KHitGroup.h:54
double fMaxEndChisq
Maximum incremental chisquare for endpoint hit.
bool fTrace
Trace flag.
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Planar surface parallel to x-axis.
KalmanFilterAlg(const fhicl::ParameterSet &pset)
Constructor.
Length_t DetHalfWidth(TPCID const &tpcid=tpc_zero) const
Returns the half width of the active volume of the specified TPC.
Double_t z
Definition: plot.C:276
A collection of KHitGroups.
constexpr auto abs(T v)
Returns the absolute value of the argument.
void cleanTrack(KGTrack &trg) const
Clean track by removing noise hits near endpoints.
void addTrack(const KHitTrack &trh)
Add track.
Definition: KGTrack.cxx:76
double getPath() const
Estimated path distance.
Definition: KHitGroup.h:60
intermediate_table::const_iterator const_iterator
FitStatus
Fit status enum.
Definition: KFitTrack.h:39
int getID() const
Unique id.
Definition: KHitBase.h:88
std::vector< double > fGTraceZMin
Graphical trace minimum z for each view.
const std::multimap< double, KHitTrack > & getTrackMap() const
KHitTrack collection, indexed by path distance.
Definition: KGTrack.h:54
std::optional< double > vec_prop(KTrack &trk, const std::shared_ptr< const Surface > &psurf, PropDirection dir, bool doDedx, TrackMatrix *prop_matrix=0, TrackError *noise_matrix=0) const
Propagate without error (long distance).
Definition: Propagator.cxx:52
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.
ublas::symmetric_matrix< double, ublas::lower, ublas::row_major, ublas::bounded_array< double, N *(N+1)/2 > > type
double fMaxPredDist
Maximum prediciton distance to accept a hit.
M::value_type trace(const M &m)
std::optional< double > noise_prop(KETrack &tre, const std::shared_ptr< const Surface > &psurf, PropDirection dir, bool doDedx, KTrack *ref=0) const
Propagate with error and noise.
Definition: Propagator.cxx:382
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:289
bool fitMomentumMS(const KGTrack &trg, const Propagator &prop, KETrack &tremom) const
Estimate track momentum using multiple scattering.
void hits()
Definition: readHits.C:15
void setChisq(double chisq)
Set chisquare.
Definition: KFitTrack.h:68
bool buildTrack(const KTrack &trk, KGTrack &trg, const Propagator &prop, const Propagator::PropDirection dir, KHitContainer &hits, bool linear) const
Make a new track.
TCanvas * c2
Definition: plot_hist.C:75
double getPath() const
Propagation distance.
Definition: KFitTrack.h:61
Kalman filter measurement class template.
ublas::vector< double, ublas::bounded_array< double, N > > type
bool syminvert(ublas::symmetric_matrix< T, TRI, L, A > &m)
double fGTraceWH
Window height.
TVirtualPad * fInfoPad
Information pad.
T get(std::string const &key) const
Definition: ParameterSet.h:314
bool extendTrack(KGTrack &trg, const Propagator &prop, KHitContainer &hits) const
Add hits to existing track.
std::optional< double > err_prop(KETrack &tre, const std::shared_ptr< const Surface > &psurf, PropDirection dir, bool doDedx, KTrack *ref=0, TrackMatrix *prop_matrix=0) const
Propagate with error, but without noise.
Definition: Propagator.cxx:343
std::vector< std::unique_ptr< TCanvas > > fCanvases
Graphical trace canvases.
const std::list< KHitGroup > & getUnsorted() const
Definition: KHitContainer.h:70
double fGTraceXMin
Graphical trace minimum x.
double fGoodPErr
Pointing error threshold for switching to free propagation.
TPaveText * fMessages
Message box.
void setStat(FitStatus stat)
Set fit status.
Definition: KFitTrack.h:69
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.
bool smoothTrack(KGTrack &trg, KGTrack *trg1, const Propagator &prop) const
Smooth track.
double fMinSampleDist
Minimum sample distance (for momentum measurement).
Kalman Filter.
Detector simulation of raw signals on wires.
double fMinSortDist
Sort low distance threshold.
const TrackVector & getVector() const
Track state vector.
Definition: KTrack.h:54
virtual void update(KETrack &tre) const =0
Update track method.
virtual double getChisq() const =0
Return incremental chisquare.
const KVector< N >::type & getMeasVector() const
Measurement vector.
Definition: KHit.h:101
TDirectory * dir
Definition: macro.C:5
constexpr double dist(const TReal *x, const TReal *y, const unsigned int dimension)
double fGTraceXMax
Graphical trace maximum x.
void sort(const KTrack &trk, bool addUnsorted, const Propagator &prop, Propagator::PropDirection dir=Propagator::UNKNOWN)
(Re)sort objects in unsorted and sorted lists.
std::vector< TVirtualPad * > fPads
View pads in current canvas.
TLegend * leg
Definition: compare.C:67
size_t numHits() const
Number of measurements in track.
Definition: KGTrack.h:57
double fGTraceWW
Window width.
ublas::matrix< double, ublas::row_major, ublas::bounded_array< double, N *M > > type
virtual bool predict(const KETrack &tre, const Propagator &prop, const KTrack *ref=0) const =0
Prediction method (return false if fail).
const std::shared_ptr< const Surface > & getSurface() const
Surface accessor.
Definition: KHitGroup.h:48
bool combineFit(const KFitTrack &trf)
Combine two tracks.
Definition: KFitTrack.cxx:65
void setPath(double path)
Set propagation distance.
Definition: KFitTrack.h:67
std::vector< double > fGTraceZMax
Graphical trace maximum z for each view.
int fMaxSamePlane
Maximum consecutive hits in same plane.
Line surface perpendicular to x-axis.
int fMaxNoiseHits
Maximum number of hits in noise cluster.
A collection of KHitTracks.
Basic Kalman filter track class, with fit information.
double fMaxPErr
Maximum pointing error for free propagation.
bool fitMomentumRange(const KGTrack &trg, KETrack &tremom) const
Estimate track momentum using range.
bool isValid() const
Validity flag.
Definition: KGTrack.h:66
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.
double getChisq() const
Fit chisquare.
Definition: KFitTrack.h:62
const KHitTrack & startTrack() const
Track at start point.
Definition: KGTrack.cxx:30
art framework interface to geometry description
const std::shared_ptr< const Surface > & getMeasSurface() const
Measurement surface.
Definition: KHitBase.h:82
int fMinLHits
Minimum number of hits to turn off linearized propagation.
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
decltype(auto) constexpr empty(T &&obj)
ADL-aware version of std::empty.
Definition: StdUtils.h:109
bool isValid() const
Test if track is valid.
Definition: KTrack.cxx:81
PropDirection
Propagation direction enum.
Definition: Propagator.h:94
FitStatus getStat() const
Fit status.
Definition: KFitTrack.h:63
int getPlane() const
Plane index.
Definition: KHitGroup.h:51
int fPlane
Preferred view plane.
const std::list< KHitGroup > & getSorted() const
Definition: KHitContainer.h:69