| File: | genfit2/code2/GBL/src/GblFitterInfo2.cc |
| Warning: | line 343, column 7 Value stored to 'prevFitterInfo' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
| 1 | #include "GblFitterInfo2.h" |
| 2 | #include "MeasurementOnPlane.h" |
| 3 | #include "HMatrixUV.h" |
| 4 | #include "HMatrixV.h" |
| 5 | #include "HMatrixU.h" |
| 6 | |
| 7 | //! Namespace for the general broken lines package |
| 8 | namespace genfit { |
| 9 | |
| 10 | GblFitterInfo2::GblFitterInfo2() : AbsFitterInfo(), jacobian_(TMatrixD(5, 5)), noise_(TMatrixDSym(5)) { |
| 11 | reset(); |
| 12 | } |
| 13 | |
| 14 | GblFitterInfo2::GblFitterInfo2(const TrackPoint* trackPoint, const AbsTrackRep* rep) : AbsFitterInfo(trackPoint, rep), jacobian_(TMatrixD(5, 5)), noise_(TMatrixDSym(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("GblFitterInfo2: Representation state is not 5D", __LINE__20, __FILE__"genfit2/code2/GBL/src/GblFitterInfo2.cc"); |
| 21 | rep->getForwardJacobianAndNoise(jacobian_, noise_); |
| 22 | |
| 23 | } |
| 24 | |
| 25 | GblFitterInfo2::GblFitterInfo2(const TrackPoint* trackPoint, const AbsTrackRep* rep, StateOnPlane& referenceState) : AbsFitterInfo(trackPoint, rep), jacobian_(TMatrixD(5, 5)) { |
| 26 | reset(); |
| 27 | |
| 28 | // Retrieve jacobian from rep by default (may not be used ... at 1st point) |
| 29 | unsigned int dim = rep->getDim(); |
| 30 | if (dim != 5) |
| 31 | throw new genfit::Exception("GblFitterInfo2: Representation state is not 5D", __LINE__31, __FILE__"genfit2/code2/GBL/src/GblFitterInfo2.cc"); |
| 32 | rep->getForwardJacobianAndNoise(jacobian_, noise_); |
| 33 | |
| 34 | setReferenceState(referenceState); |
| 35 | } |
| 36 | |
| 37 | void GblFitterInfo2::reset(unsigned int repDim) { |
| 38 | label_ = 0; |
| 39 | measType_ = EMeasurementType::noMeas; |
| 40 | refPrediction_.ResizeTo(repDim); |
| 41 | refPrediction_.Zero(); |
| 42 | |
| 43 | TMatrixDSym zero(repDim); |
| 44 | zero.Zero(); |
| 45 | |
| 46 | fwdStateCorrection_.ResizeTo(repDim); |
| 47 | fwdCov_.ResizeTo(zero); |
| 48 | bwdStateCorrection_.ResizeTo(repDim); |
| 49 | bwdCov_.ResizeTo(zero); |
| 50 | diffFwdBwdPrediction_.ResizeTo(repDim); |
| 51 | diffFwdBwdPrediction_.Zero(); |
| 52 | fwdPrediction_.ResizeTo(repDim); |
| 53 | bwdPrediction_.ResizeTo(repDim); |
| 54 | fwdCov_.SetMatrixArray(zero.GetMatrixArray()); |
| 55 | bwdCov_.SetMatrixArray(zero.GetMatrixArray()); |
| 56 | |
| 57 | measUsed_.clear(); |
| 58 | measDownWeight_.clear(); |
| 59 | measurement_.clear(); |
| 60 | measCov_.clear(); |
| 61 | hMatrix_.ResizeTo(HMatrixUV().getMatrix()); |
| 62 | hMatrix_ = HMatrixUV().getMatrix(); |
| 63 | } |
| 64 | |
| 65 | void GblFitterInfo2::setReferenceState(StateOnPlane& referenceState) { |
| 66 | updateMeasurementAndPlane(referenceState); |
| 67 | |
| 68 | refPrediction_ = referenceState.getState(); |
| 69 | fwdPrediction_ = referenceState.getState(); |
| 70 | bwdPrediction_ = referenceState.getState(); |
| 71 | |
| 72 | //TODO: reset even if already fitted? |
| 73 | fittedStateBwd_.reset(new MeasuredStateOnPlane(referenceState.getState(), trackPoint_->getTrack()->getCovSeed(), sharedPlane_, rep_, referenceState.getAuxInfo())); |
| 74 | fittedStateFwd_.reset(new MeasuredStateOnPlane(referenceState.getState(), trackPoint_->getTrack()->getCovSeed(), sharedPlane_, rep_, referenceState.getAuxInfo())); |
| 75 | |
| 76 | } |
| 77 | |
| 78 | |
| 79 | gbl::GblPoint GblFitterInfo2::constructGblPoint(unsigned int index, bool allowAmbiguities) { |
| 80 | // All meas/scat info is added from genfit data again (to cope with possible RecoHit update) |
| 81 | |
| 82 | // Create GBL point with current jacobian |
| 83 | gbl::GblPoint thePoint(jacobian_); |
| 84 | |
| 85 | // Scatterer (all except first point) |
| 86 | if (index > 0) { |
| 87 | TMatrixDSym scatCov(4); |
| 88 | // assuming no correlation between energy and scattering variance |
| 89 | noise_.GetSub(1,4,scatCov, "S"); |
| 90 | thePoint.addThickScatterer(getKinksAndSteps(), scatCov.Invert()); |
| 91 | } |
| 92 | |
| 93 | // (implemented type of) Measurement |
| 94 | if (measType_ != EMeasurementType::noMeas) { |
| 95 | unsigned int numMeasAllowed = allowAmbiguities ? measUsed_.size() : 1; |
| 96 | for (unsigned int iMeas = 0; iMeas < numMeasAllowed; iMeas++) { |
| 97 | TVectorD aRes(measurement_.at(iMeas) - hMatrix_ * fwdPrediction_); |
| 98 | TMatrixDSym aCov(measCov_.at(iMeas)); |
| 99 | if (measType_ < EMeasurementType::uvMeas) { |
| 100 | // (1D) U or V measurement |
| 101 | TVectorD aResidual(2); |
| 102 | aResidual(static_cast<int>(measType_)) = aRes(0); |
| 103 | aResidual(1-static_cast<int>(measType_)) = 0.; |
| 104 | TVectorD aPrecision(2); |
| 105 | aPrecision(static_cast<int>(measType_)) = 1. / aCov(0,0); |
| 106 | aPrecision(1-static_cast<int>(measType_)) = 0.; // disable other measurement |
| 107 | thePoint.addMeasurement(aResidual, aPrecision); |
| 108 | } else { |
| 109 | // (at least 2D) U,V measurement, arbitrary precision matrix |
| 110 | thePoint.addMeasurement(aRes, aCov.Invert()); |
| 111 | } |
| 112 | // init downweight |
| 113 | measDownWeight_.at(measUsed_.at(iMeas)) = 1.; |
| 114 | |
| 115 | // Derivatives |
| 116 | ICalibrationParametersDerivatives* globals = nullptr; |
| 117 | if ((globals = dynamic_cast<ICalibrationParametersDerivatives*>(trackPoint_->getRawMeasurement(0)) )) { |
| 118 | StateOnPlane sop(getFittedState().getState(), sharedPlane_, rep_, getFittedState(true).getAuxInfo()); |
| 119 | std::pair<std::vector<int>, TMatrixD> labelsAndMatrix = globals->globalDerivatives(&sop); |
| 120 | std::vector<int> labels = labelsAndMatrix.first; |
| 121 | TMatrixD derivs = labelsAndMatrix.second; |
| 122 | |
| 123 | if (derivs.GetNcols() > 0 && !labels.empty() && (unsigned int)derivs.GetNcols() == labels.size()) { |
| 124 | thePoint.addGlobals(labels, derivs); |
| 125 | } |
| 126 | TMatrixD locals = globals->localDerivatives(&sop); |
| 127 | if (locals.GetNcols() > 0) { |
| 128 | thePoint.addLocals(locals); |
| 129 | GblFitStatus* gblfs = dynamic_cast<GblFitStatus*>(trackPoint_->getTrack()->getFitStatus(rep_)); |
| 130 | if (gblfs) { |
| 131 | if (gblfs->getMaxLocalFitParams() < locals.GetNcols()) |
| 132 | gblfs->setMaxLocalFitParams(locals.GetNcols()); |
| 133 | } |
| 134 | } |
| 135 | } |
| 136 | } |
| 137 | } |
| 138 | |
| 139 | return thePoint; |
| 140 | |
| 141 | } |
| 142 | |
| 143 | void GblFitterInfo2::updateMeasurementAndPlane(const StateOnPlane & sop) { |
| 144 | if (!trackPoint_) |
| 145 | return; |
| 146 | if (!trackPoint_->hasRawMeasurements()) { |
| 147 | //no measurement, plane is supposed to come from state (is perpendicular) |
| 148 | setPlane(sop.getPlane()); |
| 149 | return; |
| 150 | } |
| 151 | std::vector<MeasurementOnPlane*> allMeas = trackPoint_->getRawMeasurement(0)->constructMeasurementsOnPlane(sop); |
| 152 | |
| 153 | // initial call ? |
| 154 | bool initial = measDownWeight_.empty(); |
| 155 | if (initial) { |
| 156 | // reset GBL downweights |
| 157 | for (unsigned int i = 0; i < allMeas.size(); i++) {measDownWeight_.push_back(0.);} |
| 158 | // get most probable measurement |
| 159 | unsigned int imop = 0; |
| 160 | double maxWeight = allMeas.at(0)->getWeight(); |
| 161 | for (unsigned int i = 1; i < allMeas.size(); i++) |
| 162 | if (allMeas.at(i)->getWeight() > maxWeight) { |
| 163 | imop = i; |
| 164 | maxWeight = allMeas.at(i)->getWeight(); |
| 165 | } |
| 166 | // use that |
| 167 | measUsed_.push_back(imop); |
| 168 | } |
| 169 | |
| 170 | unsigned int imop = measUsed_.at(0); |
| 171 | hMatrix_.ResizeTo(allMeas.at(imop)->getHMatrix()->getMatrix()); |
| 172 | hMatrix_ = allMeas.at(imop)->getHMatrix()->getMatrix(); |
| 173 | // (implemented) measurement type(s) |
| 174 | if (HMatrixU().getMatrix() == hMatrix_) { measType_ = EMeasurementType::uMeas; } |
| 175 | else if (HMatrixV().getMatrix() == hMatrix_) { measType_ = EMeasurementType::vMeas; } |
| 176 | else if (HMatrixUV().getMatrix() == hMatrix_) { measType_ = EMeasurementType::uvMeas; } |
| 177 | else { |
| 178 | //no implemented measurement type, plane is supposed to come from state (is perpendicular) |
| 179 | #ifdef DEBUG |
| 180 | std::cout << " no implemented meas type " << std::endl; |
| 181 | #endif |
| 182 | setPlane(sop.getPlane()); |
| 183 | return; |
| 184 | } |
| 185 | |
| 186 | // check for ambiguities |
| 187 | if (initial) { |
| 188 | // ambiguity (CDC) ? |
| 189 | if ((allMeas.size() == 2) && (measType_ == EMeasurementType::uMeas)) { |
| 190 | // use that too |
| 191 | measUsed_.push_back(1-imop); |
| 192 | } |
| 193 | } |
| 194 | |
| 195 | // measurements to be used |
| 196 | for (unsigned int i = 0; i < measUsed_.size(); i++) { |
| 197 | unsigned int index = measUsed_.at(i); |
| 198 | if (initial) { |
| 199 | measurement_.push_back(allMeas.at(index)->getState()); |
| 200 | measCov_.push_back(allMeas.at(index)->getCov()); |
| 201 | } else { |
| 202 | measurement_.at(i) = allMeas.at(index)->getState(); |
| 203 | measCov_.at(i) = allMeas.at(index)->getCov(); |
| 204 | } |
| 205 | } |
| 206 | |
| 207 | setPlane(allMeas.at(imop)->getPlane()); |
| 208 | |
| 209 | for (unsigned int imeas = 0; imeas < allMeas.size(); imeas++) |
| 210 | delete allMeas[imeas]; |
| 211 | allMeas.clear(); |
| 212 | |
| 213 | } |
| 214 | |
| 215 | void GblFitterInfo2::resolveAmbiguities(gbl::GblTrajectory& traj, gbl::GblPoint& point, unsigned int &nResolved, unsigned int &nSwapped) { |
| 216 | if (!traj.isValid()) |
| 217 | return; |
| 218 | |
| 219 | // Implemented: one or two u or v measurement (-> max 2 meas. -> size of mResiduals etc) |
| 220 | |
| 221 | // (1 or 2) 1D u or v measurement |
| 222 | if (measType_ < EMeasurementType::uvMeas) { |
| 223 | unsigned int numMRes = 0; |
| 224 | TVectorD mResiduals(2), mMeasErrors(2), mResErrors(2), mDownWeights(2); |
| 225 | if (0 != traj.getMeasResults(label_, numMRes, mResiduals, mMeasErrors, mResErrors, mDownWeights)) |
| 226 | throw genfit::Exception(" NO measurement results ", __LINE__226, __FILE__"genfit2/code2/GBL/src/GblFitterInfo2.cc"); |
| 227 | |
| 228 | // resolve side |
| 229 | if (numMRes == 2) { |
| 230 | if (fabs(mResiduals[0]) <= fabs(mResiduals[1])) |
| 231 | (point.getMeasBegin() + 1)->setEnabled(false); // disable measurement 1 |
| 232 | else { |
| 233 | point.getMeasBegin()->setEnabled(false); // disable measurement 0 |
| 234 | // swap measurements (to restore order) |
| 235 | unsigned int imop = measUsed_.at(0); |
| 236 | measUsed_.at(0) = 1-imop; |
| 237 | measUsed_.at(1) = imop; |
| 238 | nSwapped++; |
| 239 | #ifdef DEBUG |
| 240 | std::cout << " swap " << label_ << " " << mResiduals[0] << " " << mResiduals[1] << std::endl; |
| 241 | #endif |
| 242 | } |
| 243 | nResolved++; |
| 244 | } |
| 245 | } |
| 246 | return; |
| 247 | } |
| 248 | |
| 249 | void GblFitterInfo2::updateDownweights(gbl::GblTrajectory& traj) { |
| 250 | if (!traj.isValid()) |
| 251 | return; |
| 252 | |
| 253 | // Implemented: single uv or one or two u or v measurement (-> max 2 meas. -> size of mResiduals etc) |
| 254 | unsigned int numMRes = 0; |
| 255 | TVectorD mResiduals(2), mMeasErrors(2), mResErrors(2), mDownWeights(2); |
| 256 | if (0 != traj.getMeasResults(label_, numMRes, mResiduals, mMeasErrors, mResErrors, mDownWeights)) |
| 257 | throw genfit::Exception(" NO measurement results ", __LINE__257, __FILE__"genfit2/code2/GBL/src/GblFitterInfo2.cc"); |
| 258 | |
| 259 | // (1 or 2) 1D u or v measurement |
| 260 | if (measType_ < EMeasurementType::uvMeas) { |
| 261 | for (unsigned int i = 0; i < numMRes; ++i) { |
| 262 | measDownWeight_.at(measUsed_.at(i)) = mDownWeights(i); |
| 263 | } |
| 264 | // reset weights for discarded measurements |
| 265 | for (unsigned int i = numMRes; i < measUsed_.size(); ++i) { |
| 266 | measDownWeight_.at(measUsed_.at(i)) = 0.; |
| 267 | } |
| 268 | } |
| 269 | // (single) 2D uv measurement (2 downweights, use product) |
| 270 | if (measType_ == EMeasurementType::uvMeas && numMRes == 2) { |
| 271 | measDownWeight_.at(measUsed_.at(0)) = mDownWeights(0)*mDownWeights(1); |
| 272 | } |
| 273 | } |
| 274 | |
| 275 | void GblFitterInfo2::updateFitResults(gbl::GblTrajectory& traj) { |
| 276 | if (!traj.isValid()) |
| 277 | return; |
| 278 | |
| 279 | unsigned int label = label_; |
| 280 | |
| 281 | // Now update the the fitter info |
| 282 | // |
| 283 | //------------------------------------------------- |
| 284 | // Backward/forward prediction (residual) (differs at scatterers) AFTER GBL fit |
| 285 | |
| 286 | // getResults respects size of ROOT matrices while converting from (internal) EIGEN3 matrices |
| 287 | // forward prediction |
| 288 | if (0 != traj.getResults(label, fwdStateCorrection_, fwdCov_)) |
| 289 | throw genfit::Exception(" NO forward results ", __LINE__289, __FILE__"genfit2/code2/GBL/src/GblFitterInfo2.cc"); |
| 290 | |
| 291 | // backward prediction |
| 292 | if (0 != traj.getResults(-1 * label, bwdStateCorrection_, bwdCov_)) |
| 293 | throw genfit::Exception(" NO backward results ", __LINE__293, __FILE__"genfit2/code2/GBL/src/GblFitterInfo2.cc"); |
| 294 | |
| 295 | diffFwdBwdPrediction_ += fwdStateCorrection_ - bwdStateCorrection_; |
| 296 | fwdPrediction_ += fwdStateCorrection_; // This is the update! |
| 297 | bwdPrediction_ += bwdStateCorrection_; // This is the update! |
| 298 | |
| 299 | fittedStateFwd_.reset( new MeasuredStateOnPlane(fwdPrediction_, fwdCov_, sharedPlane_, rep_, getFittedState(true).getAuxInfo()) ); |
| 300 | fittedStateBwd_.reset( new MeasuredStateOnPlane(bwdPrediction_, bwdCov_, sharedPlane_, rep_, getFittedState(true).getAuxInfo()) ); |
| 301 | |
| 302 | //---------------------------------------------- |
| 303 | |
| 304 | #ifdef DEBUG |
| 305 | // Residuals (in plane) |
| 306 | // Implemented: single uv or one or two u or v measurement (-> max 2 meas. -> size of mResiduals etc) |
| 307 | unsigned int numMRes = 0; |
| 308 | TVectorD mResiduals(2), mMeasErrors(2), mResErrors(2), mDownWeights(2); |
| 309 | if (0 != traj.getMeasResults(label, numMRes, mResiduals, mMeasErrors, mResErrors, mDownWeights)) |
| 310 | throw genfit::Exception(" NO measurement results ", __LINE__310, __FILE__"genfit2/code2/GBL/src/GblFitterInfo2.cc"); |
| 311 | |
| 312 | // Kinks and steps |
| 313 | unsigned int numKRes = 0; |
| 314 | TVectorD kResiduals(4), kMeasErrors(4), kResErrors(4), kDownWeights(4); |
| 315 | if (0 != traj.getScatResults(label, numKRes, kResiduals, kMeasErrors, kResErrors, kDownWeights)) |
| 316 | throw genfit::Exception(" NO scattering results ", __LINE__316, __FILE__"genfit2/code2/GBL/src/GblFitterInfo2.cc"); |
| 317 | |
| 318 | std::cout << std::endl; |
| 319 | std::cout << " GBL residuals at label " << label << std::endl; |
| 320 | for (unsigned int i = 0; i < numMRes; ++i) { |
| 321 | std::cout << " meas " << std::setw(15) << mResiduals(i) << " " << std::setw(15) << mMeasErrors(i) << " " |
| 322 | << std::setw(15) << mResErrors(i) << " " << std::setw(15) << mDownWeights(i) << std::endl; |
| 323 | } |
| 324 | for (unsigned int i = 0; i < numKRes; ++i) { |
| 325 | std::cout << " scat " << std::setw(15) << kResiduals(i) << " " << std::setw(15) << kMeasErrors(i) << " " |
| 326 | << std::setw(15) << kResErrors(i) << " " << std::setw(15) << kDownWeights(i) << std::endl; |
| 327 | } |
| 328 | std::cout << std::endl; |
| 329 | #endif |
| 330 | |
| 331 | //------------------------------------------------- |
| 332 | } |
| 333 | |
| 334 | void GblFitterInfo2::recalculateJacobian(GblFitterInfo2* prevFitterInfo) |
| 335 | { |
| 336 | // Invalidates errors and corrections from last iteration |
| 337 | // (will be defined in different plane). But fitted state and residual is ok. |
| 338 | |
| 339 | if (!prevFitterInfo) { |
| 340 | jacobian_.UnitMatrix(); |
| 341 | noise_.Zero(); |
| 342 | |
| 343 | prevFitterInfo = this; |
Value stored to 'prevFitterInfo' is never read | |
| 344 | return; |
| 345 | } |
| 346 | |
| 347 | //TODO |
| 348 | //updateMeasurementAndPlane(StateOnPlane(fwdPrediction_, sharedPlane_, rep_)); |
| 349 | // |
| 350 | |
| 351 | // Take forward state from previous fitter info, |
| 352 | // its (maybe updated) plane |
| 353 | // and our rep |
| 354 | StateOnPlane prevState(prevFitterInfo->getFittedState(true).getState(), prevFitterInfo->getPlane(), rep_, getFittedState(true).getAuxInfo()); |
| 355 | |
| 356 | if (hasMeasurements()) { |
| 357 | SharedPlanePtr newPlane = trackPoint_->getRawMeasurement(0)->constructPlane(prevState); |
| 358 | rep_->extrapolateToPlane(prevState, newPlane, false, true); |
| 359 | } else { |
| 360 | rep_->extrapolateToPlane(prevState, sharedPlane_, false, true); |
| 361 | } |
| 362 | |
| 363 | rep_->getForwardJacobianAndNoise(jacobian_, noise_); |
| 364 | // Now update meas data |
| 365 | updateMeasurementAndPlane(prevState); |
| 366 | |
| 367 | // |
| 368 | // Extrapolate predictions to new plane |
| 369 | // |
| 370 | bwdPrediction_ = prevState.getState(); |
| 371 | fwdPrediction_ = bwdPrediction_ + diffFwdBwdPrediction_; |
| 372 | fittedStateFwd_.reset(new MeasuredStateOnPlane(fwdPrediction_, fwdCov_, sharedPlane_, rep_, getFittedState(true).getAuxInfo())); |
| 373 | fittedStateBwd_.reset(new MeasuredStateOnPlane(bwdPrediction_, bwdCov_, sharedPlane_, rep_, getFittedState(true).getAuxInfo())); |
| 374 | // |
| 375 | } |
| 376 | |
| 377 | |
| 378 | const MeasuredStateOnPlane& GblFitterInfo2::getFittedState(bool afterKink) const { |
| 379 | // ALways biased from GBL (global fit!) |
| 380 | |
| 381 | if (!fittedStateFwd_ || !fittedStateBwd_) { |
| 382 | //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 |
| 383 | fittedStateFwd_.reset(new MeasuredStateOnPlane(fwdPrediction_, fwdCov_, sharedPlane_, rep_)); |
| 384 | fittedStateBwd_.reset(new MeasuredStateOnPlane(bwdPrediction_, bwdCov_, sharedPlane_, rep_)); |
| 385 | } |
| 386 | |
| 387 | if (afterKink) { |
| 388 | return *fittedStateFwd_; |
| 389 | } |
| 390 | else { |
| 391 | return *fittedStateBwd_; |
| 392 | } |
| 393 | |
| 394 | } |
| 395 | |
| 396 | MeasurementOnPlane GblFitterInfo2::getResidual(unsigned int iMeas, bool, bool) const { |
| 397 | // only measurement errors implemented: |
| 398 | //TODO: 1D covariance + more dimensions of residuals |
| 399 | |
| 400 | if (hasMeasurements()){ |
| 401 | // this means only for reference state before gbl fit, this way will be used |
| 402 | TVectorD res( measurement_.at(iMeas) - hMatrix_ * fwdPrediction_ ); |
| 403 | return MeasurementOnPlane(res, measCov_.at(iMeas), sharedPlane_, rep_, trackPoint_->getRawMeasurement(0)->constructHMatrix(getRep())); |
| 404 | } |
| 405 | TMatrixDSym zeroCov(2); |
| 406 | zeroCov.Zero(); |
| 407 | TVectorD zeroRes(2); |
| 408 | zeroRes.Zero(); |
| 409 | // Else return 0's or whatever |
| 410 | //TODO: or throw? |
| 411 | return MeasurementOnPlane(zeroRes, zeroCov, sharedPlane_, rep_, new HMatrixUV()); |
| 412 | |
| 413 | } // calculate residual |
| 414 | |
| 415 | TVectorD GblFitterInfo2::getKinksAndSteps() const { |
| 416 | // return kinks and steps residuals = measurement(=0.) - prediction |
| 417 | TVectorD kinksAndSteps(4); |
| 418 | for (unsigned int i = 0; i < 4; ++i) { kinksAndSteps(i)= -diffFwdBwdPrediction_(i+1); } |
| 419 | return kinksAndSteps; |
| 420 | } |
| 421 | |
| 422 | MeasurementOnPlane GblFitterInfo2::getMeasurement(unsigned int iMeas) const{ |
| 423 | return MeasurementOnPlane(measurement_.at(iMeas), measCov_.at(iMeas), sharedPlane_, rep_, hasMeasurements() ? trackPoint_->getRawMeasurement(0)->constructHMatrix(rep_) : new HMatrixUV() ); |
| 424 | } |
| 425 | |
| 426 | bool GblFitterInfo2::checkConsistency(const genfit::PruneFlags*) const { |
| 427 | //TODO |
| 428 | return true; |
| 429 | } |
| 430 | |
| 431 | GblFitterInfo2* GblFitterInfo2::clone() const { |
| 432 | |
| 433 | GblFitterInfo2* retVal = new GblFitterInfo2(this->getTrackPoint(), this->getRep()); |
| 434 | |
| 435 | retVal->setPlane(sharedPlane_); |
| 436 | |
| 437 | retVal->measType_ = measType_; |
| 438 | retVal->label_ = label_; |
| 439 | retVal->jacobian_ = jacobian_; |
| 440 | retVal->noise_ = noise_; |
| 441 | retVal->bwdStateCorrection_ = bwdStateCorrection_; |
| 442 | retVal->fwdStateCorrection_ = fwdStateCorrection_; |
| 443 | retVal->bwdCov_ = bwdCov_; |
| 444 | retVal->fwdCov_ = fwdCov_; |
| 445 | retVal->diffFwdBwdPrediction_ = diffFwdBwdPrediction_; |
| 446 | retVal->fwdPrediction_ = fwdPrediction_; |
| 447 | retVal->bwdPrediction_ = bwdPrediction_; |
| 448 | retVal->refPrediction_ = refPrediction_; |
| 449 | retVal->measUsed_ = measUsed_; |
| 450 | retVal->measDownWeight_ = measDownWeight_; |
| 451 | retVal->measurement_ = measurement_; |
| 452 | retVal->measCov_ = measCov_; |
| 453 | retVal->hMatrix_.ResizeTo(hMatrix_); |
| 454 | retVal->hMatrix_ = hMatrix_; |
| 455 | |
| 456 | return retVal; |
| 457 | } |
| 458 | |
| 459 | void GblFitterInfo2::Print(const Option_t*) const { |
| 460 | //TODO |
| 461 | std::cout << "=============================================================================================" << std::endl; |
| 462 | std::cout << " >>> GblFitterInfo2 " << std::endl; |
| 463 | std::cout << " ************** " << std::endl; |
| 464 | |
| 465 | std::cout << " rep: " << rep_ << ", trackpoint: " << trackPoint_ << ", plane: " << sharedPlane_.get() << std::endl; |
| 466 | sharedPlane_->Print(); |
| 467 | std::cout << std::endl; |
| 468 | |
| 469 | std::cout << "===========================================================================================" << std::endl; |
| 470 | std::cout << " | PREDICTIONS | REFERENCE | Corrections from last iteration |" << std::endl; |
| 471 | std::cout << " | (+)prediction | (-)prediction | state | (+)correction | (-) correction |" << std::endl; |
| 472 | std::cout << "-------------------------------------------------------------------------------------------" << std::endl; |
| 473 | |
| 474 | for (int i = 0; i <5; i++) { |
| 475 | std::cout << std::left; |
| 476 | std::cout << " "; |
| 477 | if (i==0) |
| 478 | std::cout << "q/p"; |
| 479 | if (i==1) |
| 480 | std::cout << "u' "; |
| 481 | if (i==2) |
| 482 | std::cout << "v' "; |
| 483 | if (i==3) |
| 484 | std::cout << "u "; |
| 485 | if (i==4) |
| 486 | std::cout << "v "; |
| 487 | std::cout << " | " |
| 488 | << std::setw(12) << fwdPrediction_(i) << " | " |
| 489 | << std::setw(12) << bwdPrediction_(i) << " | " |
| 490 | << std::setw(12) << refPrediction_(i) << " | " |
| 491 | << std::setw(12) << fwdStateCorrection_(i) << " | " |
| 492 | << std::setw(12) << bwdStateCorrection_(i) << " | " << std::endl; |
| 493 | } |
| 494 | std::cout << "===========================================================================================" << std::endl; |
| 495 | |
| 496 | TVectorD residual(getKinksAndSteps()); |
| 497 | std::cout << " | Scat. res. u' | Scat. res. v' | Scat. res. u | Scat. res. v |" << std::endl; |
| 498 | std::cout << "--------------------------------------------------------------------------" << std::endl; |
| 499 | std::cout << " | " |
| 500 | << std::setw(12) << residual(0) << " | " |
| 501 | << std::setw(12) << residual(1) << " | " |
| 502 | << std::setw(12) << residual(2) << " | " |
| 503 | << std::setw(12) << residual(3) << " | " |
| 504 | << std::endl; |
| 505 | |
| 506 | std::cout << "=============================================================================================" << std::endl; |
| 507 | std::cout << "Measurement type (u,v,uv,none): " << static_cast<int>(measType_) << std::endl; |
| 508 | for (unsigned int iMeas = 0; iMeas < measurement_.size(); iMeas++) { |
| 509 | std::cout << "Measurement : "; measurement_.at(iMeas).Print(); |
| 510 | std::cout << "H Matrix : "; hMatrix_.Print(); |
| 511 | std::cout << "Measurement covariance : "; measCov_.at(iMeas).Print(); |
| 512 | std::cout << "Measurement residuals : "; getResidual(iMeas).getState().Print(); |
| 513 | std::cout << "Measurement down-weight: " << measDownWeight_.at(iMeas) << std::endl; |
| 514 | } |
| 515 | std::cout << "Jacobian: "; jacobian_.Print(); |
| 516 | std::cout << "Noise : "; noise_.Print(); |
| 517 | std::cout << "Backward covariance: "; bwdCov_.Print(); |
| 518 | std::cout << "Forward covariance : "; fwdCov_.Print(); |
| 519 | |
| 520 | std::cout << "=============================================================================================" << std::endl; |
| 521 | |
| 522 | } |
| 523 | |
| 524 | |
| 525 | } // end of namespace genfit |