genfit::DAF Class Reference

Determinstic Annealing Filter (DAF) implementation. More...

#include <DAF.h>

Inheritance diagram for genfit::DAF:

Collaboration diagram for genfit::DAF:

Public Member Functions
	DAF (bool useRefKalman=true, double deltaWeight=1e-3, double deltaPval=1e-3)
	Create DAF. Per default, use KalmanFitterRefTrack as fitter.
	DAF (AbsKalmanFitter *kalman, double deltaWeight=1e-3, double deltaPval=1e-3)
	Create DAF. Use the provided AbsKalmanFitter as fitter.
	~DAF ()
void	processTrackWithRep (Track *tr, const AbsTrackRep *rep, bool resortHits=false)
	Process a track using the DAF.
void	setProbCut (const double prob_cut)
	Set the probability cut for the weight calculation for the hits.
void	addProbCut (const double prob_cut, const int measDim)
	Set the probability cut for the weight calculation for the hits for a specific measurement dimensionality.
void	setBetas (double b1, double b2=-1, double b3=-1., double b4=-1., double b5=-1., double b6=-1., double b7=-1., double b8=-1., double b9=-1., double b10=-1.)
	Configure the annealing scheme.
void	setAnnealingScheme (double bStart, double bFinal, unsigned int nSteps)
	Configure the annealing scheme.
void	setMaxIterations (unsigned int n)
	Set the maximum number of iterations.
void	setConvergenceDeltaWeight (double delta)
	If all weights change less than delta between two iterations, the fit is regarded as converged.
AbsKalmanFitter *	getKalman () const
virtual void	setMaxFailedHits (int val)
virtual void	setDebugLvl (unsigned int lvl=1)
Public Member Functions inherited from genfit::AbsKalmanFitter
	AbsKalmanFitter (unsigned int maxIterations=4, double deltaPval=1e-3, double blowUpFactor=1e3)
virtual	~AbsKalmanFitter ()
void	getChiSquNdf (const Track *tr, const AbsTrackRep *rep, double &bChi2, double &fChi2, double &bNdf, double &fNdf) const
double	getChiSqu (const Track *tr, const AbsTrackRep *rep, int direction=-1) const
double	getNdf (const Track *tr, const AbsTrackRep *rep, int direction=-1) const
double	getRedChiSqu (const Track *tr, const AbsTrackRep *rep, int direction=-1) const
double	getPVal (const Track *tr, const AbsTrackRep *rep, int direction=-1) const
eMultipleMeasurementHandling	getMultipleMeasurementHandling () const
virtual void	setMinIterations (unsigned int n)
	Set the minimum number of iterations.
void	setDeltaPval (double deltaPval)
	Set Convergence criterion.
void	setRelChi2Change (double relChi2Change)
void	setMultipleMeasurementHandling (eMultipleMeasurementHandling mmh)
	How should multiple measurements be handled?
bool	isTrackPrepared (const Track *tr, const AbsTrackRep *rep) const
bool	isTrackFitted (const Track *tr, const AbsTrackRep *rep) const
bool	canIgnoreWeights () const
	returns if the fitter can ignore the weights and handle the MeasurementOnPlanes as if they had weight 1.
Public Member Functions inherited from genfit::AbsFitter
	AbsFitter ()
virtual	~AbsFitter ()
void	processTrack (Track *, bool resortHits=true)

Private Member Functions
	DAF (const DAF &)
DAF &	operator= (genfit::DAF const &)
bool	calcWeights (Track *trk, const AbsTrackRep *rep, double beta)
	Calculate and set the weights for the next fitting pass. Return if convergence is met. The convergence criterium is the largest absolute change of all weights.

Private Attributes
double	deltaWeight_
std::vector< double >	betas_
std::map< int, double >	chi2Cuts_
boost::scoped_ptr< AbsKalmanFitter >	kalman_

Additional Inherited Members
Protected Member Functions inherited from genfit::AbsKalmanFitter
const std::vector< MeasurementOnPlane * >	getMeasurements (const KalmanFitterInfo *fi, const TrackPoint *tp, int direction) const
	get the measurementsOnPlane taking the multipleMeasurementHandling_ into account
Protected Attributes inherited from genfit::AbsKalmanFitter
unsigned int	minIterations_
	Minimum number of iterations to attempt. Forward and backward are counted as one iteration.
unsigned int	maxIterations_
	Maximum number of iterations to attempt. Forward and backward are counted as one iteration.
double	deltaPval_
	Convergence criterion.
double	relChi2Change_
double	blowUpFactor_
	Blow up the covariance of the forward (backward) fit by this factor before seeding the backward (forward) fit.
eMultipleMeasurementHandling	multipleMeasurementHandling_
	How to handle if there are multiple MeasurementsOnPlane.
int	maxFailedHits_
Protected Attributes inherited from genfit::AbsFitter
unsigned int	debugLvl_

Detailed Description

Determinstic Annealing Filter (DAF) implementation.

Author

Christian Höppner (Technische Universität München, original author)

Karl Bicker (Technische Universität München)

The DAF is an iterative Kalman filter with annealing. It is capable of fitting tracks which are contaminated with noise hits. The algorithm is taken from the references R. Fruehwirth & A. Strandlie, Computer Physics Communications 120 (1999) 197-214 and CERN thesis: Dissertation by Matthias Winkler.

The weights which were assigned to the hits by the DAF are accessible in the MeasurementOnPlane objects in the KalmanFitterInfo objects.

Definition at line 48 of file DAF.h.

Constructor & Destructor Documentation

DAF() [1/3]

genfit::DAF::DAF ( const DAF &  )

private

DAF() [2/3]

genfit::DAF::DAF	(	bool	useRefKalman = true,
		double	deltaWeight = 1e-3,
		double	deltaPval = 1e-3
	)

Create DAF. Per default, use KalmanFitterRefTrack as fitter.

Parameters

useRefKalman

If false, use KalmanFitter as fitter.

Definition at line 42 of file DAF.cc.

  : AbsKalmanFitter(10, deltaPval), deltaWeight_(deltaWeight)
{
  if (useRefKalman) {
    kalman_.reset(new KalmanFitterRefTrack());
    static_cast<KalmanFitterRefTrack*>(kalman_.get())->setRefitAll();
  }
  else
    kalman_.reset(new KalmanFitter());
 
  kalman_->setMultipleMeasurementHandling(weightedAverage);
  kalman_->setMaxIterations(1);
 
  setAnnealingScheme(100, 0.1, 5); // also sets maxIterations_
  setProbCut(0.01);
}

DAF() [3/3]

genfit::DAF::DAF	(	AbsKalmanFitter *	kalman,
		double	deltaWeight = 1e-3,
		double	deltaPval = 1e-3
	)

Create DAF. Use the provided AbsKalmanFitter as fitter.

Definition at line 59 of file DAF.cc.

  : AbsKalmanFitter(10, deltaPval), deltaWeight_(deltaWeight)
{
  kalman_.reset(kalman);
  kalman_->setMultipleMeasurementHandling(weightedAverage); // DAF makes no sense otherwise
  kalman_->setMaxIterations(1);
 
  if (dynamic_cast<KalmanFitterRefTrack*>(kalman_.get()) != NULL) {
    static_cast<KalmanFitterRefTrack*>(kalman_.get())->setRefitAll();
  }
 
  setAnnealingScheme(100, 0.1, 5); // also sets maxIterations_
  setProbCut(0.01);
}

~DAF()

genfit::DAF::~DAF ( )

inline

Definition at line 67 of file DAF.h.

{};

Member Function Documentation

addProbCut()

void genfit::DAF::addProbCut	(	const double	prob_cut,
		const int	measDim
	)

Set the probability cut for the weight calculation for the hits for a specific measurement dimensionality.

Definition at line 179 of file DAF.cc.

                                                            {
  if ( prob_cut > 1.0 || prob_cut < 0.0){
    Exception exc("DAF::addProbCut prob_cut is not between 0 and 1",__LINE__,__FILE__);
    exc.setFatal();
    throw exc;
  }
  if ( measDim < 1){
    Exception exc("DAF::addProbCut measDim must be > 0",__LINE__,__FILE__);
    exc.setFatal();
    throw exc;
  }
  chi2Cuts_[measDim] = ROOT::Math::chisquared_quantile_c( prob_cut, measDim);
}

calcWeights()

bool genfit::DAF::calcWeights ( Track *  trk, const AbsTrackRep *  rep, double  beta  )

private

Calculate and set the weights for the next fitting pass. Return if convergence is met. The convergence criterium is the largest absolute change of all weights.

Definition at line 254 of file DAF.cc.

                                                                    {
 
  if (debugLvl_ > 0) {
    std::cout<<"DAF::calcWeights \n";
  }
 
  bool converged(true);
  double maxAbsChange(0);
 
  const std::vector< TrackPoint* >& trackPoints = tr->getPointsWithMeasurement();
  for (std::vector< TrackPoint* >::const_iterator tp = trackPoints.begin(); tp != trackPoints.end(); ++tp) {
    if (! (*tp)->hasFitterInfo(rep)) {
      continue;
    }
    AbsFitterInfo* fi = (*tp)->getFitterInfo(rep);
    if (dynamic_cast<KalmanFitterInfo*>(fi) == NULL){
      Exception exc("DAF::getWeights ==> can only use KalmanFitterInfos",__LINE__,__FILE__);
      throw exc;
    }
    KalmanFitterInfo* kfi = static_cast<KalmanFitterInfo*>(fi);
 
    if (kfi->areWeightsFixed()) {
      if (debugLvl_ > 0) {
        std::cout<<"weights are fixed, continue \n";
      }
      continue;
    }
 
    unsigned int nMeas = kfi->getNumMeasurements();
 
    std::vector<double> phi(nMeas, 0.);
    double phi_sum = 0;
    double phi_cut = 0;
    for(unsigned int j=0; j<nMeas; j++) {
 
      try{
        const MeasurementOnPlane& residual = kfi->getResidual(j, true, true);
        const TVectorD& resid(residual.getState());
        TMatrixDSym Vinv(residual.getCov());
        double detV;
        tools::invertMatrix(Vinv, &detV); // can throw an Exception
        int hitDim = resid.GetNrows();
        // Needed for normalization, special cases for the two common cases,
        // shouldn't matter, but the original code made some efforts to make
        // this calculation faster, and it's not complex ...
        double twoPiN = 2.*M_PI;
        if (hitDim == 2)
          twoPiN *= twoPiN;
        else if (hitDim > 2)
          twoPiN = pow(twoPiN, hitDim);
 
        double chi2 = Vinv.Similarity(resid);
        if (debugLvl_ > 1) {
          std::cout<<"chi2 = " << chi2 << "\n";
        }
 
        double norm = 1./sqrt(twoPiN * detV);
 
        phi[j] = norm*exp(-0.5*chi2/beta);
        phi_sum += phi[j];
        //std::cerr << "hitDim " << hitDim << " fchi2Cuts[hitDim] " << fchi2Cuts[hitDim] << std::endl;
        double cutVal = chi2Cuts_[hitDim];
        assert(cutVal>1.E-6);
        //the following assumes that in the competing hits could have different V otherwise calculation could be simplified
        phi_cut += norm*exp(-0.5*cutVal/beta);
      }
      catch(Exception& e) {
        std::cerr << e.what();
        e.info();
      }
    }
 
    for(unsigned int j=0; j<nMeas; j++) {
      double weight = phi[j]/(phi_sum+phi_cut);
 
      // check convergence
      double absChange(fabs(weight - kfi->getMeasurementOnPlane(j)->getWeight()));
      if (converged && absChange > deltaWeight_) {
        converged = false;
        if (absChange > maxAbsChange)
          maxAbsChange = absChange;
      }
 
      if (debugLvl_ > 0) {
        if (debugLvl_ > 1 || absChange > deltaWeight_) {
          std::cout<<"\t old weight: " << kfi->getMeasurementOnPlane(j)->getWeight();
          std::cout<<"\t new weight: " << weight;
        }
      }
 
      kfi->getMeasurementOnPlane(j)->setWeight(weight);
    }
  }
 
  if (debugLvl_ > 0) {
    std::cout << "\t  ";
    std::cout << "max abs weight change = " << maxAbsChange << "\n";
  }
 
  return converged;
}

getKalman()

AbsKalmanFitter * genfit::DAF::getKalman ( ) const

inline

Definition at line 103 of file DAF.h.

{return kalman_.get();}

operator=()

DAF & genfit::DAF::operator= ( genfit::DAF const &  )

private

processTrackWithRep()

void genfit::DAF::processTrackWithRep ( Track *  tr, const AbsTrackRep *  rep, bool  resortHits = false  )

virtual

Process a track using the DAF.

Implements genfit::AbsFitter.

Definition at line 75 of file DAF.cc.

                                                                                {
 
  if (debugLvl_ > 0) {
    std::cout<<"DAF::processTrack //////////////////////////////////////////////////////////////// \n";
  }
 
  KalmanFitStatus* status = 0;
  bool oneLastIter = false;
 
  double lastPval = -1;
 
  for(unsigned int iBeta = 0;; ++iBeta) {
 
    if (debugLvl_ > 0) {
      std::cout<<"DAF::processTrack, trackRep  " << rep << ", iteration " << iBeta+1 << ", beta = " << betas_.at(iBeta) << "\n";
    }
 
    kalman_->processTrackWithRep(tr, rep, resortHits);
 
    status = static_cast<KalmanFitStatus*>(tr->getFitStatus(rep));
    status->setIsFittedWithDaf();
    status->setNumIterations(iBeta+1);
 
 
    // check break conditions
 
    if (! status->isFitted()){
      if (debugLvl_ > 0) {
        std::cout << "DAF::Kalman could not fit!\n";
      }
      status->setIsFitted(false);
      break;
    }
 
    if( oneLastIter == true){
      if (debugLvl_ > 0) {
        std::cout << "DAF::break after one last iteration\n";
      }
      status->setIsFitConvergedFully(status->getNFailedPoints() == 0);
      status->setIsFitConvergedPartially();
      break;
    }
 
    if(iBeta >= maxIterations_-1){
      status->setIsFitConvergedFully(false);
      status->setIsFitConvergedPartially(false);
      if (debugLvl_ > 0) {
        std::cout << "DAF::number of max iterations reached!\n";
      }
      break;
    }
 
 
    // get and update weights
    bool converged(false);
    try{
      converged = calcWeights(tr, rep, betas_.at(iBeta));
      if (!converged && iBeta >= minIterations_-1 &&
          status->getBackwardPVal() != 0 && fabs(lastPval - status->getBackwardPVal()) < this->deltaPval_) {
        if (debugLvl_ > 0) {
          std::cout << "converged by Pval = " << status->getBackwardPVal() << " even though weights changed at iBeta = " << iBeta << std::endl;
        }
        converged = true;
      }
      lastPval = status->getBackwardPVal();
    } catch(Exception& e) {
      std::cerr<<e.what();
      e.info();
      //std::cerr << "calc weights failed" << std::endl;
      //mini_trk->getTrackRep(0)->setStatusFlag(1);
      status->setIsFitted(false);
      status->setIsFitConvergedFully(false);
      status->setIsFitConvergedPartially(false);
      break;
    }
 
    // check if converged
    if (iBeta >= minIterations_-1 && converged) {
      if (debugLvl_ > 0) {
        std::cout << "DAF::convergence reached in iteration " << iBeta+1 << " -> Do one last iteration with updated weights.\n";
      }
      oneLastIter = true;
      status->setIsFitConvergedFully(status->getNFailedPoints() == 0);
      status->setIsFitConvergedPartially();
    }
 
  } // end loop over betas
 
 
  if (status->getForwardPVal() == 0. &&
      status->getBackwardPVal() == 0.) {
    status->setIsFitConvergedFully(false);
    status->setIsFitConvergedPartially(false);
  }
 
}

setAnnealingScheme()

void genfit::DAF::setAnnealingScheme	(	double	bStart,
		double	bFinal,
		unsigned int	nSteps
	)

Configure the annealing scheme.

Set a start and end temperature and the number of steps. A logarithmic sequence of temperatures will be calculated. Also sets minIterations_ and maxIterations_.

Definition at line 230 of file DAF.cc.

                                                                              {
  // The betas are calculated as a geometric series that takes nSteps
  // steps to go from bStart to bFinal.
  assert(bStart > bFinal);
  assert(bFinal > 1.E-10);
  assert(nSteps > 1);
 
  minIterations_ = nSteps;
  maxIterations_ = nSteps + 4;
 
  betas_.clear();
 
  for (unsigned int i=0; i<nSteps; ++i) {
    betas_.push_back(bStart * pow(bFinal / bStart, i / (nSteps - 1.)));
  }
 
  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.
 
  /*for (unsigned int i=0; i<betas_.size(); ++i) {
    std::cout<< betas_.at(i) << ", ";
  }*/
}

setBetas()

void genfit::DAF::setBetas	(	double	b1,
		double	b2 = -1,
		double	b3 = -1.,
		double	b4 = -1.,
		double	b5 = -1.,
		double	b6 = -1.,
		double	b7 = -1.,
		double	b8 = -1.,
		double	b9 = -1.,
		double	b10 = -1.
	)

Configure the annealing scheme.

In the current implementation you need to provide at least one temperature and not more then ten temperatures. Also sets minIterations_ and maxIterations_.

Definition at line 194 of file DAF.cc.

                                                                                                                      {
  betas_.clear();
  assert(b1>0);betas_.push_back(b1);
  if(b2>0){
    assert(b2<=b1);betas_.push_back(b2);
    if(b3>=0.) {
      assert(b3<=b2);betas_.push_back(b3);
      if(b4>=0.) {
        assert(b4<=b3);betas_.push_back(b4);
        if(b5>=0.) {
          assert(b5<=b4);betas_.push_back(b5);
          if(b6>=0.) {
            assert(b6<=b5);betas_.push_back(b6);
            if(b7>=0.) {
              assert(b7<=b6);betas_.push_back(b7);
              if(b8>=0.) {
                assert(b8<=b7);betas_.push_back(b8);
                if(b9>=0.) {
                  assert(b9<=b8);betas_.push_back(b9);
                  if(b10>=0.) {
                    assert(b10<=b9);betas_.push_back(b10);
                  }
                }
              }
            }
          }
        }
      }
    }
  }
  minIterations_ = betas_.size();
  maxIterations_ = betas_.size() + 4;
  betas_.resize(maxIterations_,betas_.back()); //make sure main loop has a maximum of maxIterations_ and also make sure the last beta value is used for if more iterations are needed then the ones set by the user.
}

setConvergenceDeltaWeight()

void genfit::DAF::setConvergenceDeltaWeight ( double  delta )

inline

If all weights change less than delta between two iterations, the fit is regarded as converged.

Definition at line 101 of file DAF.h.

{deltaWeight_ = delta;}

setDebugLvl()

virtual void genfit::DAF::setDebugLvl ( unsigned int  lvl = 1 )

inlinevirtual

Reimplemented from genfit::AbsFitter.

Definition at line 107 of file DAF.h.

{AbsFitter::setDebugLvl(lvl); if (lvl > 1) getKalman()->setDebugLvl(lvl-1);}

setMaxFailedHits()

virtual void genfit::DAF::setMaxFailedHits ( int  val )

inlinevirtual

Reimplemented from genfit::AbsKalmanFitter.

Definition at line 105 of file DAF.h.

{getKalman()->setMaxFailedHits(val);}

setMaxIterations()

void genfit::DAF::setMaxIterations ( unsigned int  n )

inlinevirtual

Set the maximum number of iterations.

Reimplemented from genfit::AbsKalmanFitter.

Definition at line 98 of file DAF.h.

{maxIterations_ = n; betas_.resize(maxIterations_,betas_.back());}

setProbCut()

void genfit::DAF::setProbCut

(

const double

prob_cut

)

Set the probability cut for the weight calculation for the hits.

By default the cut values for measurements of dimensionality from 1 to 5 are calculated. If you what to have cut values for an arbitrary measurement dimensionality use addProbCut(double prob_cut, int maxDim);

Definition at line 173 of file DAF.cc.

                                         {
  for ( int i = 1; i != 6; ++i){
    addProbCut(prob_cut, i);
  }
}

Member Data Documentation

betas_

std::vector<double> genfit::DAF::betas_

private

Definition at line 119 of file DAF.h.

chi2Cuts_

std::map<int,double> genfit::DAF::chi2Cuts_

private

Definition at line 120 of file DAF.h.

deltaWeight_

double genfit::DAF::deltaWeight_

private

Definition at line 118 of file DAF.h.

kalman_

boost::scoped_ptr<AbsKalmanFitter> genfit::DAF::kalman_

private

Definition at line 122 of file DAF.h.

---

The documentation for this class was generated from the following files:

• genfit/fitters/include/DAF.h
• genfit/fitters/src/DAF.cc