Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GblFitterInfo.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file 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 
78  void GblFitterInfo::setReferenceState(StateOnPlane& referenceState) {
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 
111  gbl::GblPoint GblFitterInfo::constructGblPoint() {
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 
195  void GblFitterInfo::updateMeasurementAndPlane(const StateOnPlane & sop) {
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 
246  void GblFitterInfo::updateFitResults(gbl::GblTrajectory& traj) {
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 
346  void GblFitterInfo::recalculateJacobian(GblFitterInfo* prevFitterInfo)
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
360  SharedPlanePtr oldPlane(new genfit::DetPlane(*sharedPlane_));
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 
453  MeasurementOnPlane GblFitterInfo::getKink() const {
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 
473  MeasurementOnPlane GblFitterInfo::getMeasurement() const{
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 
482  GblFitterInfo* GblFitterInfo::clone() const {
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