GFKalman Class Reference

Generic Kalman Filter implementation. More...

#include <GFKalman.h>

Inheritance diagram for GFKalman:

Public Member Functions
	GFKalman ()
	~GFKalman ()
void	operator() (GFTrack *track)
	Operator for use with STL. More...
void	operator() (std::pair< int, GFTrack * > tr)
	Operator for use with STL. More...
void	setLazy (Int_t)
	Switch lazy error handling. More...
void	setNumIterations (Int_t i)
	Set number of iterations for Kalman Filter. More...
void	processTrack (GFTrack *trk)
	Performs fit on a GFTrack. More...
void	fittingPass (GFTrack *, int dir)
	Performs fit on a GFTrack beginning with the current hit. More...
double	getChi2Hit (GFAbsRecoHit *, GFAbsTrackRep *)
	Calculates chi2 of a given hit with respect to a given track representation. More...
void	setInitialDirection (int d)
	Sets the inital direction of the track fit (1 for inner to outer, or -1 for outer to inner). The standard is 1 and is set in the ctor. More...
void	setBlowUpFactor (double f)
	Set the blowup factor (see blowUpCovs() ) More...

Private Member Functions
void	processHit (GFTrack *, int, int, int)
	One Kalman step. More...
void	switchDirection (GFTrack *trk)
	Used to switch between forward and backward filtering. More...
TMatrixT< double >	calcGain (const TMatrixT< double > &cov, const TMatrixT< double > &HitCov, const TMatrixT< double > &H)
	Calculate Kalman Gain. More...
double	chi2Increment (const TMatrixT< double > &r, const TMatrixT< double > &H, const TMatrixT< double > &cov, const TMatrixT< double > &V)
	this returns the reduced chi2 increment for a hit More...
void	blowUpCovs (GFTrack *trk)
	this is needed to blow up the covariance matrix before a fitting pass drops off-diagonal elements and blows up diagonal by blowUpFactor More...

Private Attributes
int	fInitialDirection
Int_t	fNumIt
double	fBlowUpFactor
bool	fSmooth

Detailed Description

Generic Kalman Filter implementation.

Author

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

Sebastian Neubert (Technische Universität München, original author)

The Kalman Filter operates on genfit GFTrack objects. It is an implementation of the Kalman Filter algebra that uses the genfit interface classes GFAbsRecoHit and GFAbsTrackRep in order to be independent from the specific detector geometry and the details of the track parameterization / track extraoplation engine.

The Kalman Filter can use hits from several detectors in a single fit to estimate the parameters of several track representations in parallel.

Definition at line 50 of file GFKalman.h.

Constructor & Destructor Documentation

GFKalman::GFKalman

(

)

Definition at line 36 of file GFKalman.cxx.

:fInitialDirection(1),fNumIt(3),fBlowUpFactor(500.){;}

GFKalman::~GFKalman

(

)

Definition at line 38 of file GFKalman.cxx.

{;}

Member Function Documentation

void GFKalman::blowUpCovs ( GFTrack *  trk )

private

this is needed to blow up the covariance matrix before a fitting pass drops off-diagonal elements and blows up diagonal by blowUpFactor

Definition at line 128 of file GFKalman.cxx.

References GFTrack::blowUpCovs(), and fBlowUpFactor.

Referenced by processTrack().

                                     {
  trk->blowUpCovs(fBlowUpFactor);
}

TMatrixT< double > GFKalman::calcGain ( const TMatrixT< double > &  cov, const TMatrixT< double > &  HitCov, const TMatrixT< double > &  H  )

private

Calculate Kalman Gain.

Definition at line 327 of file GFKalman.cxx.

References GFTools::invertMatrix().

Referenced by processHit().

                                                                   {

  // calculate covsum (V + HCH^T)
  TMatrixT<double> covsum1(cov,TMatrixT<double>::kMultTranspose,H);
  TMatrixT<double> covsum(H,TMatrixT<double>::kMult,covsum1);

  covsum+=HitCov;
  
  // invert
  TMatrixT<double> covsumInv;
  GFTools::invertMatrix(covsum,covsumInv);
  
  // calculate gain
  TMatrixT<double> gain1(H,TMatrixT<double>::kTransposeMult,covsumInv);
  TMatrixT<double> gain(cov,TMatrixT<double>::kMult,gain1);

  return gain;
}

double GFKalman::chi2Increment ( const TMatrixT< double > &  r, const TMatrixT< double > &  H, const TMatrixT< double > &  cov, const TMatrixT< double > &  V  )

private

this returns the reduced chi2 increment for a hit

Definition at line 178 of file GFKalman.cxx.

References GFTools::invertMatrix(), r, R, GFException::setFatal(), and GFException::setMatrices().

Referenced by getChi2Hit(), and processHit().

                                                                                   {

  // residuals covariances:R=(V - HCH^T)
  TMatrixT<double> R(V);
  TMatrixT<double> covsum1(cov,TMatrixT<double>::kMultTranspose,H);
  TMatrixT<double> covsum(H,TMatrixT<double>::kMult,covsum1);

  R-=covsum;

  // chisq= r^TR^(-1)r
  TMatrixT<double> Rinv;
  GFTools::invertMatrix(R,Rinv);
  
  
  TMatrixT<double> residTranspose(r);
  residTranspose.T();
  TMatrixT<double> chisq=residTranspose*(Rinv*r);
  assert(chisq.GetNoElements()==1);

  if(TMath::IsNaN(chisq[0][0])){
        GFException exc("chi2 is nan",__LINE__,__FILE__);
        exc.setFatal();
        std::vector< TMatrixT<double> > matrices;
        matrices.push_back(r);
        matrices.push_back(V);
        matrices.push_back(R);
        matrices.push_back(cov);
        exc.setMatrices("r, V, R, cov",matrices);
    throw exc;
  }

  return chisq[0][0];
}

void GFKalman::fittingPass	(	GFTrack *	trk,
		int	dir
	)

Performs fit on a GFTrack beginning with the current hit.

Definition at line 133 of file GFKalman.cxx.

References GFTrack::addFailedHit(), GFBookkeeping::clearFailedHits(), GFTrack::getBK(), GFTrack::getNextHitToFit(), GFTrack::getNumHits(), GFTrack::getNumReps(), GFAbsTrackRep::getStatusFlag(), GFTrack::getTrackRep(), GFException::info(), GFException::isFatal(), nhits, processHit(), GFAbsTrackRep::setChiSqu(), GFAbsTrackRep::setNDF(), GFTrack::setNextHitToFit(), GFAbsTrackRep::setStatusFlag(), and GFException::what().

Referenced by processTrack().

                                                {
  //loop over hits
  unsigned int nhits=trk->getNumHits();
  unsigned int starthit=trk->getNextHitToFit();

  int nreps=trk->getNumReps();
  int ihit=(int)starthit;
  
  for(int irep=0; irep<nreps; ++irep) {
    GFAbsTrackRep* arep=trk->getTrackRep(irep);
    if(arep->getStatusFlag()==0) {
      //clear chi2 sum and ndf sum in track reps
      arep->setChiSqu(0.);
      arep->setNDF(0);
      //clear failedHits and outliers
      trk->getBK(irep)->clearFailedHits();
    }
  }

  while((ihit<(int)nhits && direction==1) || (ihit>-1 && direction==-1)){
    //    GFAbsRecoHit* ahit=trk->getHit(ihit);
    // loop over reps
    for(int irep=0; irep<nreps; ++irep){
          GFAbsTrackRep* arep=trk->getTrackRep(irep);
          if(arep->getStatusFlag()==0) { 
            try {
              processHit(trk,ihit,irep,direction);
            }
            catch(GFException& e) {
              trk->addFailedHit(irep,ihit);
              std::cerr << e.what();
              e.info();
              if(e.isFatal()) {
                arep->setStatusFlag(1);
                continue; // go to next rep immediately
              }
            }   
          }
    }// end loop over reps
    ihit+=direction;
  }// end loop over hits
  trk->setNextHitToFit(ihit-2*direction);
  //trk->printGFBookkeeping();
}

double GFKalman::getChi2Hit	(	GFAbsRecoHit *	hit,
		GFAbsTrackRep *	rep
	)

Calculates chi2 of a given hit with respect to a given track representation.

Definition at line 215 of file GFKalman.cxx.

References chi2Increment(), GFAbsTrackRep::extrapolate(), GFAbsRecoHit::getDetPlane(), GFAbsTrackRep::getDim(), GFAbsRecoHit::getHitCov(), GFAbsRecoHit::getHMatrix(), r, and GFAbsRecoHit::residualVector().

{
  // get prototypes for matrices
  int repDim=rep->getDim();
  TMatrixT<double> state(repDim,1);
  TMatrixT<double> cov(repDim,repDim);;
  GFDetPlane pl=hit->getDetPlane(rep);
  rep->extrapolate(pl,state,cov);


  TMatrixT<double> H = hit->getHMatrix(rep);
  // get hit covariances  
  TMatrixT<double> V=hit->getHitCov(pl);

  TMatrixT<double> r=hit->residualVector(rep,state,pl);
  assert(r.GetNrows()>0);

  //this is where chi2 is calculated
  double chi2 = chi2Increment(r,H,cov,V);

  return chi2/r.GetNrows();
}

void GFKalman::operator() ( GFTrack *  track )

inline

Operator for use with STL.

This operator allows to use the std::foreach algorithm with an STL container o GFTrack* objects.

Definition at line 65 of file GFKalman.h.

References processTrack().

{processTrack(track);}

void GFKalman::operator() ( std::pair< int, GFTrack * >  tr )

inline

Operator for use with STL.

This operator allows to use the std::foreach algorithm with an STL container o GFTrack* objects.

Definition at line 72 of file GFKalman.h.

References processTrack().

{processTrack(tr.second);}

void GFKalman::processHit ( GFTrack *  tr, int  ihit, int  irep, int  direction  )

private

One Kalman step.

Performs

• Extrapolation to detector plane of the hit
• Calculation of residual and Kalman Gain
• Update of track representation state and chi2

Definition at line 240 of file GFKalman.cxx.

References GFAbsTrackRep::addChiSqu(), GFAbsTrackRep::addNDF(), calcGain(), chi2Increment(), COVEXC, GFAbsTrackRep::extrapolate(), fSmooth, GFAbsTrackRep::getAuxInfo(), GFTrack::getBK(), GFAbsTrackRep::getCov(), GFAbsRecoHit::getDetPlane(), GFAbsTrackRep::getDim(), GFTrack::getHit(), GFAbsRecoHit::getHitCov(), GFAbsRecoHit::getHMatrix(), GFAbsTrackRep::getReferencePlane(), GFTrack::getRepAtHit(), GFAbsTrackRep::getState(), GFTrack::getTrackRep(), GFAbsTrackRep::hasAuxInfo(), hit, r, res, GFAbsRecoHit::residualVector(), GFAbsTrackRep::setData(), GFBookkeeping::setDetPlane(), GFBookkeeping::setMatrix(), GFTrack::setRepAtHit(), and update().

Referenced by fittingPass().

                                                                 {
  GFAbsRecoHit* hit = tr->getHit(ihit);
  GFAbsTrackRep* rep = tr->getTrackRep(irep);

  // get prototypes for matrices
  int repDim = rep->getDim();
  TMatrixT<double> state(repDim,1);
  TMatrixT<double> cov(repDim,repDim);;
  GFDetPlane pl;
  /* do an extrapolation, if the trackrep irep is not given
   * at this ihit position. This will usually be the case, but
   * not if the fit turnes around
   */
  //std::cout << "fitting " << ihit << std::endl;
  if(ihit!=tr->getRepAtHit(irep)){
    //std::cout << "not same" << std::endl;
    // get the (virtual) detector plane
    pl=hit->getDetPlane(rep);
    //do the extrapolation
    rep->extrapolate(pl,state,cov);
  }
  else{
    //std::cout << "same" << std::endl;
    //    return;
    pl = rep->getReferencePlane();
    state = rep->getState();
    cov = rep->getCov();
  }
  
  if(cov[0][0]<1.E-50){
    GFException exc(COVEXC,__LINE__,__FILE__);
    throw exc;
  }
  
  //  TMatrixT<double> origcov=rep->getCov();
  TMatrixT<double> H=hit->getHMatrix(rep);
  // get hit covariances  
  TMatrixT<double> V=hit->getHitCov(pl);
  // calculate kalman gain ------------------------------
  TMatrixT<double> Gain(calcGain(cov,V,H));
  TMatrixT<double> res=hit->residualVector(rep,state,pl);
  // calculate update -----------------------------------
  TMatrixT<double> update=Gain*res;
  state+=update; // prediction overwritten!
  cov-=Gain*(H*cov);

  if(fSmooth) {
    if(direction == 1) {
          tr->getBK(irep)->setMatrix("fUpSt",ihit,state);
          tr->getBK(irep)->setMatrix("fUpCov",ihit,cov);
          if(rep->hasAuxInfo()) tr->getBK(irep)->setMatrix("fAuxInfo",ihit,*(rep->getAuxInfo(pl)));
          tr->getBK(irep)->setDetPlane("fPl",ihit,pl);
        } else {
          tr->getBK(irep)->setMatrix("bUpSt",ihit,state);
          tr->getBK(irep)->setMatrix("bUpCov",ihit,cov);
          if(rep->hasAuxInfo()) tr->getBK(irep)->setMatrix("bAuxInfo",ihit,*(rep->getAuxInfo(pl)));
          tr->getBK(irep)->setDetPlane("bPl",ihit,pl);
        }
  }

  // calculate filtered chisq
  // filtered residual
  TMatrixT<double> r=hit->residualVector(rep,state,pl);
  double chi2 = chi2Increment(r,H,cov,V);
  int ndf = r.GetNrows();
  rep->addChiSqu( chi2 );
  rep->addNDF( ndf );

  /*
  if(direction==1){
    tr->getBK(irep)->setNumber("fChi2",ihit,chi2/ndf);
  }
  else{
    tr->getBK(irep)->setNumber("bChi2",ihit,chi2/ndf);
  }
  */

  // if we survive until here: update TrackRep
  //rep->setState(state);
  //rep->setCov(cov);
  //rep->setReferencePlane(pl);
  rep->setData(state,pl,&cov);
  tr->setRepAtHit(irep,ihit);
}

void GFKalman::processTrack

(

GFTrack *

trk

)

Performs fit on a GFTrack.

The hits are processed in the order in which they are stored in the GFTrack object. Sorting of hits in space has to be done before!

Definition at line 40 of file GFKalman.cxx.

References blowUpCovs(), GFBookkeeping::bookGFDetPlanes(), GFBookkeeping::bookMatrices(), GFTrack::clearRepAtHit(), fInitialDirection, fittingPass(), fNumIt, fSmooth, GFTrack::getBK(), GFAbsTrackRep::getCov(), GFBookkeeping::getMatrixKeys(), GFTrack::getNumHits(), GFTrack::getNumReps(), GFAbsTrackRep::getReferencePlane(), GFTrack::getSmoothing(), GFAbsTrackRep::getState(), GFTrack::getTrackRep(), GFAbsTrackRep::hasAuxInfo(), i, GFTrack::setNextHitToFit(), GFBookkeeping::setNhits(), and switchDirection().

Referenced by PndHypKalmanTask::Exec(), PndHypDKalmanTask::Exec(), PndLmdKalmanTask::Exec(), PndRecoKalmanFit::Fit(), operator()(), and GFDaf::processTrack().

                                       {

  fSmooth = trk->getSmoothing();

  if(fSmooth) {
        int nreps = trk->getNumReps();
    for(int i=0; i<nreps; i++) {

                std::vector<std::string> mat_keys = trk->getBK(i)->getMatrixKeys();
                bool already_there = false;
                for(unsigned int j=0; j<mat_keys.size(); j++) {
                        if(mat_keys.at(j) == "fUpSt") already_there = true;
                }
                if(already_there) continue;

                trk->getBK(i)->setNhits(trk->getNumHits());
                trk->getBK(i)->bookMatrices("fUpSt");
                trk->getBK(i)->bookMatrices("fUpCov");
                trk->getBK(i)->bookMatrices("bUpSt");
                trk->getBK(i)->bookMatrices("bUpCov");
                trk->getBK(i)->bookGFDetPlanes("fPl");
                trk->getBK(i)->bookGFDetPlanes("bPl");
                if(trk->getTrackRep(i)->hasAuxInfo()) {
                        trk->getBK(i)->bookMatrices("fAuxInfo");
                        trk->getBK(i)->bookMatrices("bAuxInfo");
                }
    }
  }

  int direction=fInitialDirection;
  assert(direction==1 || direction==-1);
  //  trk->clearGFBookkeeping();
  trk->clearRepAtHit();

  /*why is there a factor of two here (in the for statement)?
    Because we consider one full iteration to be one back and
    one forth fitting pass */
  for(int ipass=0; ipass<2*fNumIt; ipass++){
    if(ipass>0) blowUpCovs(trk);
    if(direction==1){
      trk->setNextHitToFit(0);
    }
    else {
      trk->setNextHitToFit(trk->getNumHits()-1);
    }
    fittingPass(trk,direction);
    
    //save first and last plane,state&cov after the fitting pass
    if(direction==1){//forward at last hit
      int nreps=trk->getNumReps();
      for(int i=0; i<nreps; ++i){
        trk->getTrackRep(i)->
          setLastPlane( trk->getTrackRep(i)->getReferencePlane() );
        trk->getTrackRep(i)->
          setLastState( trk->getTrackRep(i)->getState() );
        trk->getTrackRep(i)->
          setLastCov( trk->getTrackRep(i)->getCov() );
      }
    }
    else{//backward at first hit
      int nreps=trk->getNumReps();
      for(int i=0; i<nreps; ++i){
        trk->getTrackRep(i)->
          setFirstPlane( trk->getTrackRep(i)->getReferencePlane() );
        trk->getTrackRep(i)->
          setFirstState( trk->getTrackRep(i)->getState() );
        trk->getTrackRep(i)->
          setFirstCov( trk->getTrackRep(i)->getCov() );
      }
    }

    //switch direction of fitting and also inside all the reps
    if(direction==1) direction=-1;
    else direction=1;
    switchDirection(trk);
  }
  
  return;
}

void GFKalman::setBlowUpFactor ( double  f )

inline

Set the blowup factor (see blowUpCovs() )

Definition at line 111 of file GFKalman.h.

References f, and fBlowUpFactor.

{fBlowUpFactor=f;}

void GFKalman::setInitialDirection ( int  d )

inline

Sets the inital direction of the track fit (1 for inner to outer, or -1 for outer to inner). The standard is 1 and is set in the ctor.

Definition at line 107 of file GFKalman.h.

References d, and fInitialDirection.

{fInitialDirection=d;}

void GFKalman::setLazy ( Int_t  )

inline

Switch lazy error handling.

This is a historically left-over method and shall be deleted some time

Definition at line 80 of file GFKalman.h.

{std::cerr<<"Using outdates setLazy method of class GFKalman:"<<std::endl;}

void GFKalman::setNumIterations ( Int_t  i )

inline

Set number of iterations for Kalman Filter.

One iteration is one forward pass plus one backward pass

Definition at line 86 of file GFKalman.h.

References fNumIt, and i.

Referenced by PndHypDKalmanTask::Exec(), PndLmdKalmanTask::Exec(), GFDaf::GFDaf(), and PndRecoKalmanFit::Init().

{fNumIt=i;}

void GFKalman::switchDirection ( GFTrack *  trk )

private

Used to switch between forward and backward filtering.

Definition at line 121 of file GFKalman.cxx.

References GFTrack::getNumReps(), GFTrack::getTrackRep(), i, and GFAbsTrackRep::switchDirection().

Referenced by processTrack().

                                     {
  int nreps=trk->getNumReps();
  for(int i=0; i<nreps; ++i){
    trk->getTrackRep(i)->switchDirection();
  }
}

Member Data Documentation

double GFKalman::fBlowUpFactor

private

Definition at line 147 of file GFKalman.h.

Referenced by blowUpCovs(), and setBlowUpFactor().

int GFKalman::fInitialDirection

private

Definition at line 145 of file GFKalman.h.

Referenced by processTrack(), and setInitialDirection().

Int_t GFKalman::fNumIt

private

Definition at line 146 of file GFKalman.h.

Referenced by processTrack(), and setNumIterations().

bool GFKalman::fSmooth

private

Definition at line 148 of file GFKalman.h.

Referenced by processHit(), and processTrack().

---

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

• GFKalman.h
• GFKalman.cxx