GBL/src/GblFitterInfo.cc

0001 #include "GblFitterInfo.h"
0002 #include "MeasurementOnPlane.h"
0003 #include "HMatrixUV.h"
0004 #include "HMatrixV.h"
0005 #include "HMatrixU.h"
0006
0007 //! Namespace for the general broken lines package
0008 namespace genfit {
0009
0010   GblFitterInfo::GblFitterInfo() : AbsFitterInfo(), jacobian_(TMatrixD(5, 5)) {
0011     reset();
0012   }
0013
0014   GblFitterInfo::GblFitterInfo(const TrackPoint* trackPoint, const AbsTrackRep* rep) : AbsFitterInfo(trackPoint, rep), jacobian_(TMatrixD(5, 5)) {
0015     reset();
0016
0017     // Retrieve jacobian from rep by default (may not be used ... at 1st point)
0018     unsigned int dim = rep->getDim();
0019     if (dim != 5)
0020       throw new genfit::Exception("GblFitterInfo: Representation state is not 5D", __LINE__, __FILE__);
0021     TMatrixDSym noise(dim, dim);
0022     TVectorD dState(dim);
0023     rep->getForwardJacobianAndNoise(jacobian_, noise, dState);
0024
0025   }
0026
0027   GblFitterInfo::GblFitterInfo(const TrackPoint* trackPoint, const AbsTrackRep* rep, StateOnPlane& referenceState) : AbsFitterInfo(trackPoint, rep), jacobian_(TMatrixD(5, 5)) {
0028     reset();
0029
0030     // Retrieve jacobian from rep by default (may not be used ... at 1st point)
0031     unsigned int dim = rep->getDim();
0032     if (dim != 5)
0033       throw new genfit::Exception("GblFitterInfo: Representation state is not 5D", __LINE__, __FILE__);
0034     TMatrixDSym noise(dim, dim);
0035     TVectorD dState(dim);
0036     rep->getForwardJacobianAndNoise(jacobian_, noise, dState);
0037
0038     setReferenceState(referenceState);
0039   }
0040
0041   void GblFitterInfo::reset(unsigned int measurementDim, unsigned int repDim) {
0042     refPrediction_.ResizeTo(repDim);
0043     refPrediction_.Zero();
0044
0045     measResiduals_.ResizeTo(measurementDim);
0046     measResiduals_.Zero();
0047     kinkResiduals_.ResizeTo(measurementDim);
0048     kinkResiduals_.Zero();
0049     measResidualErrors_.ResizeTo(measurementDim);
0050     measResidualErrors_.Zero();
0051     kinkResidualErrors_.ResizeTo(measurementDim);
0052     kinkResidualErrors_.Zero();
0053     measDownWeights_.ResizeTo(measurementDim);
0054     measDownWeights_.Zero();
0055     kinkDownWeights_.ResizeTo(measurementDim);
0056     kinkDownWeights_.Zero();
0057
0058     TMatrixDSym zero(repDim);
0059     zero.Zero();
0060
0061     fwdStateCorrection_.ResizeTo(repDim);
0062     fwdCov_.ResizeTo(zero);
0063     bwdStateCorrection_.ResizeTo(repDim);
0064     bwdCov_.ResizeTo(zero);
0065     fwdPrediction_.ResizeTo(repDim);
0066     bwdPrediction_.ResizeTo(repDim);
0067     fwdCov_.SetMatrixArray(zero.GetMatrixArray());
0068     bwdCov_.SetMatrixArray(zero.GetMatrixArray());
0069
0070     measurement_.ResizeTo(2);
0071     measurement_.Zero();
0072     measCov_.ResizeTo(TMatrixDSym(2));
0073     measCov_.Zero();
0074     hMatrix_.ResizeTo(HMatrixUV().getMatrix());
0075     hMatrix_ = HMatrixUV().getMatrix();
0076   }
0077
0078   void GblFitterInfo::setReferenceState(StateOnPlane& referenceState) {
0079     updateMeasurementAndPlane(referenceState);
0080
0081     refPrediction_ = referenceState.getState();
0082     fwdPrediction_ = referenceState.getState();
0083     bwdPrediction_ = referenceState.getState();
0084
0085     //TODO: reset even if already fitted?
0086     fittedStateBwd_.reset(new MeasuredStateOnPlane(referenceState.getState(), trackPoint_->getTrack()->getCovSeed(), sharedPlane_, rep_, referenceState.getAuxInfo()));
0087     fittedStateFwd_.reset(new MeasuredStateOnPlane(referenceState.getState(), trackPoint_->getTrack()->getCovSeed(), sharedPlane_, rep_, referenceState.getAuxInfo()));
0088
0089   }
0090
0091   void GblFitterInfo::setJacobian(TMatrixD jacobian) {
0092     jacobian_.ResizeTo(jacobian);
0093     jacobian_ = jacobian;
0094   }
0095
0096   TMatrixDSym GblFitterInfo::getCovariance(double variance, TVector3 trackDirection, SharedPlanePtr measurementPlane) const {
0097
0098     double c1 = trackDirection.Dot(measurementPlane->getU());
0099     double c2 = trackDirection.Dot(measurementPlane->getV());
0100
0101     TMatrixDSym scatCov(2);
0102     scatCov(0, 0) = 1. - c2 * c2;
0103     scatCov(1, 1) = 1. - c1 * c1;
0104     scatCov(0, 1) = c1 * c2;
0105     scatCov(1, 0) = c1 * c2;
0106     scatCov *= variance * variance / (1. - c1 * c1 - c2 * c2) / (1. - c1 * c1 - c2 * c2) ;
0107
0108     return scatCov;
0109   }
0110
0111   gbl::GblPoint GblFitterInfo::constructGblPoint() {
0112     // All meas/scat info is added from genfit data again (to cope with possible RecoHit update)
0113
0114     // Create GBL point with current jacobian
0115     gbl::GblPoint thePoint(jacobian_);
0116
0117     //NOTE: 3rd update and update anytime GblPoint is requested
0118     // mostly likely will update with reference as on 2nd update
0119     StateOnPlane sop(getFittedState().getState(), sharedPlane_, rep_, getFittedState(true).getAuxInfo());
0120
0121
0122     // Scatterer
0123     if (trackPoint_->hasThinScatterer()) {
0124       if (!hasMeasurements()) {
0125         //double scatVariance = trackPoint_->getMaterialInfo()->getMaterial().getDensity();
0126         //TVectorD kinkPrec(2);
0127         //kinkPrec(0) = 1./scatVariance/scatVariance; kinkPrec(1) = 1./scatVariance/scatVariance;
0128         //thePoint.addScatterer(getKinks(), kinkPrec);
0129         //TODO: if state at scatterer is updated, the direction of track might not be perpendicular anymore
0130         // plane does not change at pure scatterer
0131         TMatrixDSym kinkCov = getCovariance(trackPoint_->getMaterialInfo()->getMaterial().density, sop.getDir(), sop.getPlane());
0132         thePoint.addScatterer(getKinks(), kinkCov.Invert());
0133       } else {
0134         TMatrixDSym kinkCov = getCovariance(trackPoint_->getMaterialInfo()->getMaterial().density, sop.getDir(), sop.getPlane());
0135         thePoint.addScatterer(getKinks(), kinkCov.Invert());
0136       }
0137     }
0138
0139
0140
0141     // Measurement
0142     if (hasMeasurements()) {
0143       MeasuredStateOnPlane measurement = getResidual(0, true, true);
0144       TVectorD aResiduals(measurement.getState());
0145       TMatrixDSym aPrecision(measurement.getCov().Invert());
0146       if (HMatrixU().getMatrix() == hMatrix_) {
0147         double res = aResiduals(0);
0148         double prec = aPrecision(0, 0);
0149         aResiduals.ResizeTo(2);
0150         aPrecision.ResizeTo(TMatrixDSym(2));
0151         aResiduals.Zero();
0152         aResiduals(0) = res;
0153         aPrecision.Zero();
0154         aPrecision(0, 0) = prec;
0155       }
0156       if (HMatrixV().getMatrix() == hMatrix_) {
0157         double res = aResiduals(0);
0158         double prec = aPrecision(0, 0);
0159         aResiduals.ResizeTo(2);
0160         aPrecision.ResizeTo(TMatrixDSym(2));
0161         aResiduals.Zero();
0162         aResiduals(1) = res;
0163         aPrecision.Zero();
0164         aPrecision(1, 1) = prec;
0165       }
0166       // always 2D, U/V set by precision (=0 to disable)
0167       thePoint.addMeasurement(aResiduals, aPrecision);
0168     }
0169
0170     // Derivatives
0171     ICalibrationParametersDerivatives* globals = nullptr;
0172     if (hasMeasurements() && (globals = dynamic_cast<ICalibrationParametersDerivatives*>(trackPoint_->getRawMeasurement(0)) )) {
0173       std::pair<std::vector<int>, TMatrixD> labelsAndMatrix = globals->globalDerivatives(&sop);
0174       std::vector<int> labels = labelsAndMatrix.first;
0175       TMatrixD derivs = labelsAndMatrix.second;
0176
0177       if (derivs.GetNcols() > 0 && !labels.empty() && (unsigned int)derivs.GetNcols() == labels.size()) {
0178         thePoint.addGlobals(labels, derivs);
0179       }
0180       TMatrixD locals = globals->localDerivatives(&sop);
0181       if (locals.GetNcols() > 0) {
0182         thePoint.addLocals(locals);
0183         GblFitStatus* gblfs = dynamic_cast<GblFitStatus*>(trackPoint_->getTrack()->getFitStatus(rep_));
0184         if (gblfs) {
0185           if (gblfs->getMaxLocalFitParams() < locals.GetNcols())
0186             gblfs->setMaxLocalFitParams(locals.GetNcols());
0187         }
0188       }
0189     }
0190
0191     return thePoint;
0192
0193   }
0194
0195   void GblFitterInfo::updateMeasurementAndPlane(const StateOnPlane & sop) {
0196     if (!trackPoint_)
0197       return;
0198     if (!trackPoint_->hasRawMeasurements()) {
0199       //no measurement, plane is supposed to come from state (is perpendicular)
0200       setPlane(sop.getPlane());
0201       return;
0202     }
0203     std::vector<MeasurementOnPlane*> allMeas = trackPoint_->getRawMeasurement(0)->constructMeasurementsOnPlane(sop);
0204
0205     /*
0206     double normMin(9.99E99);
0207     unsigned int imop(0);
0208     const AbsHMatrix* H = allMeas[0]->getHMatrix();
0209     for (unsigned int i=0; i<allMeas.size(); ++i) {
0210       if (*(allMeas[i]->getHMatrix()) != *H){
0211         Exception e("GblFitterInfo::updateMeasurementAndPlane: Cannot compare measurements with different H-Matrices.", __LINE__,__FILE__);
0212         e.setFatal();
0213         throw e;
0214       }
0215
0216       TVectorD res = allMeas[i]->getState() - H->Hv(sop.getState());
0217       double norm = sqrt(res.Norm2Sqr());
0218       if (norm < normMin) {
0219         normMin = norm;
0220         imop = i;
0221       }
0222     }
0223     */
0224     unsigned int imop = 0;
0225     double maxWeight = allMeas.at(0)->getWeight();
0226     for (unsigned int i = 0; i < allMeas.size(); i++)
0227       if (allMeas.at(i)->getWeight() > maxWeight)
0228         imop = i;
0229
0230     measurement_.ResizeTo(allMeas.at(imop)->getState());
0231     measurement_ = allMeas.at(imop)->getState();
0232     measCov_.ResizeTo(allMeas.at(imop)->getCov());
0233     measCov_ = allMeas.at(imop)->getCov();
0234     hMatrix_.ResizeTo(allMeas.at(imop)->getHMatrix()->getMatrix());
0235     hMatrix_ = allMeas.at(imop)->getHMatrix()->getMatrix();
0236
0237     setPlane(allMeas.at(imop)->getPlane());
0238
0239     for (unsigned int imeas = 0; imeas < allMeas.size(); imeas++)
0240       delete allMeas[imeas];
0241     allMeas.clear();
0242
0243
0244   }
0245
0246   void GblFitterInfo::updateFitResults(gbl::GblTrajectory& traj) {
0247     if (!traj.isValid())
0248       return;
0249
0250     // Deduce our position in the trajectory
0251     //----------------------------------------------
0252     unsigned int label = 0;
0253     // Loop over track and find index of this fitter info
0254     genfit::Track* trk = trackPoint_->getTrack();
0255     for (unsigned int ip = 0; ip < trk->getNumPoints(); ip++) {
0256       // Indexing of GblFitter info in track (each gives one point to trajectory at trackpoint position)
0257       if (dynamic_cast<GblFitterInfo*>( trk->getPoint(ip)->getFitterInfo(rep_) )) {
0258         // First fitter info has label 1 (for point in GBL traj)
0259         label++;
0260         // We found itself = we have our label
0261         if (trk->getPoint(ip)->getFitterInfo(rep_) == this)
0262           break;
0263       }
0264     }
0265     if (label == 0)
0266       throw genfit::Exception("GblFitterInfo: fitter info did not found itself in track to update", __LINE__, __FILE__);
0267     if (label > traj.getNumPoints())
0268       throw genfit::Exception("GblFitterInfo: Deduced point label not valid", __LINE__, __FILE__);
0269     //----------------------------------------------
0270
0271     // Residuals (in plane)
0272     unsigned int numMRes = 2;
0273     TVectorD mResiduals(2), mMeasErrors(2), mResErrors(2), mDownWeights(2);
0274     if (0 != traj.getMeasResults(label, numMRes, mResiduals, mMeasErrors, mResErrors, mDownWeights))
0275       throw genfit::Exception(" NO measurement results ", __LINE__, __FILE__);
0276
0277     // Kinks
0278     unsigned int numKRes = 2;
0279     TVectorD kResiduals(2), kMeasErrors(2), kResErrors(2), kDownWeights(2);
0280     if (0 != traj.getScatResults(label, numKRes, kResiduals, kMeasErrors, kResErrors, kDownWeights))
0281       throw genfit::Exception(" NO scattering results ", __LINE__, __FILE__);
0282
0283     // Check for local derivatives
0284     int nLocals = 0;
0285     GblFitStatus* gblfs = dynamic_cast<GblFitStatus*>(trackPoint_->getTrack()->getFitStatus(rep_));
0286     if (gblfs)
0287       nLocals = gblfs->getMaxLocalFitParams();
0288
0289     // Predictions (different at scatterer)
0290     TVectorD bwdUpdate(5 + nLocals), fwdUpdate(5 + nLocals);
0291     TMatrixDSym bwdCov(5 + nLocals), fwdCov(5 + nLocals);
0292
0293     // forward prediction
0294     if (0 != traj.getResults(label, fwdUpdate, fwdCov))
0295       throw genfit::Exception(" NO forward results ", __LINE__, __FILE__);
0296
0297     // backward prediction
0298     if (0 != traj.getResults(-1 * label, bwdUpdate, bwdCov))
0299       throw genfit::Exception(" NO backward results ", __LINE__, __FILE__);
0300
0301     if (nLocals > 0) {
0302       TVectorD _bwdUpdate(5 + nLocals), _fwdUpdate(5 + nLocals);
0303       TMatrixDSym _bwdCov(5 + nLocals), _fwdCov(5 + nLocals);
0304       _bwdUpdate = bwdUpdate;
0305       _fwdUpdate = fwdUpdate;
0306       _bwdCov = bwdCov;
0307       _fwdCov = fwdCov;
0308       bwdUpdate.ResizeTo(5);
0309       fwdUpdate.ResizeTo(5);
0310       bwdCov.ResizeTo(TMatrixDSym(5));
0311       fwdCov.ResizeTo(TMatrixDSym(5));
0312       for (int i = 0; i < 5; i++) {
0313         bwdUpdate(i) = _bwdUpdate(i);
0314         fwdUpdate(i) = _fwdUpdate(i);
0315         for (int j = 0; j < 5; j++) {
0316           bwdCov(i, j) = _bwdCov(i, j);
0317           fwdCov(i, j) = _fwdCov(i, j);
0318         }
0319       }
0320     }
0321     // Now update the the fitter info
0322     //
0323     //-------------------------------------------------
0324     // Backward/forward prediction (residual) (differs at scatterers) AFTER GBL fit
0325     bwdStateCorrection_ = bwdUpdate;
0326     fwdStateCorrection_ = fwdUpdate;
0327     bwdCov_ = bwdCov;
0328     fwdCov_ = fwdCov;
0329     fwdPrediction_ += fwdStateCorrection_; // This is the update!
0330     bwdPrediction_ += bwdStateCorrection_; // This is the update!
0331
0332     fittedStateFwd_.reset( new MeasuredStateOnPlane(fwdPrediction_, fwdCov_, sharedPlane_, rep_, getFittedState(true).getAuxInfo()) );
0333     fittedStateBwd_.reset( new MeasuredStateOnPlane(bwdPrediction_, bwdCov_, sharedPlane_, rep_, getFittedState(true).getAuxInfo()) );
0334
0335     // Set scattering/measurement residual data
0336     kinkResiduals_ = kResiduals;
0337     measResiduals_ = mResiduals;
0338     kinkResidualErrors_ = kResErrors;
0339     measResidualErrors_ = mResErrors;
0340     measDownWeights_ = mDownWeights;
0341     kinkDownWeights_ = kDownWeights;
0342
0343     //-------------------------------------------------
0344   }
0345
0346   void GblFitterInfo::recalculateJacobian(GblFitterInfo* prevFitterInfo)
0347   {
0348     // Invalidates errors and corrections from last iteration
0349     // (will be defined in different plane). But fitted state and residual is ok.
0350
0351     if (!prevFitterInfo) {
0352       jacobian_.UnitMatrix();
0353
0354       //TODO: For 1st plane this has no sense
0355       //return;
0356       prevFitterInfo = this;
0357     }
0358
0359     //TODO
0360     SharedPlanePtr oldPlane(new genfit::DetPlane(*sharedPlane_));
0361
0362     //updateMeasurementAndPlane(StateOnPlane(fwdPrediction_, sharedPlane_, rep_));
0363     //
0364
0365     TMatrixDSym noise;
0366     TVectorD dstate;
0367     // Take forward state from previous fitter info,
0368     // its (maybe updated) plane
0369     // and our rep
0370     StateOnPlane prevState(prevFitterInfo->getFittedState(true).getState(), prevFitterInfo->getPlane(), rep_, getFittedState(true).getAuxInfo());
0371
0372     if (hasMeasurements()) {
0373       SharedPlanePtr newPlane = trackPoint_->getRawMeasurement(0)->constructPlane(prevState);
0374       rep_->extrapolateToPlane(prevState, newPlane, false, true);
0375     } else {
0376       rep_->extrapolateToPlane(prevState, sharedPlane_, false, true);
0377     }
0378
0379     rep_->getForwardJacobianAndNoise(jacobian_, noise, dstate);
0380     // Now update meas data
0381     updateMeasurementAndPlane(prevState);
0382
0383     //
0384     // Extrap predictions to new plane
0385     //
0386     StateOnPlane oldFwdState(fwdPrediction_, oldPlane, rep_, getFittedState(true).getAuxInfo());
0387     StateOnPlane oldBwdState(bwdPrediction_, oldPlane, rep_, getFittedState(true).getAuxInfo());
0388     rep_->extrapolateToPlane(oldFwdState, sharedPlane_);
0389     rep_->extrapolateToPlane(oldBwdState, sharedPlane_);
0390     fwdPrediction_ = oldFwdState.getState();
0391     bwdPrediction_ = oldBwdState.getState();
0392     fittedStateBwd_.reset(new MeasuredStateOnPlane(fwdPrediction_, fwdCov_, sharedPlane_, rep_, getFittedState(true).getAuxInfo()));
0393     fittedStateFwd_.reset(new MeasuredStateOnPlane(bwdPrediction_, bwdCov_, sharedPlane_, rep_, getFittedState(true).getAuxInfo()));
0394     //
0395   }
0396
0397
0398   const MeasuredStateOnPlane& GblFitterInfo::getFittedState(bool afterKink) const {
0399     // ALways biased from GBL (global fit!)
0400
0401     if (!fittedStateFwd_ || !fittedStateBwd_) {
0402       //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
0403       fittedStateFwd_.reset(new MeasuredStateOnPlane(fwdPrediction_, fwdCov_, sharedPlane_, rep_));
0404       fittedStateBwd_.reset(new MeasuredStateOnPlane(bwdPrediction_, bwdCov_, sharedPlane_, rep_));
0405     }
0406
0407     if (afterKink) {
0408       return *fittedStateFwd_;
0409     }
0410     else {
0411       return *fittedStateBwd_;
0412     }
0413
0414   }
0415
0416   MeasurementOnPlane GblFitterInfo::getResidual(unsigned int, bool, bool onlyMeasurementErrors) const {
0417     // Return error of resduals (=0, none, for reference track, no errors known before fit (seed covariance might not be correct)
0418     TMatrixDSym localCovariance(2);
0419     localCovariance.Zero();
0420     // combined meas+fit errors:
0421     //TODO: 1D covariance + more dimensions of residuals
0422
0423     localCovariance(0, 0) = measResidualErrors_(0) * measResidualErrors_(0);
0424     localCovariance(1, 1) = measResidualErrors_(1) * measResidualErrors_(1);
0425
0426     if (HMatrixU().getMatrix() == hMatrix_) {
0427       localCovariance.ResizeTo(TMatrixDSym(1));
0428       localCovariance(0, 0) = measResidualErrors_(0) * measResidualErrors_(0);
0429     }
0430     if (HMatrixV().getMatrix() == hMatrix_) {
0431       localCovariance.ResizeTo(TMatrixDSym(1));
0432       localCovariance(0, 0) = measResidualErrors_(1) * measResidualErrors_(1);
0433     }
0434
0435     if (hasMeasurements()){
0436       // this means only for reference state before gbl fit, this way will be used
0437       TVectorD res( measurement_ - hMatrix_ * getFittedState(true).getState() );
0438
0439       if (onlyMeasurementErrors) {
0440         localCovariance.ResizeTo(measCov_);
0441         localCovariance = measCov_;
0442       }
0443       return MeasurementOnPlane(res, localCovariance, sharedPlane_, rep_, trackPoint_->getRawMeasurement(0)->constructHMatrix(getRep()));
0444     }
0445     TVectorD zeroRes(2);
0446     zeroRes.Zero();
0447     // Else return 0's or whatever
0448     //TODO: or throw?
0449     return MeasurementOnPlane(zeroRes, localCovariance, sharedPlane_, rep_, new HMatrixUV());
0450
0451   } // calculate residual, track and measurement errors are added if onlyMeasurementErrors is false
0452
0453   MeasurementOnPlane GblFitterInfo::getKink() const {
0454     TMatrixDSym localCovariance(2);
0455     localCovariance.Zero();
0456     localCovariance(0, 0) = kinkResidualErrors_(0) * kinkResidualErrors_(0);
0457     localCovariance(1, 1) = kinkResidualErrors_(1) * kinkResidualErrors_(1);
0458     return MeasurementOnPlane(getKinks(), localCovariance, sharedPlane_, rep_, new genfit::HMatrixUV());
0459
0460   }
0461
0462   TVectorD GblFitterInfo::getKinks() const {
0463     TVectorD kinks(2);
0464     kinks.Zero();
0465
0466     TVectorD stateDiff(getFittedState(true).getState() - getFittedState(false).getState());
0467     kinks(0) = -stateDiff(1);
0468     kinks(1) = -stateDiff(2);
0469
0470     return kinks;
0471   }
0472
0473   MeasurementOnPlane GblFitterInfo::getMeasurement() const{
0474     return MeasurementOnPlane(measurement_, measCov_, sharedPlane_, rep_, hasMeasurements() ? trackPoint_->getRawMeasurement(0)->constructHMatrix(rep_) : new HMatrixUV() );
0475   }
0476
0477   bool GblFitterInfo::checkConsistency(const genfit::PruneFlags*) const {
0478     //TODO
0479     return true;
0480   }
0481
0482   GblFitterInfo* GblFitterInfo::clone() const {
0483
0484     GblFitterInfo* retVal = new GblFitterInfo(this->getTrackPoint(), this->getRep());
0485
0486     retVal->setPlane(sharedPlane_);
0487
0488     retVal->jacobian_ = jacobian_;
0489     retVal->measResiduals_ = measResiduals_;
0490     retVal->measResidualErrors_ = measResidualErrors_;
0491     retVal->kinkResiduals_ = kinkResiduals_;
0492     retVal->kinkResidualErrors_ = kinkResidualErrors_;
0493     retVal->measDownWeights_ = measDownWeights_;
0494     retVal->kinkDownWeights_ = kinkDownWeights_;
0495     retVal->bwdStateCorrection_ = bwdStateCorrection_;
0496     retVal->fwdStateCorrection_ = fwdStateCorrection_;
0497     retVal->bwdCov_ = bwdCov_;
0498     retVal->fwdCov_ = fwdCov_;
0499     retVal->fwdPrediction_ = fwdPrediction_;
0500     retVal->bwdPrediction_ = bwdPrediction_;
0501     retVal->refPrediction_ = refPrediction_;
0502     retVal->measurement_.ResizeTo(measurement_);
0503     retVal->measCov_.ResizeTo(measCov_);
0504     retVal->hMatrix_.ResizeTo(hMatrix_);
0505     retVal->measurement_ = measurement_;
0506     retVal->measCov_ = measCov_;
0507     retVal->hMatrix_ = hMatrix_;
0508
0509     return retVal;
0510   }
0511
0512   void GblFitterInfo::Print(const Option_t*) const {
0513     //TODO
0514     std::cout << "=============================================================================================" << std::endl;
0515     std::cout << " >>>  GblFitterInfo " << std::endl;
0516     std::cout << "      ************* " << std::endl;
0517
0518     std::cout << "                      rep: " << rep_ << ", trackpoint: " << trackPoint_ << ", plane: " << sharedPlane_.get() << std::endl;
0519     sharedPlane_->Print();
0520     std::cout << std::endl;
0521
0522     std::cout << "=============================================================================================" << std::endl;
0523     std::cout << "     |          PREDICTIONS            |   REFERENCE    |  Corrections from last iteration  |" << std::endl;
0524     std::cout << "     | (+)prediction  | (-)prediction  |     state      |  (+)correction |  (-) correction  |" << std::endl;
0525     std::cout << "---------------------------------------------------------------------------------------------" << std::endl;
0526
0527     for (int i = 0; i <5; i++) {
0528       std::cout << std::left;
0529       std::cout << " ";
0530       if (i==0)
0531         std::cout << "q/p";
0532       if (i==1)
0533         std::cout << "u' ";
0534       if (i==2)
0535         std::cout << "v' ";
0536       if (i==3)
0537         std::cout << "u  ";
0538       if (i==4)
0539         std::cout << "v  ";
0540       std::cout << " | "
0541       << std::setw(12) << fwdPrediction_(i)              << "   | "
0542       << std::setw(12) << bwdPrediction_(i)              << "   | "
0543       << std::setw(12) << refPrediction_(i)              << "   | "
0544       << std::setw(12) << fwdStateCorrection_(i)           << "   | "
0545       << std::setw(12) << bwdStateCorrection_(i)           << std::endl;
0546     }
0547     std::cout << "=============================================================================================" << std::endl;
0548
0549     if (hasMeasurements()) {
0550       std::cout << "     | Meas. residual |     measurement - prediction    |   Down-weight  |   Fit+meas Err.  |" << std::endl;
0551       std::cout << "     |                |                                 |                |   -diagonaliz.   |" << std::endl;
0552       std::cout << "---------------------------------------------------------------------------------------------" << std::endl;
0553
0554       TVectorD residual = getResidual().getState();
0555       if (residual.GetNoElements()<2) {
0556         residual.ResizeTo(2);
0557         residual.Zero();
0558         if (trackPoint_->getRawMeasurement(0)->constructHMatrix(getRep())->isEqual(genfit::HMatrixU()))
0559           residual(0) = getResidual().getState()(0);
0560         else
0561           residual(1) = getResidual().getState()(0);
0562
0563       }
0564       std::cout << " u   | "
0565       << std::setw(12) << measResiduals_(0)      << "   |  "
0566       << std::setw(12) << residual(0)            << "                   |  "
0567       << std::setw(12) << measDownWeights_(0)    << "  |  "
0568       << std::setw(12) << measResidualErrors_(0)
0569       << std::endl;
0570
0571       std::cout << " v   | "
0572       << std::setw(12) << measResiduals_(1)      << "   |  "
0573       << std::setw(12) << residual(1)            << "                   |  "
0574       << std::setw(12) << measDownWeights_(1)    << "  |  "
0575       << std::setw(12) << measResidualErrors_(1)
0576       << std::endl;
0577
0578       std::cout << "---------------------------------------------------------------------------------------------" << std::endl;
0579     }
0580
0581     std::cout << "     | Kink residual  |  Residual of slope difference   |   Down-weight  |   Fit Kink Err.  |" << std::endl;
0582     std::cout << "     | -diagonalized  |   - (   (+)pred - (-)pred   )   |                |   -diagonaliz.   |" << std::endl;
0583     std::cout << "---------------------------------------------------------------------------------------------" << std::endl;
0584
0585     std::cout << " u'  | "
0586     << std::setw(12) << kinkResiduals_(0)        << "   |  "
0587     << std::setw(12) << getKinks()(0)            << "                   |  "
0588     << std::setw(12) << kinkDownWeights_(0)      << "  |  "
0589     << std::setw(12) << kinkResidualErrors_(0)
0590     << std::endl;
0591
0592     std::cout << " v'  | "
0593     << std::setw(12) << kinkResiduals_(1)        << "   |  "
0594     << std::setw(12) << getKinks()(1)            << "                   |  "
0595     << std::setw(12) << kinkDownWeights_(1)      << "  |  "
0596     << std::setw(12) << kinkResidualErrors_(1)
0597     << std::endl;
0598     std::cout << "=============================================================================================" << std::endl;
0599     std::cout << "Measurement: "; measurement_.Print();
0600
0601     std::cout << "H Matrix: "; hMatrix_.Print();
0602
0603     std::cout << "Measurement covariance: "; measCov_.Print();
0604
0605     std::cout << "Jacobian: "; jacobian_.Print();
0606     std::cout << "Backward covariance: "; bwdCov_.Print();
0607     std::cout << "Forward covariance : "; fwdCov_.Print();
0608
0609     std::cout << "=============================================================================================" << std::endl;
0610
0611   }
0612
0613
0614 } // end of namespace genfit