GFDaf Class Reference

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

#include <GFDaf.h>

Inheritance diagram for GFDaf:

Public Member Functions
	GFDaf ()
	~GFDaf ()
void	processTrack (GFTrack *trk)
	Process a track using the DAF. More...
const std::vector< std::vector < std::vector< double > > >	getWeights ()
	Return the weights present after the track was processed. More...
void	setProbCut (double prob_cut)
	Set the probabilty cut for the weight calculation for the hits. More...
void	setBetas (double b1, double b2, 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. More...

Private Member Functions
std::vector< GFDafHit * >	initHitsWeights (GFTrack *trk)
	Initialize the GFDafHit and their weights before the fit. More...
std::vector< std::vector < double > >	calcWeights (GFTrack *trk, double beta)
	Calculate the weights for the next fitting pass. More...
void	copySmoothing (GFTrack *source, GFTrack *target, int target_ire)
	Copy the smoothing matrices from the source track to the target. More...
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	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 Attributes
std::vector< std::vector < std::vector< double > > >	fWeights
std::vector< double >	fBeta
std::map< int, double >	fchi2Cuts

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.

Definition at line 46 of file GFDaf.h.

Constructor & Destructor Documentation

GFDaf::GFDaf

(

)

Definition at line 21 of file GFDaf.cxx.

References setBetas(), GFKalman::setNumIterations(), and setProbCut().

             : GFKalman::GFKalman() {

        setBetas(81.,8.,4.,1.,1.,1.);
        setProbCut(0.01);

        GFKalman::setNumIterations(1);

};

GFDaf::~GFDaf ( )

inline

Definition at line 50 of file GFDaf.h.

{ };

Member Function Documentation

std::vector< std::vector< double > > GFDaf::calcWeights ( GFTrack *  trk, double  beta  )

private

Calculate the weights for the next fitting pass.

Definition at line 80 of file GFDaf.cxx.

References exp(), fchi2Cuts, GFAbsRecoHit::getDetPlane(), GFDafHit::getHit(), GFTrack::getHit(), GFAbsRecoHit::getHitCoord(), GFAbsRecoHit::getHitCov(), GFDafHit::getNumHits(), GFTrack::getNumHits(), GFTools::getSmoothedPos(), GFTrack::getTrackRep(), i, GFException::info(), GFTools::invertMatrix(), m, phi, Pi, sqrt(), and GFException::what().

Referenced by processTrack().

                                                                          {

        std::vector<std::vector<double> > ret_val;

        for(unsigned int i=0; i<trk->getNumHits(); i++) {

                GFDafHit* eff_hit = static_cast<GFDafHit*>(trk->getHit(i));
                std::vector<double> phi;
                double phi_sum = 0;
                double phi_cut = 0;

                std::vector<double> weights;

                TMatrixT<double> x_smoo(GFTools::getSmoothedPos(trk, 0, i));
                if(x_smoo.GetNrows() == 0) {
                        weights.assign(eff_hit->getNumHits(),0.5);
                        std::cout<<"Assumed weight 0.5!!"<<std::endl;
                        ret_val.push_back(weights);
                        continue;
                }

                for(unsigned int j=0; j<eff_hit->getNumHits(); j++) {
                        GFAbsRecoHit* real_hit = eff_hit->getHit(j);
                        GFDetPlane pl;
                        try{
                                pl = real_hit->getDetPlane(trk->getTrackRep(0));
                        } catch(GFException& e) {
                                std::cerr<<e.what();
                                e.info();
                        }
                        TMatrixT<double> m(real_hit->getHitCoord(pl));
                        TMatrixT<double> V( beta * real_hit->getHitCov(pl));
                        TMatrixT<double> resid(m - x_smoo);
                        TMatrixT<double> resid_T(resid);
                        resid_T.T();
                        double detV = V.Determinant();
                        TMatrixT<double> Vinv;
                        GFTools::invertMatrix(V,Vinv);

                        phi.push_back((1./(pow(2.*TMath::Pi(),V.GetNrows()/2.)*sqrt(detV)))*exp(-0.5*(resid_T * Vinv * resid)[0][0]));
                        phi_sum += phi.at(j);

                        double cutVal = fchi2Cuts[V.GetNrows()];
                        assert(cutVal>1.E-6);
                        phi_cut += 1./(2.*TMath::Pi()*sqrt(detV))*exp(-0.5*cutVal/beta);

                }

                for(unsigned int j=0; j<eff_hit->getNumHits(); j++) {

//std::cout<<"Calculated weight (beta = "<<beta<<", hit "<<i<<") "<<phi.at(j)<<"/("<<phi_sum<<"+"<<phi_cut<<") = "<<phi.at(j)/(phi_sum+phi_cut)<<std::endl;
                        weights.push_back(phi.at(j)/(phi_sum+phi_cut));

                }

                ret_val.push_back(weights);

        }

        return ret_val;

};

void GFDaf::copySmoothing ( GFTrack *  source, GFTrack *  target, int  target_ire  )

private

Copy the smoothing matrices from the source track to the target.

Definition at line 238 of file GFDaf.cxx.

References GFBookkeeping::bookGFDetPlanes(), GFBookkeeping::bookMatrices(), GFTrack::getBK(), GFBookkeeping::getDetPlane(), GFTrack::getHit(), GFBookkeeping::getMatrix(), GFBookkeeping::getMatrixKeys(), GFDafHit::getNumHits(), GFTrack::getNumHits(), GFTrack::getNumReps(), GFTrack::getTrackRep(), GFAbsTrackRep::hasAuxInfo(), i, GFBookkeeping::setDetPlane(), GFBookkeeping::setMatrix(), and GFBookkeeping::setNhits().

Referenced by processTrack().

                                                                           {

        for(unsigned int i=0; i<target->getNumReps(); i++) {

                std::vector<std::string> mat_keys = target->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;

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

        int hit_count = 0;

        for(unsigned int pl_i=0; pl_i<source->getNumHits(); pl_i++) {

                GFDafHit* eff_hit = static_cast<GFDafHit*>(source->getHit(pl_i));

                for(unsigned int hit_i=0; hit_i<eff_hit->getNumHits(); hit_i++) {

                        TMatrixT<double> fUpSt, fUpCov, bUpSt, bUpCov, fAuxInfo, bAuxInfo;
                        GFDetPlane fPl, bPl;

                        source->getBK(0)->getMatrix("fUpSt",hit_i,fUpSt);
                        source->getBK(0)->getMatrix("fUpCov",hit_i,fUpCov);
                        source->getBK(0)->getDetPlane("fPl",hit_i,fPl);
                        source->getBK(0)->getMatrix("bUpSt",hit_i,bUpSt);
                        source->getBK(0)->getMatrix("bUpCov",hit_i,bUpCov);
                        source->getBK(0)->getDetPlane("bPl",hit_i,bPl);

                        if(source->getTrackRep(0)->hasAuxInfo()) {
                                source->getBK(0)->getMatrix("fAuxInfo",hit_i,fAuxInfo);
                                source->getBK(0)->getMatrix("bAuxInfo",hit_i,bAuxInfo);
                                target->getBK(target_irep)->setMatrix("fAuxInfo",hit_count,fAuxInfo);
                                target->getBK(target_irep)->setMatrix("bAuxInfo",hit_count,bAuxInfo);
                        }

                        target->getBK(target_irep)->setMatrix("fUpSt",hit_count,fUpSt);
                        target->getBK(target_irep)->setMatrix("fUpCov",hit_count,fUpCov);
                        target->getBK(target_irep)->setDetPlane("fPl",hit_count,fPl);
                        target->getBK(target_irep)->setMatrix("bUpSt",hit_count,bUpSt);
                        target->getBK(target_irep)->setMatrix("bUpCov",hit_count,bUpCov);
                        target->getBK(target_irep)->setDetPlane("bPl",hit_count,bPl);

                        hit_count++;

                }

        }

}

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

inherited

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, GFKalman::processHit(), GFAbsTrackRep::setChiSqu(), GFAbsTrackRep::setNDF(), GFTrack::setNextHitToFit(), GFAbsTrackRep::setStatusFlag(), and GFException::what().

Referenced by GFKalman::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  )

inherited

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

Definition at line 215 of file GFKalman.cxx.

References GFKalman::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();
}

const std::vector<std::vector<std::vector<double> > > GFDaf::getWeights ( )

inline

Return the weights present after the track was processed.

The DAF uses special effective hits defined in the class GFDafHit. A GFDafHit is a wrappe class and contains all the real hits from one plane. The structure of the return vector of getWeights allows to reconstruct in what way the hits were grouped: the outermost vector represents the track representation, there is one entry per track representation. The middle vector represents the effective hits, and the innermost vector contains the real hits contained in the corresponding effective hit.

Definition at line 66 of file GFDaf.h.

References fWeights.

{ return fWeights; };

std::vector< GFDafHit * > GFDaf::initHitsWeights ( GFTrack *  trk )

private

Initialize the GFDafHit and their weights before the fit.

Definition at line 203 of file GFDaf.cxx.

References fWeights, GFTrack::getHit(), GFTrack::getHitsByPlane(), GFTrack::getNumReps(), hits, and i.

Referenced by processTrack().

                                                        {

        std::vector< std::vector<int>* > planes;
        trk->getHitsByPlane(planes);
        int nPlanes = planes.size();

        std::vector<GFDafHit*> eff_hits;

        for(int i=0; i<nPlanes; i++) {

                std::vector<GFAbsRecoHit*> hits;

                for(unsigned int j=0; j<planes.at(i)->size(); j++) {
                        hits.push_back(trk->getHit(planes.at(i)->at(j)) );
                }

                GFDafHit* eff_hit = new GFDafHit(hits);
                eff_hits.push_back(eff_hit);

        }

        for(unsigned int i=0; i<trk->getNumReps(); i++) {
                std::vector<std::vector<double> > rep_weights;
                for(unsigned int j=0; j<eff_hits.size(); j++) {
                        std::vector<double> single_weights;
                        single_weights.assign(eff_hits.at(j)->getNumHits(),1.);
                        rep_weights.push_back(single_weights);
                }
                fWeights.push_back(rep_weights);
        }

        return eff_hits;

}

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

inlineinherited

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 GFKalman::processTrack().

{processTrack(track);}

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

inlineinherited

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 GFKalman::processTrack().

{processTrack(tr.second);}

void GFDaf::processTrack

(

GFTrack *

trk

)

Process a track using the DAF.

Definition at line 30 of file GFDaf.cxx.

References GFTrack::addHit(), GFTrack::addTrackRep(), calcWeights(), copySmoothing(), fBeta, fWeights, GFTrack::getHit(), GFTrack::getNumHits(), GFTrack::getNumReps(), GFTrack::getSmoothing(), GFAbsTrackRep::getStatusFlag(), GFTrack::getTrackRep(), hit, i, GFException::info(), initHitsWeights(), GFKalman::processTrack(), GFTrack::releaseTrackReps(), GFDafHit::setBlowUp(), GFTrack::setSmoothing(), GFAbsTrackRep::setStatusFlag(), GFDafHit::setWeights(), and GFException::what().

Referenced by PndRecoDafFit::Fit().

                                     {

        std::vector<GFDafHit*> eff_hits = initHitsWeights(trk);

        GFTrack* mini_trk = new GFTrack();

        for(unsigned int j=0; j<eff_hits.size(); j++) mini_trk->addHit(eff_hits.at(j));

        mini_trk->setSmoothing();

        for(unsigned int i=0; i<trk->getNumReps(); i++) {

                if(trk->getTrackRep(i)->getStatusFlag()!=0) continue;

                mini_trk->addTrackRep(trk->getTrackRep(i));

                for(unsigned int iBeta=0; iBeta<fBeta.size(); iBeta++) {

                        for(unsigned int j=0; j<mini_trk->getNumHits(); j++) {
                                GFDafHit* hit = static_cast<GFDafHit*>(mini_trk->getHit(j));
                                hit->setWeights(fWeights.at(i).at(j));
                                hit->setBlowUp(fBeta.at(iBeta));
                        }

                        GFKalman::processTrack(mini_trk);

                        if(mini_trk->getTrackRep(0)->getStatusFlag() != 0) break;

                        if(iBeta != fBeta.size()-1 ) 
                                try{
                                        fWeights.at(i) = calcWeights(mini_trk, fBeta.at(iBeta));
                                } catch(GFException& e) {
                                        std::cerr<<e.what();
                                        e.info();
                                        mini_trk->getTrackRep(0)->setStatusFlag(1);
                                        break;
                                }

                }

                if(trk->getSmoothing()) copySmoothing(mini_trk, trk, i);

                mini_trk->releaseTrackReps();

        }

        delete mini_trk;

};

void GFDaf::setBetas	(	double	b1,
		double	b2,
		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 two temperatures. The maximum would ten tempertatures.

Definition at line 174 of file GFDaf.cxx.

References fBeta.

Referenced by GFDaf().

                                                                                                                        {
        assert(b1>0);fBeta.push_back(b1);
        assert(b2>0 && b2<=b1);fBeta.push_back(b2);
        if(b3>=0.) {
                assert(b3<=b2);fBeta.push_back(b3);
                if(b4>=0.) {
                        assert(b4<=b3);fBeta.push_back(b4);
                        if(b5>=0.) {
                                assert(b5<=b4);fBeta.push_back(b5);
                                if(b6>=0.) {
                                        assert(b6<=b5);fBeta.push_back(b6);
                                        if(b7>=0.) {
                                                assert(b7<=b6);fBeta.push_back(b7);
                                                if(b8>=0.) {
                                                        assert(b8<=b7);fBeta.push_back(b8);
                                                        if(b9>=0.) {
                                                                assert(b9<=b8);fBeta.push_back(b9);
                                                                if(b10>=0.) {
                                                                        assert(b10<=b9);fBeta.push_back(b10);
                                                                }
                                                        }
                                                }
                                        }
                                }
                        }
                }
        }
}

void GFKalman::setBlowUpFactor ( double  f )

inlineinherited

Set the blowup factor (see blowUpCovs() )

Definition at line 111 of file GFKalman.h.

References f, and GFKalman::fBlowUpFactor.

{fBlowUpFactor=f;}

void GFKalman::setInitialDirection ( int  d )

inlineinherited

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 GFKalman::fInitialDirection.

{fInitialDirection=d;}

void GFKalman::setLazy ( Int_t  )

inlineinherited

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 )

inlineinherited

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 GFKalman::fNumIt, and i.

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

{fNumIt=i;}

void GFDaf::setProbCut

(

double

prob_cut

)

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

Currently supported are the values 0.01 0.005, and 0.001. The corresponding chi2 cuts for different hits dimensionalities are hardcoded in the implementation because I did not yet figure out how to calculate them. Please feel very welcome to change the implementtion if you know how to do it.

Definition at line 143 of file GFDaf.cxx.

References fabs(), fchi2Cuts, and GFException::setFatal().

Referenced by GFDaf().

                                {

        if(fabs(val-0.01)<1.E-10){
                fchi2Cuts[1] = 6.63;
                fchi2Cuts[2] = 9.21;
                fchi2Cuts[3] = 11.34;
                fchi2Cuts[4] = 13.23;
                fchi2Cuts[5] = 15.09;
        }
        else   if(fabs(val-0.005)<1.E-10){
                fchi2Cuts[1] = 7.88;
                fchi2Cuts[2] = 10.60;
                fchi2Cuts[3] = 12.84;
                fchi2Cuts[4] = 14.86;
                fchi2Cuts[5] = 16.75;
        }
        else   if(fabs(val-0.001)<1.E-10){
                fchi2Cuts[1] = 10.83;
                fchi2Cuts[2] = 13.82;
                fchi2Cuts[3] = 16.27;
                fchi2Cuts[4] = 18.47;
                fchi2Cuts[5] = 20.51;
        }
        else{
                GFException exc("GFDafsetProbCut() value is not supported",__LINE__,__FILE__);
                exc.setFatal();
                throw exc;
        }

}

Member Data Documentation

std::vector<double> GFDaf::fBeta

private

Definition at line 100 of file GFDaf.h.

Referenced by processTrack(), and setBetas().

std::map<int,double> GFDaf::fchi2Cuts

private

Definition at line 101 of file GFDaf.h.

Referenced by calcWeights(), and setProbCut().

std::vector<std::vector<std::vector<double> > > GFDaf::fWeights

private

Definition at line 99 of file GFDaf.h.

Referenced by getWeights(), initHitsWeights(), and processTrack().

---

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

• GFDaf.h
• GFDaf.cxx