fitters/src/DAF.cc

0001 /* Copyright 2008-2013, Technische Universitaet Muenchen, Ludwig-Maximilians-Universität München
0002    Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch & Tobias Schlüter
0003
0004    This file is part of GENFIT.
0005
0006    GENFIT is free software: you can redistribute it and/or modify
0007    it under the terms of the GNU Lesser General Public License as published
0008    by the Free Software Foundation, either version 3 of the License, or
0009    (at your option) any later version.
0010
0011    GENFIT is distributed in the hope that it will be useful,
0012    but WITHOUT ANY WARRANTY; without even the implied warranty of
0013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0014    GNU Lesser General Public License for more details.
0015
0016    You should have received a copy of the GNU Lesser General Public License
0017    along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
0018 */
0019
0020 #include "DAF.h"
0021 #include "Exception.h"
0022 #include "IO.h"
0023 #include "KalmanFitterInfo.h"
0024 #include "KalmanFitter.h"
0025 #include "KalmanFitterRefTrack.h"
0026 #include "KalmanFitStatus.h"
0027 #include "Tools.h"
0028 #include "Track.h"
0029 #include "TrackPoint.h"
0030
0031 #include <assert.h>
0032 #include <cmath>
0033
0034 //root stuff
0035 #include <TBuffer.h>
0036 #include <TMath.h>
0037 #include <Math/QuantFuncMathCore.h>
0038
0039
0040
0041 namespace genfit {
0042
0043 DAF::DAF(bool useRefKalman, double deltaPval, double deltaWeight)
0044   : AbsKalmanFitter(10, deltaPval), deltaWeight_(deltaWeight)
0045 {
0046   if (useRefKalman) {
0047     kalman_.reset(new KalmanFitterRefTrack());
0048     static_cast<KalmanFitterRefTrack*>(kalman_.get())->setRefitAll();
0049   }
0050   else
0051     kalman_.reset(new KalmanFitter());
0052
0053   kalman_->setMultipleMeasurementHandling(weightedAverage);
0054   kalman_->setMaxIterations(1);
0055
0056   setAnnealingScheme(100, 0.1, 5); // also sets maxIterations_
0057   setProbCut(0.001);
0058 }
0059
0060 DAF::DAF(AbsKalmanFitter* kalman, double deltaPval, double deltaWeight)
0061   : AbsKalmanFitter(10, deltaPval), deltaWeight_(deltaWeight)
0062 {
0063   kalman_.reset(kalman);
0064   kalman_->setMultipleMeasurementHandling(weightedAverage); // DAF makes no sense otherwise
0065   kalman_->setMaxIterations(1);
0066
0067   if (dynamic_cast<KalmanFitterRefTrack*>(kalman_.get()) != nullptr) {
0068     static_cast<KalmanFitterRefTrack*>(kalman_.get())->setRefitAll();
0069   }
0070
0071   setAnnealingScheme(100, 0.1, 5); // also sets maxIterations_
0072   setProbCut(0.01);
0073 }
0074
0075
0076 void DAF::processTrackWithRep(Track* tr, const AbsTrackRep* rep, bool resortHits) {
0077
0078   if (debugLvl_ > 0) {
0079     debugOut<<"DAF::processTrack //////////////////////////////////////////////////////////////// \n";
0080   }
0081
0082   KalmanFitStatus* status = 0;
0083   bool oneLastIter = false;
0084
0085   double lastPval = -1;
0086
0087   for(unsigned int iBeta = 0;; ++iBeta) {
0088
0089     if (debugLvl_ > 0) {
0090       debugOut<<"DAF::processTrack, trackRep  " << rep << ", iteration " << iBeta+1 << ", beta = " << betas_.at(iBeta) << "\n";
0091     }
0092
0093     kalman_->processTrackWithRep(tr, rep, resortHits);
0094
0095     status = static_cast<KalmanFitStatus*>(tr->getFitStatus(rep));
0096     status->setIsFittedWithDaf();
0097     status->setNumIterations(iBeta+1);
0098
0099
0100     // check break conditions
0101
0102     if (! status->isFitted()){
0103       if (debugLvl_ > 0) {
0104         debugOut << "DAF::Kalman could not fit!\n";
0105       }
0106       status->setIsFitted(false);
0107       break;
0108     }
0109
0110     if( oneLastIter == true){
0111       if (debugLvl_ > 0) {
0112         debugOut << "DAF::break after one last iteration\n";
0113       }
0114       status->setIsFitConvergedFully(status->getNFailedPoints() == 0);
0115       status->setIsFitConvergedPartially();
0116       break;
0117     }
0118
0119     if(iBeta >= maxIterations_-1){
0120       status->setIsFitConvergedFully(false);
0121       status->setIsFitConvergedPartially(false);
0122       if (debugLvl_ > 0) {
0123         debugOut << "DAF::number of max iterations reached!\n";
0124       }
0125       break;
0126     }
0127
0128
0129     // get and update weights
0130     bool converged(false);
0131     try{
0132       converged = calcWeights(tr, rep, betas_.at(iBeta));
0133       if (!converged && iBeta >= minIterations_-1 &&
0134           status->getBackwardPVal() != 0 && fabs(lastPval - status->getBackwardPVal()) < this->deltaPval_) {
0135         if (debugLvl_ > 0) {
0136           debugOut << "converged by Pval = " << status->getBackwardPVal() << " even though weights changed at iBeta = " << iBeta << std::endl;
0137         }
0138         converged = true;
0139       }
0140       lastPval = status->getBackwardPVal();
0141     } catch(Exception& e) {
0142       errorOut<<e.what();
0143       e.info();
0144       //errorOut << "calc weights failed" << std::endl;
0145       //mini_trk->getTrackRep(0)->setStatusFlag(1);
0146       status->setIsFitted(false);
0147       status->setIsFitConvergedFully(false);
0148       status->setIsFitConvergedPartially(false);
0149       break;
0150     }
0151
0152     // check if converged
0153     if (iBeta >= minIterations_-1 && converged) {
0154       if (debugLvl_ > 0) {
0155         debugOut << "DAF::convergence reached in iteration " << iBeta+1 << " -> Do one last iteration with updated weights.\n";
0156       }
0157       oneLastIter = true;
0158       status->setIsFitConvergedFully(status->getNFailedPoints() == 0);
0159       status->setIsFitConvergedPartially();
0160     }
0161
0162   } // end loop over betas
0163
0164
0165   if (status->getForwardPVal() == 0. &&
0166       status->getBackwardPVal() == 0.) {
0167     status->setIsFitConvergedFully(false);
0168     status->setIsFitConvergedPartially(false);
0169   }
0170
0171 }
0172
0173
0174 void DAF::setProbCut(const double prob_cut){
0175   for ( int i = 1; i < 7; ++i){
0176     addProbCut(prob_cut, i);
0177   }
0178 }
0179
0180 void DAF::addProbCut(const double prob_cut, const int measDim){
0181   if ( prob_cut > 1.0 || prob_cut < 0.0){
0182     Exception exc("DAF::addProbCut prob_cut is not between 0 and 1",__LINE__,__FILE__);
0183     exc.setFatal();
0184     throw exc;
0185   }
0186   if ( measDim < 1 || measDim > 6 ){
0187     Exception exc("DAF::addProbCut measDim must be in the range [1,6]",__LINE__,__FILE__);
0188     exc.setFatal();
0189     throw exc;
0190   }
0191   chi2Cuts_[measDim] = ROOT::Math::chisquared_quantile_c( prob_cut, measDim);
0192 }
0193
0194 void DAF::setAnnealingScheme(double bStart, double bFinal, unsigned int nSteps) {
0195   // The betas are calculated as a geometric series that takes nSteps
0196   // steps to go from bStart to bFinal.
0197   assert(bStart > bFinal);
0198   assert(bFinal > 1.E-10);
0199   assert(nSteps > 1);
0200
0201   minIterations_ = nSteps;
0202   maxIterations_ = nSteps + 4;
0203
0204   betas_.clear();
0205
0206   for (unsigned int i=0; i<nSteps; ++i) {
0207     betas_.push_back(bStart * pow(bFinal / bStart, i / (nSteps - 1.)));
0208   }
0209
0210   betas_.resize(maxIterations_,betas_.back()); //make sure main loop has a maximum of 10 iterations and also make sure the last beta value is used for if more iterations are needed then the ones set by the user.
0211
0212   /*for (unsigned int i=0; i<betas_.size(); ++i) {
0213     debugOut<< betas_.at(i) << ", ";
0214   }*/
0215 }
0216
0217
0218 bool DAF::calcWeights(Track* tr, const AbsTrackRep* rep, double beta) {
0219
0220   if (debugLvl_ > 0) {
0221     debugOut<<"DAF::calcWeights \n";
0222   }
0223
0224   bool converged(true);
0225   double maxAbsChange(0);
0226
0227   const std::vector< TrackPoint* >& trackPoints = tr->getPointsWithMeasurement();
0228   for (std::vector< TrackPoint* >::const_iterator tp = trackPoints.begin(); tp != trackPoints.end(); ++tp) {
0229     if (! (*tp)->hasFitterInfo(rep)) {
0230       continue;
0231     }
0232     AbsFitterInfo* fi = (*tp)->getFitterInfo(rep);
0233     if (dynamic_cast<KalmanFitterInfo*>(fi) == nullptr){
0234       Exception exc("DAF::getWeights ==> can only use KalmanFitterInfos",__LINE__,__FILE__);
0235       throw exc;
0236     }
0237     KalmanFitterInfo* kfi = static_cast<KalmanFitterInfo*>(fi);
0238
0239     if (kfi->areWeightsFixed()) {
0240       if (debugLvl_ > 0) {
0241         debugOut<<"weights are fixed, continue \n";
0242       }
0243       continue;
0244     }
0245
0246     unsigned int nMeas = kfi->getNumMeasurements();
0247
0248     std::vector<double> phi(nMeas, 0.);
0249     double phi_sum = 0;
0250     double phi_cut = 0;
0251     for(unsigned int j=0; j<nMeas; j++) {
0252
0253       try{
0254         const MeasurementOnPlane& residual = kfi->getResidual(j, true, true);
0255         const TVectorD& resid(residual.getState());
0256         TMatrixDSym Vinv(residual.getCov());
0257         double detV;
0258         tools::invertMatrix(Vinv, &detV); // can throw an Exception
0259         int hitDim = resid.GetNrows();
0260         // Needed for normalization, special cases for the two common cases,
0261         // shouldn't matter, but the original code made some efforts to make
0262         // this calculation faster, and it's not complex ...
0263         double twoPiN = 2.*M_PI;
0264         if (hitDim == 2)
0265           twoPiN *= twoPiN;
0266         else if (hitDim > 2)
0267           twoPiN = pow(twoPiN, hitDim);
0268
0269         double chi2 = Vinv.Similarity(resid);
0270         if (debugLvl_ > 1) {
0271           debugOut<<"chi2 = " << chi2 << "\n";
0272         }
0273
0274     // The common factor beta is eliminated.
0275         double norm = 1./sqrt(twoPiN * detV);
0276
0277         phi[j] = norm*exp(-0.5*chi2/beta);
0278     phi_sum += phi[j];
0279         //errorOut << "hitDim " << hitDim << " fchi2Cuts[hitDim] " << fchi2Cuts[hitDim] << std::endl;
0280         double cutVal = chi2Cuts_[hitDim];
0281         assert(cutVal>1.E-6);
0282         //the following assumes that in the competing hits could have different V otherwise calculation could be simplified
0283         phi_cut += norm*exp(-0.5*cutVal/beta);
0284       }
0285       catch(Exception& e) {
0286         errorOut << e.what();
0287         e.info();
0288       }
0289     }
0290
0291     for(unsigned int j=0; j<nMeas; j++) {
0292       double weight = phi[j]/(phi_sum+phi_cut);
0293       //debugOut << phi_sum << " " << phi_cut << " " << weight << std::endl;
0294
0295       // check convergence
0296       double absChange(fabs(weight - kfi->getMeasurementOnPlane(j)->getWeight()));
0297       if (converged && absChange > deltaWeight_) {
0298         converged = false;
0299         if (absChange > maxAbsChange)
0300           maxAbsChange = absChange;
0301       }
0302
0303       if (debugLvl_ > 0) {
0304         if (debugLvl_ > 1 || absChange > deltaWeight_) {
0305           debugOut<<"\t old weight: " << kfi->getMeasurementOnPlane(j)->getWeight();
0306           debugOut<<"\t new weight: " << weight;
0307         }
0308       }
0309
0310       kfi->getMeasurementOnPlane(j)->setWeight(weight);
0311     }
0312   }
0313
0314   if (debugLvl_ > 0) {
0315     debugOut << "\t  ";
0316     debugOut << "max abs weight change = " << maxAbsChange << "\n";
0317   }
0318
0319   return converged;
0320 }
0321
0322
0323 // Customized from generated Streamer.
0324 void DAF::Streamer(TBuffer &R__b)
0325 {
0326    // Stream an object of class genfit::DAF.
0327
0328    //This works around a msvc bug and should be harmless on other platforms
0329    typedef ::genfit::DAF thisClass;
0330    UInt_t R__s, R__c;
0331    if (R__b.IsReading()) {
0332       Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v) { }
0333       //This works around a msvc bug and should be harmless on other platforms
0334       typedef genfit::AbsKalmanFitter baseClass0;
0335       baseClass0::Streamer(R__b);
0336       R__b >> deltaWeight_;
0337       // weights_ are only of intermediate use -> not saved
0338       {
0339          std::vector<double> &R__stl =  betas_;
0340          R__stl.clear();
0341          int R__i, R__n;
0342          R__b >> R__n;
0343          R__stl.reserve(R__n);
0344          for (R__i = 0; R__i < R__n; R__i++) {
0345             double R__t;
0346             R__b >> R__t;
0347             R__stl.push_back(R__t);
0348          }
0349       }
0350       if (R__v == 1) {
0351      // Old versions kept the chi2Cuts_ in a map.
0352          // We ignore non-sensical dimensionalities when reading it again.
0353          int R__i, R__n;
0354          R__b >> R__n;
0355          for (R__i = 0; R__i < R__n; R__i++) {
0356         memset(chi2Cuts_, 0, sizeof(chi2Cuts_));
0357             int R__t;
0358             R__b >> R__t;
0359             double R__t2;
0360             R__b >> R__t2;
0361         if (R__t >= 1 && R__t <= 6)
0362           chi2Cuts_[R__t] = R__t2;
0363      }
0364       } else {
0365     char n_chi2Cuts;  // should be six
0366     R__b >> n_chi2Cuts;
0367     assert(n_chi2Cuts == 6);  // Cannot be different as long as sanity prevails.
0368     chi2Cuts_[0] = 0; // nonsensical.
0369     R__b.ReadFastArray(&chi2Cuts_[1], n_chi2Cuts);
0370       }
0371       AbsKalmanFitter *p;
0372       R__b >> p;
0373       kalman_.reset(p);
0374       R__b.CheckByteCount(R__s, R__c, thisClass::IsA());
0375    } else {
0376       R__c = R__b.WriteVersion(thisClass::IsA(), kTRUE);
0377       //This works around a msvc bug and should be harmless on other platforms
0378       typedef genfit::AbsKalmanFitter baseClass0;
0379       baseClass0::Streamer(R__b);
0380       R__b << deltaWeight_;
0381       {
0382          std::vector<double> &R__stl =  betas_;
0383          int R__n=int(R__stl.size());
0384          R__b << R__n;
0385          if(R__n) {
0386             std::vector<double>::iterator R__k;
0387             for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
0388             R__b << (*R__k);
0389             }
0390          }
0391       }
0392       {
0393     R__b << (char)6;  // Number of chi2Cuts_
0394     R__b.WriteFastArray(&chi2Cuts_[1], 6);
0395       }
0396       R__b << kalman_.get();
0397       R__b.SetByteCount(R__c, kTRUE);
0398    }
0399 }
0400
0401 } /* End of namespace genfit */