1 #include "GblFitterInfo.h"
2 #include "MeasurementOnPlane.h"
18 unsigned int dim = rep->
getDim();
20 throw new genfit::Exception(
"GblFitterInfo: Representation state is not 5D", __LINE__, __FILE__);
21 TMatrixDSym noise(dim, dim);
31 unsigned int dim = rep->
getDim();
33 throw new genfit::Exception(
"GblFitterInfo: Representation state is not 5D", __LINE__, __FILE__);
34 TMatrixDSym noise(dim, dim);
42 refPrediction_.ResizeTo(repDim);
43 refPrediction_.Zero();
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();
58 TMatrixDSym zero(repDim);
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());
70 measurement_.ResizeTo(2);
72 measCov_.ResizeTo(TMatrixDSym(2));
74 hMatrix_.ResizeTo(
HMatrixUV().getMatrix());
81 refPrediction_ = referenceState.getState();
82 fwdPrediction_ = referenceState.getState();
83 bwdPrediction_ = referenceState.getState();
92 jacobian_.ResizeTo(jacobian);
98 double c1 = trackDirection.Dot(measurementPlane->getU());
99 double c2 = trackDirection.Dot(measurementPlane->getV());
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) ;
124 if (!hasMeasurements()) {
131 TMatrixDSym kinkCov =
getCovariance(
trackPoint_->getMaterialInfo()->getMaterial().density, sop.getDir(), sop.getPlane());
134 TMatrixDSym kinkCov =
getCovariance(
trackPoint_->getMaterialInfo()->getMaterial().density, sop.getDir(), sop.getPlane());
142 if (hasMeasurements()) {
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));
154 aPrecision(0, 0) = prec;
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));
164 aPrecision(1, 1) = prec;
173 std::pair<std::vector<int>, TMatrixD> labelsAndMatrix = globals->
globalDerivatives(&sop);
174 std::vector<int> labels = labelsAndMatrix.first;
175 TMatrixD derivs = labelsAndMatrix.second;
177 if (derivs.GetNcols() > 0 && !labels.empty() && (
unsigned int)derivs.GetNcols() == labels.size()) {
181 if (locals.GetNcols() > 0) {
185 if (gblfs->getMaxLocalFitParams() < locals.GetNcols())
186 gblfs->setMaxLocalFitParams(locals.GetNcols());
200 setPlane(sop.getPlane());
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)
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();
237 setPlane(allMeas.at(imop)->getPlane());
239 for (
unsigned int imeas = 0; imeas < allMeas.size(); imeas++)
240 delete allMeas[imeas];
252 unsigned int label = 0;
255 for (
unsigned int ip = 0; ip < trk->getNumPoints(); ip++) {
266 throw genfit::Exception(
"GblFitterInfo: fitter info did not found itself in track to update", __LINE__, __FILE__);
268 throw genfit::Exception(
"GblFitterInfo: Deduced point label not valid", __LINE__, __FILE__);
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))
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))
287 nLocals = gblfs->getMaxLocalFitParams();
290 TVectorD bwdUpdate(5 + nLocals), fwdUpdate(5 + nLocals);
291 TMatrixDSym bwdCov(5 + nLocals), fwdCov(5 + nLocals);
294 if (0 != traj.
getResults(label, fwdUpdate, fwdCov))
298 if (0 != traj.
getResults(-1 * label, bwdUpdate, bwdCov))
302 TVectorD _bwdUpdate(5 + nLocals), _fwdUpdate(5 + nLocals);
303 TMatrixDSym _bwdCov(5 + nLocals), _fwdCov(5 + nLocals);
304 _bwdUpdate = bwdUpdate;
305 _fwdUpdate = fwdUpdate;
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);
325 bwdStateCorrection_ = bwdUpdate;
326 fwdStateCorrection_ = fwdUpdate;
329 fwdPrediction_ += fwdStateCorrection_;
330 bwdPrediction_ += bwdStateCorrection_;
336 kinkResiduals_ = kResiduals;
337 measResiduals_ = mResiduals;
338 kinkResidualErrors_ = kResErrors;
339 measResidualErrors_ = mResErrors;
340 measDownWeights_ = mDownWeights;
341 kinkDownWeights_ = kDownWeights;
351 if (!prevFitterInfo) {
352 jacobian_.UnitMatrix();
356 prevFitterInfo =
this;
372 if (hasMeasurements()) {
390 fwdPrediction_ = oldFwdState.getState();
391 bwdPrediction_ = oldBwdState.getState();
411 return *fittedStateBwd_;
418 TMatrixDSym localCovariance(2);
419 localCovariance.Zero();
423 localCovariance(0, 0) = measResidualErrors_(0) * measResidualErrors_(0);
424 localCovariance(1, 1) = measResidualErrors_(1) * measResidualErrors_(1);
426 if (
HMatrixU().getMatrix() == hMatrix_) {
427 localCovariance.ResizeTo(TMatrixDSym(1));
428 localCovariance(0, 0) = measResidualErrors_(0) * measResidualErrors_(0);
430 if (
HMatrixV().getMatrix() == hMatrix_) {
431 localCovariance.ResizeTo(TMatrixDSym(1));
432 localCovariance(0, 0) = measResidualErrors_(1) * measResidualErrors_(1);
435 if (hasMeasurements()){
437 TVectorD res( measurement_ - hMatrix_ *
getFittedState(
true).getState() );
439 if (onlyMeasurementErrors) {
440 localCovariance.ResizeTo(measCov_);
441 localCovariance = measCov_;
454 TMatrixDSym localCovariance(2);
455 localCovariance.Zero();
456 localCovariance(0, 0) = kinkResidualErrors_(0) * kinkResidualErrors_(0);
457 localCovariance(1, 1) = kinkResidualErrors_(1) * kinkResidualErrors_(1);
467 kinks(0) = -stateDiff(1);
468 kinks(1) = -stateDiff(2);
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_;
512 void GblFitterInfo::Print(
const Option_t*)
const {
514 std::cout <<
"=============================================================================================" << std::endl;
515 std::cout <<
" >>> GblFitterInfo " << std::endl;
516 std::cout <<
" ************* " << std::endl;
520 std::cout << std::endl;
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;
527 for (
int i = 0; i <5; i++) {
528 std::cout << std::left;
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;
547 std::cout <<
"=============================================================================================" << std::endl;
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;
555 if (residual.GetNoElements()<2) {
556 residual.ResizeTo(2);
565 << std::setw(12) << measResiduals_(0) <<
" | "
566 << std::setw(12) << residual(0) <<
" | "
567 << std::setw(12) << measDownWeights_(0) <<
" | "
568 << std::setw(12) << measResidualErrors_(0)
572 << std::setw(12) << measResiduals_(1) <<
" | "
573 << std::setw(12) << residual(1) <<
" | "
574 << std::setw(12) << measDownWeights_(1) <<
" | "
575 << std::setw(12) << measResidualErrors_(1)
578 std::cout <<
"---------------------------------------------------------------------------------------------" << std::endl;
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;
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)
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)
598 std::cout <<
"=============================================================================================" << std::endl;
599 std::cout <<
"Measurement: "; measurement_.Print();
601 std::cout <<
"H Matrix: "; hMatrix_.Print();
603 std::cout <<
"Measurement covariance: "; measCov_.Print();
605 std::cout <<
"Jacobian: "; jacobian_.Print();
606 std::cout <<
"Backward covariance: "; bwdCov_.Print();
607 std::cout <<
"Forward covariance : "; fwdCov_.Print();
609 std::cout <<
"=============================================================================================" << std::endl;