Belle II Software  release-05-01-25
GblFitterInfo.cc
1 #include "GblFitterInfo.h"
2 #include "MeasurementOnPlane.h"
3 #include "HMatrixUV.h"
4 #include "HMatrixV.h"
5 #include "HMatrixU.h"
6 
8 namespace genfit {
9 
10  GblFitterInfo::GblFitterInfo() : AbsFitterInfo(), jacobian_(TMatrixD(5, 5)) {
11  reset();
12  }
13 
14  GblFitterInfo::GblFitterInfo(const TrackPoint* trackPoint, const AbsTrackRep* rep) : AbsFitterInfo(trackPoint, rep), jacobian_(TMatrixD(5, 5)) {
15  reset();
16 
17  // Retrieve jacobian from rep by default (may not be used ... at 1st point)
18  unsigned int dim = rep->getDim();
19  if (dim != 5)
20  throw new genfit::Exception("GblFitterInfo: Representation state is not 5D", __LINE__, __FILE__);
21  TMatrixDSym noise(dim, dim);
22  TVectorD dState(dim);
23  rep->getForwardJacobianAndNoise(jacobian_, noise, dState);
24 
25  }
26 
27  GblFitterInfo::GblFitterInfo(const TrackPoint* trackPoint, const AbsTrackRep* rep, StateOnPlane& referenceState) : AbsFitterInfo(trackPoint, rep), jacobian_(TMatrixD(5, 5)) {
28  reset();
29 
30  // Retrieve jacobian from rep by default (may not be used ... at 1st point)
31  unsigned int dim = rep->getDim();
32  if (dim != 5)
33  throw new genfit::Exception("GblFitterInfo: Representation state is not 5D", __LINE__, __FILE__);
34  TMatrixDSym noise(dim, dim);
35  TVectorD dState(dim);
36  rep->getForwardJacobianAndNoise(jacobian_, noise, dState);
37 
38  setReferenceState(referenceState);
39  }
40 
41  void GblFitterInfo::reset(unsigned int measurementDim, unsigned int repDim) {
42  refPrediction_.ResizeTo(repDim);
43  refPrediction_.Zero();
44 
45  measResiduals_.ResizeTo(measurementDim);
46  measResiduals_.Zero();
47  kinkResiduals_.ResizeTo(measurementDim);
48  kinkResiduals_.Zero();
49  measResidualErrors_.ResizeTo(measurementDim);
50  measResidualErrors_.Zero();
51  kinkResidualErrors_.ResizeTo(measurementDim);
52  kinkResidualErrors_.Zero();
53  measDownWeights_.ResizeTo(measurementDim);
54  measDownWeights_.Zero();
55  kinkDownWeights_.ResizeTo(measurementDim);
56  kinkDownWeights_.Zero();
57 
58  TMatrixDSym zero(repDim);
59  zero.Zero();
60 
61  fwdStateCorrection_.ResizeTo(repDim);
62  fwdCov_.ResizeTo(zero);
63  bwdStateCorrection_.ResizeTo(repDim);
64  bwdCov_.ResizeTo(zero);
65  fwdPrediction_.ResizeTo(repDim);
66  bwdPrediction_.ResizeTo(repDim);
67  fwdCov_.SetMatrixArray(zero.GetMatrixArray());
68  bwdCov_.SetMatrixArray(zero.GetMatrixArray());
69 
70  measurement_.ResizeTo(2);
71  measurement_.Zero();
72  measCov_.ResizeTo(TMatrixDSym(2));
73  measCov_.Zero();
74  hMatrix_.ResizeTo(HMatrixUV().getMatrix());
75  hMatrix_ = HMatrixUV().getMatrix();
76  }
77 
79  updateMeasurementAndPlane(referenceState);
80 
81  refPrediction_ = referenceState.getState();
82  fwdPrediction_ = referenceState.getState();
83  bwdPrediction_ = referenceState.getState();
84 
85  //TODO: reset even if already fitted?
86  fittedStateBwd_.reset(new MeasuredStateOnPlane(referenceState.getState(), trackPoint_->getTrack()->getCovSeed(), sharedPlane_, rep_, referenceState.getAuxInfo()));
87  fittedStateFwd_.reset(new MeasuredStateOnPlane(referenceState.getState(), trackPoint_->getTrack()->getCovSeed(), sharedPlane_, rep_, referenceState.getAuxInfo()));
88 
89  }
90 
91  void GblFitterInfo::setJacobian(TMatrixD jacobian) {
92  jacobian_.ResizeTo(jacobian);
93  jacobian_ = jacobian;
94  }
95 
96  TMatrixDSym GblFitterInfo::getCovariance(double variance, TVector3 trackDirection, SharedPlanePtr measurementPlane) const {
97 
98  double c1 = trackDirection.Dot(measurementPlane->getU());
99  double c2 = trackDirection.Dot(measurementPlane->getV());
100 
101  TMatrixDSym scatCov(2);
102  scatCov(0, 0) = 1. - c2 * c2;
103  scatCov(1, 1) = 1. - c1 * c1;
104  scatCov(0, 1) = c1 * c2;
105  scatCov(1, 0) = c1 * c2;
106  scatCov *= variance * variance / (1. - c1 * c1 - c2 * c2) / (1. - c1 * c1 - c2 * c2) ;
107 
108  return scatCov;
109  }
110 
112  // All meas/scat info is added from genfit data again (to cope with possible RecoHit update)
113 
114  // Create GBL point with current jacobian
115  gbl::GblPoint thePoint(jacobian_);
116 
117  //NOTE: 3rd update and update anytime GblPoint is requested
118  // mostly likely will update with reference as on 2nd update
119  StateOnPlane sop(getFittedState().getState(), sharedPlane_, rep_, getFittedState(true).getAuxInfo());
120 
121 
122  // Scatterer
123  if (trackPoint_->hasThinScatterer()) {
124  if (!hasMeasurements()) {
125  //double scatVariance = trackPoint_->getMaterialInfo()->getMaterial().getDensity();
126  //TVectorD kinkPrec(2);
127  //kinkPrec(0) = 1./scatVariance/scatVariance; kinkPrec(1) = 1./scatVariance/scatVariance;
128  //thePoint.addScatterer(getKinks(), kinkPrec);
129  //TODO: if state at scatterer is updated, the direction of track might not be perpendicular anymore
130  // plane does not change at pure scatterer
131  TMatrixDSym kinkCov = getCovariance(trackPoint_->getMaterialInfo()->getMaterial().density, sop.getDir(), sop.getPlane());
132  thePoint.addScatterer(getKinks(), kinkCov.Invert());
133  } else {
134  TMatrixDSym kinkCov = getCovariance(trackPoint_->getMaterialInfo()->getMaterial().density, sop.getDir(), sop.getPlane());
135  thePoint.addScatterer(getKinks(), kinkCov.Invert());
136  }
137  }
138 
139 
140 
141  // Measurement
142  if (hasMeasurements()) {
143  MeasuredStateOnPlane measurement = getResidual(0, true, true);
144  TVectorD aResiduals(measurement.getState());
145  TMatrixDSym aPrecision(measurement.getCov().Invert());
146  if (HMatrixU().getMatrix() == hMatrix_) {
147  double res = aResiduals(0);
148  double prec = aPrecision(0, 0);
149  aResiduals.ResizeTo(2);
150  aPrecision.ResizeTo(TMatrixDSym(2));
151  aResiduals.Zero();
152  aResiduals(0) = res;
153  aPrecision.Zero();
154  aPrecision(0, 0) = prec;
155  }
156  if (HMatrixV().getMatrix() == hMatrix_) {
157  double res = aResiduals(0);
158  double prec = aPrecision(0, 0);
159  aResiduals.ResizeTo(2);
160  aPrecision.ResizeTo(TMatrixDSym(2));
161  aResiduals.Zero();
162  aResiduals(1) = res;
163  aPrecision.Zero();
164  aPrecision(1, 1) = prec;
165  }
166  // always 2D, U/V set by precision (=0 to disable)
167  thePoint.addMeasurement(aResiduals, aPrecision);
168  }
169 
170  // Derivatives
171  ICalibrationParametersDerivatives* globals = nullptr;
172  if (hasMeasurements() && (globals = dynamic_cast<ICalibrationParametersDerivatives*>(trackPoint_->getRawMeasurement(0)) )) {
173  std::pair<std::vector<int>, TMatrixD> labelsAndMatrix = globals->globalDerivatives(&sop);
174  std::vector<int> labels = labelsAndMatrix.first;
175  TMatrixD derivs = labelsAndMatrix.second;
176 
177  if (derivs.GetNcols() > 0 && !labels.empty() && (unsigned int)derivs.GetNcols() == labels.size()) {
178  thePoint.addGlobals(labels, derivs);
179  }
180  TMatrixD locals = globals->localDerivatives(&sop);
181  if (locals.GetNcols() > 0) {
182  thePoint.addLocals(locals);
183  GblFitStatus* gblfs = dynamic_cast<GblFitStatus*>(trackPoint_->getTrack()->getFitStatus(rep_));
184  if (gblfs) {
185  if (gblfs->getMaxLocalFitParams() < locals.GetNcols())
186  gblfs->setMaxLocalFitParams(locals.GetNcols());
187  }
188  }
189  }
190 
191  return thePoint;
192 
193  }
194 
196  if (!trackPoint_)
197  return;
198  if (!trackPoint_->hasRawMeasurements()) {
199  //no measurement, plane is supposed to come from state (is perpendicular)
200  setPlane(sop.getPlane());
201  return;
202  }
203  std::vector<MeasurementOnPlane*> allMeas = trackPoint_->getRawMeasurement(0)->constructMeasurementsOnPlane(sop);
204 
205  /*
206  double normMin(9.99E99);
207  unsigned int imop(0);
208  const AbsHMatrix* H = allMeas[0]->getHMatrix();
209  for (unsigned int i=0; i<allMeas.size(); ++i) {
210  if (*(allMeas[i]->getHMatrix()) != *H){
211  Exception e("GblFitterInfo::updateMeasurementAndPlane: Cannot compare measurements with different H-Matrices.", __LINE__,__FILE__);
212  e.setFatal();
213  throw e;
214  }
215 
216  TVectorD res = allMeas[i]->getState() - H->Hv(sop.getState());
217  double norm = sqrt(res.Norm2Sqr());
218  if (norm < normMin) {
219  normMin = norm;
220  imop = i;
221  }
222  }
223  */
224  unsigned int imop = 0;
225  double maxWeight = allMeas.at(0)->getWeight();
226  for (unsigned int i = 0; i < allMeas.size(); i++)
227  if (allMeas.at(i)->getWeight() > maxWeight)
228  imop = i;
229 
230  measurement_.ResizeTo(allMeas.at(imop)->getState());
231  measurement_ = allMeas.at(imop)->getState();
232  measCov_.ResizeTo(allMeas.at(imop)->getCov());
233  measCov_ = allMeas.at(imop)->getCov();
234  hMatrix_.ResizeTo(allMeas.at(imop)->getHMatrix()->getMatrix());
235  hMatrix_ = allMeas.at(imop)->getHMatrix()->getMatrix();
236 
237  setPlane(allMeas.at(imop)->getPlane());
238 
239  for (unsigned int imeas = 0; imeas < allMeas.size(); imeas++)
240  delete allMeas[imeas];
241  allMeas.clear();
242 
243 
244  }
245 
247  if (!traj.isValid())
248  return;
249 
250  // Deduce our position in the trajectory
251  //----------------------------------------------
252  unsigned int label = 0;
253  // Loop over track and find index of this fitter info
254  genfit::Track* trk = trackPoint_->getTrack();
255  for (unsigned int ip = 0; ip < trk->getNumPoints(); ip++) {
256  // Indexing of GblFitter info in track (each gives one point to trajectory at trackpoint position)
257  if (dynamic_cast<GblFitterInfo*>( trk->getPoint(ip)->getFitterInfo(rep_) )) {
258  // First fitter info has label 1 (for point in GBL traj)
259  label++;
260  // We found itself = we have our label
261  if (trk->getPoint(ip)->getFitterInfo(rep_) == this)
262  break;
263  }
264  }
265  if (label == 0)
266  throw genfit::Exception("GblFitterInfo: fitter info did not found itself in track to update", __LINE__, __FILE__);
267  if (label > traj.getNumPoints())
268  throw genfit::Exception("GblFitterInfo: Deduced point label not valid", __LINE__, __FILE__);
269  //----------------------------------------------
270 
271  // Residuals (in plane)
272  unsigned int numMRes = 2;
273  TVectorD mResiduals(2), mMeasErrors(2), mResErrors(2), mDownWeights(2);
274  if (0 != traj.getMeasResults(label, numMRes, mResiduals, mMeasErrors, mResErrors, mDownWeights))
275  throw genfit::Exception(" NO measurement results ", __LINE__, __FILE__);
276 
277  // Kinks
278  unsigned int numKRes = 2;
279  TVectorD kResiduals(2), kMeasErrors(2), kResErrors(2), kDownWeights(2);
280  if (0 != traj.getScatResults(label, numKRes, kResiduals, kMeasErrors, kResErrors, kDownWeights))
281  throw genfit::Exception(" NO scattering results ", __LINE__, __FILE__);
282 
283  // Check for local derivatives
284  int nLocals = 0;
285  GblFitStatus* gblfs = dynamic_cast<GblFitStatus*>(trackPoint_->getTrack()->getFitStatus(rep_));
286  if (gblfs)
287  nLocals = gblfs->getMaxLocalFitParams();
288 
289  // Predictions (different at scatterer)
290  TVectorD bwdUpdate(5 + nLocals), fwdUpdate(5 + nLocals);
291  TMatrixDSym bwdCov(5 + nLocals), fwdCov(5 + nLocals);
292 
293  // forward prediction
294  if (0 != traj.getResults(label, fwdUpdate, fwdCov))
295  throw genfit::Exception(" NO forward results ", __LINE__, __FILE__);
296 
297  // backward prediction
298  if (0 != traj.getResults(-1 * label, bwdUpdate, bwdCov))
299  throw genfit::Exception(" NO backward results ", __LINE__, __FILE__);
300 
301  if (nLocals > 0) {
302  TVectorD _bwdUpdate(5 + nLocals), _fwdUpdate(5 + nLocals);
303  TMatrixDSym _bwdCov(5 + nLocals), _fwdCov(5 + nLocals);
304  _bwdUpdate = bwdUpdate;
305  _fwdUpdate = fwdUpdate;
306  _bwdCov = bwdCov;
307  _fwdCov = fwdCov;
308  bwdUpdate.ResizeTo(5);
309  fwdUpdate.ResizeTo(5);
310  bwdCov.ResizeTo(TMatrixDSym(5));
311  fwdCov.ResizeTo(TMatrixDSym(5));
312  for (int i = 0; i < 5; i++) {
313  bwdUpdate(i) = _bwdUpdate(i);
314  fwdUpdate(i) = _fwdUpdate(i);
315  for (int j = 0; j < 5; j++) {
316  bwdCov(i, j) = _bwdCov(i, j);
317  fwdCov(i, j) = _fwdCov(i, j);
318  }
319  }
320  }
321  // Now update the the fitter info
322  //
323  //-------------------------------------------------
324  // Backward/forward prediction (residual) (differs at scatterers) AFTER GBL fit
325  bwdStateCorrection_ = bwdUpdate;
326  fwdStateCorrection_ = fwdUpdate;
327  bwdCov_ = bwdCov;
328  fwdCov_ = fwdCov;
329  fwdPrediction_ += fwdStateCorrection_; // This is the update!
330  bwdPrediction_ += bwdStateCorrection_; // This is the update!
331 
332  fittedStateFwd_.reset( new MeasuredStateOnPlane(fwdPrediction_, fwdCov_, sharedPlane_, rep_, getFittedState(true).getAuxInfo()) );
333  fittedStateBwd_.reset( new MeasuredStateOnPlane(bwdPrediction_, bwdCov_, sharedPlane_, rep_, getFittedState(true).getAuxInfo()) );
334 
335  // Set scattering/measurement residual data
336  kinkResiduals_ = kResiduals;
337  measResiduals_ = mResiduals;
338  kinkResidualErrors_ = kResErrors;
339  measResidualErrors_ = mResErrors;
340  measDownWeights_ = mDownWeights;
341  kinkDownWeights_ = kDownWeights;
342 
343  //-------------------------------------------------
344  }
345 
347  {
348  // Invalidates errors and corrections from last iteration
349  // (will be defined in different plane). But fitted state and residual is ok.
350 
351  if (!prevFitterInfo) {
352  jacobian_.UnitMatrix();
353 
354  //TODO: For 1st plane this has no sense
355  //return;
356  prevFitterInfo = this;
357  }
358 
359  //TODO
361 
362  //updateMeasurementAndPlane(StateOnPlane(fwdPrediction_, sharedPlane_, rep_));
363  //
364 
365  TMatrixDSym noise;
366  TVectorD dstate;
367  // Take forward state from previous fitter info,
368  // its (maybe updated) plane
369  // and our rep
370  StateOnPlane prevState(prevFitterInfo->getFittedState(true).getState(), prevFitterInfo->getPlane(), rep_, getFittedState(true).getAuxInfo());
371 
372  if (hasMeasurements()) {
373  SharedPlanePtr newPlane = trackPoint_->getRawMeasurement(0)->constructPlane(prevState);
374  rep_->extrapolateToPlane(prevState, newPlane, false, true);
375  } else {
376  rep_->extrapolateToPlane(prevState, sharedPlane_, false, true);
377  }
378 
379  rep_->getForwardJacobianAndNoise(jacobian_, noise, dstate);
380  // Now update meas data
381  updateMeasurementAndPlane(prevState);
382 
383  //
384  // Extrap predictions to new plane
385  //
386  StateOnPlane oldFwdState(fwdPrediction_, oldPlane, rep_, getFittedState(true).getAuxInfo());
387  StateOnPlane oldBwdState(bwdPrediction_, oldPlane, rep_, getFittedState(true).getAuxInfo());
388  rep_->extrapolateToPlane(oldFwdState, sharedPlane_);
389  rep_->extrapolateToPlane(oldBwdState, sharedPlane_);
390  fwdPrediction_ = oldFwdState.getState();
391  bwdPrediction_ = oldBwdState.getState();
392  fittedStateBwd_.reset(new MeasuredStateOnPlane(fwdPrediction_, fwdCov_, sharedPlane_, rep_, getFittedState(true).getAuxInfo()));
393  fittedStateFwd_.reset(new MeasuredStateOnPlane(bwdPrediction_, bwdCov_, sharedPlane_, rep_, getFittedState(true).getAuxInfo()));
394  //
395  }
396 
397 
398  const MeasuredStateOnPlane& GblFitterInfo::getFittedState(bool afterKink) const {
399  // ALways biased from GBL (global fit!)
400 
401  if (!fittedStateFwd_ || !fittedStateBwd_) {
402  //NOTE: This should be already set (from reference)! The auxInfo is being book-kept by it. If reference is not set, default auxInfo is used
403  fittedStateFwd_.reset(new MeasuredStateOnPlane(fwdPrediction_, fwdCov_, sharedPlane_, rep_));
404  fittedStateBwd_.reset(new MeasuredStateOnPlane(bwdPrediction_, bwdCov_, sharedPlane_, rep_));
405  }
406 
407  if (afterKink) {
408  return *fittedStateFwd_;
409  }
410  else {
411  return *fittedStateBwd_;
412  }
413 
414  }
415 
416  MeasurementOnPlane GblFitterInfo::getResidual(unsigned int, bool, bool onlyMeasurementErrors) const {
417  // Return error of resduals (=0, none, for reference track, no errors known before fit (seed covariance might not be correct)
418  TMatrixDSym localCovariance(2);
419  localCovariance.Zero();
420  // combined meas+fit errors:
421  //TODO: 1D covariance + more dimensions of residuals
422 
423  localCovariance(0, 0) = measResidualErrors_(0) * measResidualErrors_(0);
424  localCovariance(1, 1) = measResidualErrors_(1) * measResidualErrors_(1);
425 
426  if (HMatrixU().getMatrix() == hMatrix_) {
427  localCovariance.ResizeTo(TMatrixDSym(1));
428  localCovariance(0, 0) = measResidualErrors_(0) * measResidualErrors_(0);
429  }
430  if (HMatrixV().getMatrix() == hMatrix_) {
431  localCovariance.ResizeTo(TMatrixDSym(1));
432  localCovariance(0, 0) = measResidualErrors_(1) * measResidualErrors_(1);
433  }
434 
435  if (hasMeasurements()){
436  // this means only for reference state before gbl fit, this way will be used
437  TVectorD res( measurement_ - hMatrix_ * getFittedState(true).getState() );
438 
439  if (onlyMeasurementErrors) {
440  localCovariance.ResizeTo(measCov_);
441  localCovariance = measCov_;
442  }
443  return MeasurementOnPlane(res, localCovariance, sharedPlane_, rep_, trackPoint_->getRawMeasurement(0)->constructHMatrix(getRep()));
444  }
445  TVectorD zeroRes(2);
446  zeroRes.Zero();
447  // Else return 0's or whatever
448  //TODO: or throw?
449  return MeasurementOnPlane(zeroRes, localCovariance, sharedPlane_, rep_, new HMatrixUV());
450 
451  } // calculate residual, track and measurement errors are added if onlyMeasurementErrors is false
452 
454  TMatrixDSym localCovariance(2);
455  localCovariance.Zero();
456  localCovariance(0, 0) = kinkResidualErrors_(0) * kinkResidualErrors_(0);
457  localCovariance(1, 1) = kinkResidualErrors_(1) * kinkResidualErrors_(1);
458  return MeasurementOnPlane(getKinks(), localCovariance, sharedPlane_, rep_, new genfit::HMatrixUV());
459 
460  }
461 
462  TVectorD GblFitterInfo::getKinks() const {
463  TVectorD kinks(2);
464  kinks.Zero();
465 
466  TVectorD stateDiff(getFittedState(true).getState() - getFittedState(false).getState());
467  kinks(0) = -stateDiff(1);
468  kinks(1) = -stateDiff(2);
469 
470  return kinks;
471  }
472 
474  return MeasurementOnPlane(measurement_, measCov_, sharedPlane_, rep_, hasMeasurements() ? trackPoint_->getRawMeasurement(0)->constructHMatrix(rep_) : new HMatrixUV() );
475  }
476 
477  bool GblFitterInfo::checkConsistency(const genfit::PruneFlags*) const {
478  //TODO
479  return true;
480  }
481 
483 
484  GblFitterInfo* retVal = new GblFitterInfo(this->getTrackPoint(), this->getRep());
485 
486  retVal->setPlane(sharedPlane_);
487 
488  retVal->jacobian_ = jacobian_;
489  retVal->measResiduals_ = measResiduals_;
490  retVal->measResidualErrors_ = measResidualErrors_;
491  retVal->kinkResiduals_ = kinkResiduals_;
492  retVal->kinkResidualErrors_ = kinkResidualErrors_;
493  retVal->measDownWeights_ = measDownWeights_;
494  retVal->kinkDownWeights_ = kinkDownWeights_;
495  retVal->bwdStateCorrection_ = bwdStateCorrection_;
496  retVal->fwdStateCorrection_ = fwdStateCorrection_;
497  retVal->bwdCov_ = bwdCov_;
498  retVal->fwdCov_ = fwdCov_;
499  retVal->fwdPrediction_ = fwdPrediction_;
500  retVal->bwdPrediction_ = bwdPrediction_;
501  retVal->refPrediction_ = refPrediction_;
502  retVal->measurement_.ResizeTo(measurement_);
503  retVal->measCov_.ResizeTo(measCov_);
504  retVal->hMatrix_.ResizeTo(hMatrix_);
505  retVal->measurement_ = measurement_;
506  retVal->measCov_ = measCov_;
507  retVal->hMatrix_ = hMatrix_;
508 
509  return retVal;
510  }
511 
512  void GblFitterInfo::Print(const Option_t*) const {
513  //TODO
514  std::cout << "=============================================================================================" << std::endl;
515  std::cout << " >>> GblFitterInfo " << std::endl;
516  std::cout << " ************* " << std::endl;
517 
518  std::cout << " rep: " << rep_ << ", trackpoint: " << trackPoint_ << ", plane: " << sharedPlane_.get() << std::endl;
519  sharedPlane_->Print();
520  std::cout << std::endl;
521 
522  std::cout << "=============================================================================================" << std::endl;
523  std::cout << " | PREDICTIONS | REFERENCE | Corrections from last iteration |" << std::endl;
524  std::cout << " | (+)prediction | (-)prediction | state | (+)correction | (-) correction |" << std::endl;
525  std::cout << "---------------------------------------------------------------------------------------------" << std::endl;
526 
527  for (int i = 0; i <5; i++) {
528  std::cout << std::left;
529  std::cout << " ";
530  if (i==0)
531  std::cout << "q/p";
532  if (i==1)
533  std::cout << "u' ";
534  if (i==2)
535  std::cout << "v' ";
536  if (i==3)
537  std::cout << "u ";
538  if (i==4)
539  std::cout << "v ";
540  std::cout << " | "
541  << std::setw(12) << fwdPrediction_(i) << " | "
542  << std::setw(12) << bwdPrediction_(i) << " | "
543  << std::setw(12) << refPrediction_(i) << " | "
544  << std::setw(12) << fwdStateCorrection_(i) << " | "
545  << std::setw(12) << bwdStateCorrection_(i) << std::endl;
546  }
547  std::cout << "=============================================================================================" << std::endl;
548 
549  if (hasMeasurements()) {
550  std::cout << " | Meas. residual | measurement - prediction | Down-weight | Fit+meas Err. |" << std::endl;
551  std::cout << " | | | | -diagonaliz. |" << std::endl;
552  std::cout << "---------------------------------------------------------------------------------------------" << std::endl;
553 
554  TVectorD residual = getResidual().getState();
555  if (residual.GetNoElements()<2) {
556  residual.ResizeTo(2);
557  residual.Zero();
558  if (trackPoint_->getRawMeasurement(0)->constructHMatrix(getRep())->isEqual(genfit::HMatrixU()))
559  residual(0) = getResidual().getState()(0);
560  else
561  residual(1) = getResidual().getState()(0);
562 
563  }
564  std::cout << " u | "
565  << std::setw(12) << measResiduals_(0) << " | "
566  << std::setw(12) << residual(0) << " | "
567  << std::setw(12) << measDownWeights_(0) << " | "
568  << std::setw(12) << measResidualErrors_(0)
569  << std::endl;
570 
571  std::cout << " v | "
572  << std::setw(12) << measResiduals_(1) << " | "
573  << std::setw(12) << residual(1) << " | "
574  << std::setw(12) << measDownWeights_(1) << " | "
575  << std::setw(12) << measResidualErrors_(1)
576  << std::endl;
577 
578  std::cout << "---------------------------------------------------------------------------------------------" << std::endl;
579  }
580 
581  std::cout << " | Kink residual | Residual of slope difference | Down-weight | Fit Kink Err. |" << std::endl;
582  std::cout << " | -diagonalized | - ( (+)pred - (-)pred ) | | -diagonaliz. |" << std::endl;
583  std::cout << "---------------------------------------------------------------------------------------------" << std::endl;
584 
585  std::cout << " u' | "
586  << std::setw(12) << kinkResiduals_(0) << " | "
587  << std::setw(12) << getKinks()(0) << " | "
588  << std::setw(12) << kinkDownWeights_(0) << " | "
589  << std::setw(12) << kinkResidualErrors_(0)
590  << std::endl;
591 
592  std::cout << " v' | "
593  << std::setw(12) << kinkResiduals_(1) << " | "
594  << std::setw(12) << getKinks()(1) << " | "
595  << std::setw(12) << kinkDownWeights_(1) << " | "
596  << std::setw(12) << kinkResidualErrors_(1)
597  << std::endl;
598  std::cout << "=============================================================================================" << std::endl;
599  std::cout << "Measurement: "; measurement_.Print();
600 
601  std::cout << "H Matrix: "; hMatrix_.Print();
602 
603  std::cout << "Measurement covariance: "; measCov_.Print();
604 
605  std::cout << "Jacobian: "; jacobian_.Print();
606  std::cout << "Backward covariance: "; bwdCov_.Print();
607  std::cout << "Forward covariance : "; fwdCov_.Print();
608 
609  std::cout << "=============================================================================================" << std::endl;
610 
611  }
612 
613 
614 } // end of namespace genfit
gbl::GblTrajectory::getMeasResults
unsigned int getMeasResults(unsigned int aLabel, unsigned int &numRes, TVectorD &aResiduals, TVectorD &aMeasErrors, TVectorD &aResErrors, TVectorD &aDownWeights)
Get residuals from fit at point for measurement.
Definition: GblTrajectory.cc:661
genfit::AbsMeasurement::constructMeasurementsOnPlane
virtual std::vector< genfit::MeasurementOnPlane * > constructMeasurementsOnPlane(const StateOnPlane &state) const =0
Construct MeasurementOnPlane on plane of the state and wrt the states TrackRep.
gbl::GblTrajectory::getResults
unsigned int getResults(int aSignedLabel, TVectorD &localPar, TMatrixDSym &localCov) const
Get fit results at point.
Definition: GblTrajectory.cc:631
genfit::Exception
Exception class for error handling in GENFIT (provides storage for diagnostic information)
Definition: Exception.h:48
genfit::TrackPoint
Object containing AbsMeasurement and AbsFitterInfo objects.
Definition: TrackPoint.h:46
genfit::SharedPlanePtr
std::shared_ptr< genfit::DetPlane > SharedPlanePtr
Shared Pointer to a DetPlane.
Definition: SharedPlanePtr.h:40
genfit::AbsFitterInfo::sharedPlane_
SharedPlanePtr sharedPlane_
No ownership.
Definition: AbsFitterInfo.h:94
genfit::GblFitterInfo::getResidual
MeasurementOnPlane getResidual(unsigned int=0, bool=false, bool onlyMeasurementErrors=true) const override
Get the residual.
Definition: GblFitterInfo.cc:416
genfit::MeasuredStateOnPlane
#StateOnPlane with additional covariance matrix.
Definition: MeasuredStateOnPlane.h:39
genfit::GblFitStatus
FitStatus for use with GblFitter.
Definition: GblFitStatus.h:39
genfit::PruneFlags
Info which information has been pruned from the Track.
Definition: FitStatus.h:47
genfit::AbsTrackRep::extrapolateToPlane
virtual double extrapolateToPlane(StateOnPlane &state, const genfit::SharedPlanePtr &plane, bool stopAtBoundary=false, bool calcJacobianNoise=false) const =0
Extrapolates the state to plane, and returns the extrapolation length and, via reference,...
genfit::AbsFitterInfo
This class collects all information needed and produced by a specific AbsFitter and is specific to on...
Definition: AbsFitterInfo.h:42
genfit::AbsMeasurement::constructPlane
virtual SharedPlanePtr constructPlane(const StateOnPlane &state) const =0
Construct (virtual) detector plane (use state's AbsTrackRep).
genfit::Track::getFitStatus
FitStatus * getFitStatus(const AbsTrackRep *rep=nullptr) const
Get FitStatus for a AbsTrackRep. Per default, return FitStatus for cardinalRep.
Definition: Track.h:154
genfit::StateOnPlane
A state with arbitrary dimension defined in a DetPlane.
Definition: StateOnPlane.h:47
genfit
Defines for I/O streams used for error and debug printing.
Definition: AlignablePXDRecoHit.h:19
genfit::HMatrixUV
AbsHMatrix implementation for two-dimensional MeasurementOnPlane and RKTrackRep parameterization.
Definition: HMatrixUV.h:39
genfit::Track
Collection of TrackPoint objects, AbsTrackRep objects and FitStatus objects.
Definition: Track.h:71
gbl::GblPoint::addScatterer
void addScatterer(const TVectorD &aResiduals, const TVectorD &aPrecision)
Add a (thin) scatterer to a point.
Definition: GblPoint.cc:210
gbl::GblPoint::addGlobals
void addGlobals(const std::vector< int > &aLabels, const TMatrixD &aDerivatives)
Add global derivatives to a point.
Definition: GblPoint.cc:320
genfit::AbsTrackRep
Abstract base class for a track representation.
Definition: AbsTrackRep.h:66
gbl::GblTrajectory::getScatResults
unsigned int getScatResults(unsigned int aLabel, unsigned int &numRes, TVectorD &aResiduals, TVectorD &aMeasErrors, TVectorD &aResErrors, TVectorD &aDownWeights)
Get (kink) residuals from fit at point for scatterer.
Definition: GblTrajectory.cc:690
genfit::GblFitterInfo::GblFitterInfo
GblFitterInfo()
Constructor for ROOT I/O.
Definition: GblFitterInfo.cc:10
genfit::GblFitterInfo::getFittedState
const MeasuredStateOnPlane & getFittedState(bool afterKink=true) const override
Get the prediction at this point Always biased in GBL (global fit) There are 2 states,...
Definition: GblFitterInfo.cc:398
gbl::GblTrajectory
GBL trajectory.
Definition: GblTrajectory.h:48
genfit::HMatrixV
AbsHMatrix implementation for one-dimensional MeasurementOnPlane and RKTrackRep parameterization.
Definition: HMatrixV.h:37
genfit::GblFitterInfo::getKink
MeasurementOnPlane getKink() const
Get kink (residual) with diagonalized covariance (2D) Covariance may be zero if not yet fitted or no ...
Definition: GblFitterInfo.cc:453
genfit::GblFitterInfo::updateFitResults
void updateFitResults(gbl::GblTrajectory &traj)
Update fitter info from GBL fit results.
Definition: GblFitterInfo.cc:246
genfit::GblFitterInfo::setJacobian
void setJacobian(TMatrixD jacobian)
Set the Jacobian for further GblPoint construction.
Definition: GblFitterInfo.cc:91
genfit::GblFitterInfo::reset
void reset(unsigned int measurementDim=2, unsigned int repDim=5)
(Initial) reset of fitter info
Definition: GblFitterInfo.cc:41
genfit::GblFitterInfo::updateMeasurementAndPlane
void updateMeasurementAndPlane(const StateOnPlane &sop)
SHOULD BE USED ONLY INTERNALY! Update the plane from measurement constructed with state or take plane...
Definition: GblFitterInfo.cc:195
gbl::GblPoint::addLocals
void addLocals(const TMatrixD &aDerivatives)
Add local derivatives to a point.
Definition: GblPoint.cc:293
genfit::TrackPoint::getFitterInfo
AbsFitterInfo * getFitterInfo(const AbsTrackRep *rep=nullptr) const
Get fitterInfo for rep. Per default, use cardinal rep.
Definition: TrackPoint.cc:170
genfit::ICalibrationParametersDerivatives::globalDerivatives
virtual std::pair< std::vector< int >, TMatrixD > globalDerivatives(const genfit::StateOnPlane *sop)
Labels and derivatives of residuals (local measurement coordinates) w.r.t.
Definition: ICalibrationParametersDerivatives.h:80
genfit::GblFitterInfo::getCovariance
TMatrixDSym getCovariance(double variance, TVector3 trackDirection, SharedPlanePtr measurementPlane) const
Get scattering covariance projected into (measurement) plane.
Definition: GblFitterInfo.cc:96
genfit::ICalibrationParametersDerivatives::localDerivatives
virtual TMatrixD localDerivatives(const genfit::StateOnPlane *)
Derivatives for additional local parameters to be fitted in global calibration algorithms together wi...
Definition: ICalibrationParametersDerivatives.h:141
genfit::GblFitterInfo::recalculateJacobian
void recalculateJacobian(GblFitterInfo *prevFitterInfo)
Re-extrapolates between prevFitterInfo and this point using forward state to update the Jacobian (if ...
Definition: GblFitterInfo.cc:346
genfit::GblFitterInfo
Collects information needed and produced by a GblFitter/GBL and is specific to one AbsTrackRep of the...
Definition: GblFitterInfo.h:52
genfit::GblFitterInfo::setReferenceState
void setReferenceState(StateOnPlane &referenceState)
Set the prediction and plane (from measurement if any) You should use the user constructor instead.
Definition: GblFitterInfo.cc:78
genfit::ICalibrationParametersDerivatives
Abstract base class to establish an interface between physical representation of the detector for ali...
Definition: ICalibrationParametersDerivatives.h:44
genfit::HMatrixU
AbsHMatrix implementation for one-dimensional MeasurementOnPlane and RKTrackRep parameterization.
Definition: HMatrixU.h:37
genfit::GblFitterInfo::clone
virtual GblFitterInfo * clone() const override
Deep copy ctor for polymorphic class.
Definition: GblFitterInfo.cc:482
genfit::GblFitterInfo::getMeasurement
MeasurementOnPlane getMeasurement() const
Get the measurement on plane from stored measurement data (from last construction/update)
Definition: GblFitterInfo.cc:473
genfit::AbsTrackRep::getDim
virtual unsigned int getDim() const =0
Get the dimension of the state vector used by the track representation.
gbl::GblTrajectory::getNumPoints
unsigned int getNumPoints() const
Retrieve number of point from trajectory.
Definition: GblTrajectory.cc:294
genfit::AbsMeasurement::constructHMatrix
virtual const AbsHMatrix * constructHMatrix(const AbsTrackRep *) const =0
Returns a new AbsHMatrix object.
genfit::GblFitterInfo::constructGblPoint
gbl::GblPoint constructGblPoint()
Collect all data and create a GblPoint.
Definition: GblFitterInfo.cc:111
genfit::GblFitterInfo::getKinks
TVectorD getKinks() const
Get kink (residual) (2D) = 0 - ( (+)pred - (-)pred )
Definition: GblFitterInfo.cc:462
genfit::GblFitterInfo::fittedStateFwd_
std::unique_ptr< MeasuredStateOnPlane > fittedStateFwd_
cache
Definition: GblFitterInfo.h:268
genfit::DetPlane
Detector plane.
Definition: DetPlane.h:59
gbl::GblTrajectory::isValid
bool isValid() const
Retrieve validity of trajectory.
Definition: GblTrajectory.cc:289
genfit::AbsFitterInfo::trackPoint_
const TrackPoint * trackPoint_
Pointer to TrackPoint where the FitterInfo belongs to.
Definition: AbsFitterInfo.h:88
gbl::GblPoint::addMeasurement
void addMeasurement(const TMatrixD &aProjection, const TVectorD &aResiduals, const TVectorD &aPrecision, double minPrecision=0.)
Add a measurement to a point.
Definition: GblPoint.cc:69
gbl::GblPoint
Point on trajectory.
Definition: GblPoint.h:68
genfit::AbsTrackRep::getForwardJacobianAndNoise
virtual void getForwardJacobianAndNoise(TMatrixD &jacobian, TMatrixDSym &noise, TVectorD &deltaState) const =0
Get the jacobian and noise matrix of the last extrapolation.
genfit::MeasurementOnPlane
Measured coordinates on a plane.
Definition: MeasurementOnPlane.h:46
genfit::HMatrixUV::getMatrix
const TMatrixD & getMatrix() const override
Get the actual matrix representation.
Definition: HMatrixUV.cc:33
genfit::AbsFitterInfo::rep_
const AbsTrackRep * rep_
No ownership.
Definition: AbsFitterInfo.h:92