GFTools Namespace Reference

Namespace for various tools, mainly smoothing. More...

Functions
TMatrixT< double >	getSmoothedPos (GFTrack *trk, unsigned int irep, unsigned int ihit)
	Get smoothed track position in plane coordinates. More...
TMatrixT< double >	getSmoothedCov (GFTrack *trk, unsigned int irep, unsigned int ihit)
	Get smoothed track covariance in plane coordinates. More...
bool	getSmoothedData (GFTrack *trk, unsigned int irep, unsigned int ihit, TMatrixT< double > &smoothed_state, TMatrixT< double > &smoothed_cov)
	Get smoothed state vector and state covariance. More...
bool	getSmoothedData (GFTrack *trk, unsigned int irep, unsigned int ihit, TMatrixT< double > &smoothed_state, TMatrixT< double > &smoothed_cov, GFDetPlane &smoothing_plane)
	Get smoothed state vector, state covariance and smoothing plane. More...
bool	getSmoothedData (GFTrack *trk, unsigned int irep, unsigned int ihit, TMatrixT< double > &smoothed_state, TMatrixT< double > &smoothed_cov, GFDetPlane &smoothing_plane, TMatrixT< double > &auxInfo)
	Get smoothed state vector, state covariance and smoothing plane. More...
GFDetPlane	getSmoothingPlane (GFTrack *trk, unsigned int irep, unsigned int ihit)
	Get smoothing plane. More...
bool	getBiasedSmoothedData (GFTrack *trk, unsigned int irep, unsigned int ihit, TMatrixT< double > &smoothed_state, TMatrixT< double > &smoothed_cov)
	Get biased smoothed state vector and state covariance. More...
bool	getBiasedSmoothedData (GFTrack *trk, unsigned int irep, unsigned int ihit, TMatrixT< double > &smoothed_state, TMatrixT< double > &smoothed_cov, GFDetPlane &smoothing_plane)
	Get biased smoothed state vector, state covariance and smoothing plane. More...
bool	getBiasedSmoothedData (GFTrack *trk, unsigned int irep, unsigned int ihit, TMatrixT< double > &smoothed_state, TMatrixT< double > &smoothed_cov, GFDetPlane &smoothing_plane, TMatrixT< double > &auxInfo)
	Get biased smoothed state vector, state covariance and smoothing plane. More...
TMatrixT< double >	getBiasedSmoothedPos (GFTrack *trk, unsigned int irep, unsigned int ihit)
	Get biased smoothed track position in plane coordinates. More...
TMatrixT< double >	getBiasedSmoothedCov (GFTrack *trk, unsigned int irep, unsigned int ihit)
	Get biased smoothed track covariance in plane coordinates. More...
void	invertMatrix (const TMatrixT< double > &mat, TMatrixT< double > &inv)
	Invert a matrix, throwing GFException when inversion fails. More...
double	getSmoothedChiSqu (GFTrack *const trk, unsigned int irep, unsigned int ihit)
	Get smoothed chi2 for a specific hit (ihit). More...

Detailed Description

Namespace for various tools, mainly smoothing.

The GFTools namespace contains various functions, most of which er intended to be used with the smoothing feature of the Kalman filter. These functions allow to extract the smoothed information for any track which was fitted with smoothing enabled.

Function Documentation

TMatrixT< double > GFTools::getBiasedSmoothedCov	(	GFTrack *	trk,
		unsigned int	irep,
		unsigned int	ihit
	)

Get biased smoothed track covariance in plane coordinates.

Definition at line 57 of file GFTools.cxx.

References getBiasedSmoothedData(), GFTrack::getHit(), GFAbsRecoHit::getHMatrix(), and GFTrack::getTrackRep().

                                                                                                 {

        TMatrixT<double> smoothed_state;
        TMatrixT<double> smoothed_cov;
        GFDetPlane pl;

        GFTools::getBiasedSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov, pl);
        TMatrixT<double> H = trk->getHit(ihit)->getHMatrix(trk->getTrackRep(irep));

        TMatrixT<double> cov_tmp(smoothed_cov,TMatrixT<double>::kMultTranspose,H);
        TMatrixT<double> cov(H,TMatrixT<double>::kMult,cov_tmp);

        return cov;

}

bool GFTools::getBiasedSmoothedData	(	GFTrack *	trk,
		unsigned int	irep,
		unsigned int	ihit,
		TMatrixT< double > &	smoothed_state,
		TMatrixT< double > &	smoothed_cov
	)

Get biased smoothed state vector and state covariance.

Definition at line 81 of file GFTools.cxx.

Referenced by getBiasedSmoothedCov(), getBiasedSmoothedData(), getBiasedSmoothedPos(), and getSmoothedChiSqu().

                                                                                                                                                        {

        GFDetPlane pl;
        TMatrixT<double> auxInfo;
        return GFTools::getBiasedSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov, pl, auxInfo);

}

bool GFTools::getBiasedSmoothedData	(	GFTrack *	trk,
		unsigned int	irep,
		unsigned int	ihit,
		TMatrixT< double > &	smoothed_state,
		TMatrixT< double > &	smoothed_cov,
		GFDetPlane &	smoothing_plane
	)

Get biased smoothed state vector, state covariance and smoothing plane.

The smoothed data returned by this function includes the detector plane in which the data is calculated.

Definition at line 96 of file GFTools.cxx.

References getBiasedSmoothedData().

                                                                                                                                                                                     {

        TMatrixT<double> auxInfo;
        return GFTools::getBiasedSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov, smoothing_plane, auxInfo);

}

bool GFTools::getBiasedSmoothedData	(	GFTrack *	trk,
		unsigned int	irep,
		unsigned int	ihit,
		TMatrixT< double > &	smoothed_state,
		TMatrixT< double > &	smoothed_cov,
		GFDetPlane &	smoothing_plane,
		TMatrixT< double > &	auxInfo
	)

Get biased smoothed state vector, state covariance and smoothing plane.

The smoothed data returned by this function includes the detector plane in which the data is calculated as well as the auxillary information for this plane.

Definition at line 216 of file GFTools.cxx.

References GFAbsTrackRep::clone(), GFAbsTrackRep::extrapolate(), GFTrack::getBK(), GFBookkeeping::getDetPlane(), GFBookkeeping::getMatrix(), GFTrack::getNumHits(), GFTrack::getSmoothing(), GFTrack::getTrackRep(), GFAbsTrackRep::hasAuxInfo(), invertMatrix(), and GFAbsTrackRep::setData().

                                                                                                                                                                                                                {

        if(!trk->getSmoothing()) {
                std::cout << "Trying to get smoothed hit coordinates from a track without smoothing! Aborting..." << std::endl;
                TMatrixT<double> failed(1,1);
                return false;
        }

        if(ihit >= trk->getNumHits()) {
                std::cout << "Hit number out of bounds while trying to get smoothed hit coordinates! Aborting..." <<std::endl;
                TMatrixT<double> failed(1,1);
                failed(0,0) = 0;
                return false;
        }

        TMatrixT<double> bUpSt;
        TMatrixT<double> bUpCov;
        TMatrixT<double> bAuxInfo;
        TMatrixT<double>* bAuxInfoP;
        GFDetPlane bPl;

        GFAbsTrackRep* rep = trk->getTrackRep(irep)->clone();
        
        if(trk->getTrackRep(irep)->hasAuxInfo()) {
                trk->getBK(irep)->getMatrix("fAuxInfo",ihit,auxInfo);
        }

        if(ihit == 0) {
                trk->getBK(irep)->getMatrix("bUpSt",ihit,smoothed_state);
                trk->getBK(irep)->getMatrix("bUpCov",ihit,smoothed_cov);
                trk->getBK(irep)->getDetPlane("bPl",ihit,smoothing_plane);
                return true;
        }

        if(ihit == trk->getNumHits()-1) {
                trk->getBK(irep)->getMatrix("fUpSt",ihit,smoothed_state);
                trk->getBK(irep)->getMatrix("fUpCov",ihit,smoothed_cov);
                trk->getBK(irep)->getDetPlane("fPl",ihit,smoothing_plane);
                return true;
        }

        TMatrixT<double> fSt;
        TMatrixT<double> fCov;

        trk->getBK(irep)->getMatrix("fUpSt",ihit,fSt);
        trk->getBK(irep)->getMatrix("fUpCov",ihit,fCov);
        trk->getBK(irep)->getMatrix("bUpSt",ihit+1,bUpSt);
        trk->getBK(irep)->getMatrix("bUpCov",ihit+1,bUpCov);
        if(trk->getTrackRep(irep)->hasAuxInfo()) {
                trk->getBK(irep)->getMatrix("bAuxInfo",ihit+1,bAuxInfo);
                bAuxInfoP = &bAuxInfo;
        } else {
                bAuxInfoP = NULL;
        }
        trk->getBK(irep)->getDetPlane("fPl",ihit,smoothing_plane);
        trk->getBK(irep)->getDetPlane("bPl",ihit+1,bPl);


        TMatrixT<double> bSt;
        TMatrixT<double> bCov;

        if(bUpSt.GetNrows() == 0) return false;

        rep->setData(bUpSt,bPl,&bUpCov,bAuxInfoP);
        rep->extrapolate(smoothing_plane,bSt,bCov);

        delete rep;

        TMatrixT<double> fCovInvert;
        TMatrixT<double> bCovInvert;

        GFTools::invertMatrix(fCov, fCovInvert);
        GFTools::invertMatrix(bCov, bCovInvert);

        GFTools::invertMatrix(fCovInvert + bCovInvert, smoothed_cov);

        smoothed_state.ResizeTo(smoothed_cov * (fCovInvert*fSt + bCovInvert*bSt));
        smoothed_state = smoothed_cov * (fCovInvert*fSt + bCovInvert*bSt);

        return true;

}

TMatrixT< double > GFTools::getBiasedSmoothedPos	(	GFTrack *	trk,
		unsigned int	irep,
		unsigned int	ihit
	)

Get biased smoothed track position in plane coordinates.

Definition at line 22 of file GFTools.cxx.

References getBiasedSmoothedData(), GFTrack::getHit(), GFAbsRecoHit::getHMatrix(), GFTrack::getTrackRep(), and pos.

                                                                                                 {

TMatrixT<double> smoothed_state;
        TMatrixT<double> smoothed_cov;
        TMatrixT<double> pos;

        if(GFTools::getBiasedSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov)) {

                TMatrixT<double> H = trk->getHit(ihit)->getHMatrix(trk->getTrackRep(irep));
                H.Print();smoothed_state.Print();
                TMatrixT<double> pos_tmp(H * smoothed_state);
                pos.ResizeTo(pos_tmp);
                pos = pos_tmp;

        }
        return pos;

}

double GFTools::getSmoothedChiSqu	(	GFTrack *const	trk,
		unsigned int	irep,
		unsigned int	ihit
	)

Get smoothed chi2 for a specific hit (ihit).

This function calculates a smoothed chi2 value for a specific hit from the the so called "biased" smoothed state vector and covariance matrix. If many values from different tracks in the same layer are calculated they should be chi2 distributed with dim(m) degree of freedom, where m is the measurement vector of the hit. So for a pixel detector that measures x and y dim(m) will be 2.

Definition at line 358 of file GFTools.cxx.

References getBiasedSmoothedData(), GFTrack::getHit(), GFAbsRecoHit::getHMatrix(), GFAbsRecoHit::getRawHitCoord(), GFAbsRecoHit::getRawHitCov(), GFTrack::getTrackRep(), invertMatrix(), m, pos, R, and res.

                                                                                         {
        TMatrixT<double> smoothed_state;
        TMatrixT<double> smoothed_cov;
        GFTools::getBiasedSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov);
        TMatrixT<double> H = trk->getHit(ihit)->getHMatrix(trk->getTrackRep(irep));
        TMatrixT<double> HT(TMatrixT<double>::kTransposed, H);
        TMatrixT<double> pos = H * smoothed_state;
        TMatrixT<double> m = trk->getHit(ihit)->getRawHitCoord(); //measurement of hit
        TMatrixT<double> V = trk->getHit(ihit)->getRawHitCov(); //covariance matrix of hit
        TMatrixT<double> res = m - pos;
        TMatrixT<double> resT(TMatrixT<double>::kTransposed, res);
        TMatrixT<double> R = V - H*smoothed_cov*HT;
        TMatrixT<double> invR;
        invertMatrix(R,invR);
        TMatrixT<double> smoothedChiSqu = resT*invR*res;
        return smoothedChiSqu[0][0];
}

TMatrixT< double > GFTools::getSmoothedCov	(	GFTrack *	trk,
		unsigned int	irep,
		unsigned int	ihit
	)

Get smoothed track covariance in plane coordinates.

Definition at line 41 of file GFTools.cxx.

References GFTrack::getHit(), GFAbsRecoHit::getHMatrix(), getSmoothedData(), and GFTrack::getTrackRep().

                                                                                           {

        TMatrixT<double> smoothed_state;
        TMatrixT<double> smoothed_cov;
        GFDetPlane pl;

        GFTools::getSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov, pl);
        TMatrixT<double> H = trk->getHit(ihit)->getHMatrix(trk->getTrackRep(irep));

        TMatrixT<double> cov_tmp(smoothed_cov,TMatrixT<double>::kMultTranspose,H);
        TMatrixT<double> cov(H,TMatrixT<double>::kMult,cov_tmp);

        return cov;

}

bool GFTools::getSmoothedData	(	GFTrack *	trk,
		unsigned int	irep,
		unsigned int	ihit,
		TMatrixT< double > &	smoothed_state,
		TMatrixT< double > &	smoothed_cov
	)

Get smoothed state vector and state covariance.

Definition at line 73 of file GFTools.cxx.

Referenced by getSmoothedCov(), getSmoothedData(), and getSmoothedPos().

                                                                                                                                                  {

        GFDetPlane pl;
        TMatrixT<double> auxInfo;
        return GFTools::getSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov, pl, auxInfo);

}

bool GFTools::getSmoothedData	(	GFTrack *	trk,
		unsigned int	irep,
		unsigned int	ihit,
		TMatrixT< double > &	smoothed_state,
		TMatrixT< double > &	smoothed_cov,
		GFDetPlane &	smoothing_plane
	)

Get smoothed state vector, state covariance and smoothing plane.

The smoothed data returned by this function includes the detector plane in which the data is calculated.

Definition at line 89 of file GFTools.cxx.

References getSmoothedData().

                                                                                                                                                                               {

        TMatrixT<double> auxInfo;
        return GFTools::getSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov, smoothing_plane, auxInfo);

}

bool GFTools::getSmoothedData	(	GFTrack *	trk,
		unsigned int	irep,
		unsigned int	ihit,
		TMatrixT< double > &	smoothed_state,
		TMatrixT< double > &	smoothed_cov,
		GFDetPlane &	smoothing_plane,
		TMatrixT< double > &	auxInfo
	)

Get smoothed state vector, state covariance and smoothing plane.

The smoothed data returned by this function includes the detector plane in which the data is calculated as well as the auxillary information for this plane. The smoothed covariance matrix at hit i is:

\[ C_{k,\mbox{smooth}} = (C_{k,k-1,\rightarrow}^{-1}+C_{k,k+1,\leftarrow}^{-1})^{-1}. \\ \]

The smoothed state vector at hit i is:

\[ x_{k,\mbox{smooth}} = C_{k,\mbox{smooth}} \cdot ((C_{k,k-1,\rightarrow}^{-1}\cdot x_{k,k-1,\rightarrow} + C_{k,k+1,\leftarrow}^{-1}\cdot x_{k,k+1,\leftarrow}) \]

The index k,k-1 denotes that the state vector or covariance matrix contains the information of all hits 1,2,....,k-1 and is extrapolated to hit k. The arrow pointing to the right indicates that the information is saved during the forward fit. The index k,k+1 indicates the the state vector or covariance matrix contains the information of all hits N,N-1,...,k+1 and is extrapolated to hit k. The left-poiting arrow indicates that the information is saved during the backward fit.

Definition at line 103 of file GFTools.cxx.

References GFAbsTrackRep::clone(), GFAbsTrackRep::extrapolate(), GFTrack::getBK(), GFBookkeeping::getDetPlane(), GFBookkeeping::getMatrix(), GFTrack::getNumHits(), GFTrack::getSmoothing(), GFTrack::getTrackRep(), GFAbsTrackRep::hasAuxInfo(), invertMatrix(), and GFAbsTrackRep::setData().

                                                                                                                                                                                                          {

        if(!trk->getSmoothing()) {
                std::cout << "Trying to get smoothed hit coordinates from a track without smoothing! Aborting..." << std::endl;
                TMatrixT<double> failed(1,1);
                return false;
        }

        if(ihit >= trk->getNumHits()) {
                std::cout << "Hit number out of bounds while trying to get smoothed hit coordinates! Aborting..." <<std::endl;
                TMatrixT<double> failed(1,1);
                failed(0,0) = 0;
                return false;
        }

        TMatrixT<double> fUpSt;
        TMatrixT<double> fUpCov;
        TMatrixT<double> fAuxInfo;
        TMatrixT<double>* fAuxInfoP;
        TMatrixT<double> bUpSt;
        TMatrixT<double> bUpCov;
        TMatrixT<double> bAuxInfo;
        TMatrixT<double>* bAuxInfoP;

        GFDetPlane fPl;
        GFDetPlane bPl;

        GFAbsTrackRep* rep = trk->getTrackRep(irep)->clone();

        if(trk->getTrackRep(irep)->hasAuxInfo()) {
                trk->getBK(irep)->getMatrix("fAuxInfo",ihit,auxInfo);
        }

        if(ihit == 0) {
                trk->getBK(irep)->getMatrix("bUpSt",ihit+1,bUpSt);
                trk->getBK(irep)->getMatrix("bUpCov",ihit+1,bUpCov);
                trk->getBK(irep)->getDetPlane("fPl",ihit,smoothing_plane);
                trk->getBK(irep)->getDetPlane("bPl",ihit+1,bPl);
                if(trk->getTrackRep(irep)->hasAuxInfo()) {
                        trk->getBK(irep)->getMatrix("bAuxInfo",ihit+1,bAuxInfo);
                        bAuxInfoP = &bAuxInfo;
                } else {
                        bAuxInfoP = NULL;
                }
                if(bUpSt.GetNrows() == 0) return false;
                rep->setData(bUpSt,bPl,&bUpCov,bAuxInfoP);
                rep->extrapolate(smoothing_plane,smoothed_state,smoothed_cov);
                return true;
        }

        if(ihit == trk->getNumHits()-1) {
                trk->getBK(irep)->getMatrix("fUpSt",ihit-1,fUpSt);
                trk->getBK(irep)->getMatrix("fUpCov",ihit-1,fUpCov);
                trk->getBK(irep)->getDetPlane("fPl",ihit-1,fPl);
                trk->getBK(irep)->getDetPlane("fPl",ihit,smoothing_plane);
                if(trk->getTrackRep(irep)->hasAuxInfo()) {
                        trk->getBK(irep)->getMatrix("fAuxInfo",ihit-1,fAuxInfo);
                        fAuxInfoP = &fAuxInfo;
                } else {
                        fAuxInfoP = NULL;
                }
                if(fUpSt.GetNrows() == 0) return false;
                rep->setData(fUpSt,fPl,&fUpCov,fAuxInfoP);
                rep->extrapolate(smoothing_plane,smoothed_state,smoothed_cov);
                return true;
        }

        trk->getBK(irep)->getMatrix("fUpSt",ihit-1,fUpSt);
        trk->getBK(irep)->getMatrix("fUpCov",ihit-1,fUpCov);
        trk->getBK(irep)->getMatrix("bUpSt",ihit+1,bUpSt);
        trk->getBK(irep)->getMatrix("bUpCov",ihit+1,bUpCov);
        if(trk->getTrackRep(irep)->hasAuxInfo()) {
                trk->getBK(irep)->getMatrix("fAuxInfo",ihit-1,fAuxInfo);
                trk->getBK(irep)->getMatrix("bAuxInfo",ihit+1,bAuxInfo);
                fAuxInfoP = &fAuxInfo;
                bAuxInfoP = &bAuxInfo;
        } else {
                fAuxInfoP = bAuxInfoP = NULL;
        }
        trk->getBK(irep)->getDetPlane("fPl",ihit-1,fPl);
        trk->getBK(irep)->getDetPlane("fPl",ihit,smoothing_plane);
        trk->getBK(irep)->getDetPlane("bPl",ihit+1,bPl);


        TMatrixT<double> fSt;
        TMatrixT<double> fCov;
        TMatrixT<double> bSt;
        TMatrixT<double> bCov;

        if(fUpSt.GetNrows() == 0 || bUpSt.GetNrows() == 0) return false;

        rep->setData(fUpSt,fPl,&fUpCov,fAuxInfoP);
        rep->extrapolate(smoothing_plane,fSt,fCov);
        rep->setData(bUpSt,bPl,&bUpCov,bAuxInfoP);
        rep->extrapolate(smoothing_plane,bSt,bCov);

        delete rep;

        TMatrixT<double> fCovInvert;
        TMatrixT<double> bCovInvert;

        GFTools::invertMatrix(fCov, fCovInvert);
        GFTools::invertMatrix(bCov, bCovInvert);

        GFTools::invertMatrix(fCovInvert + bCovInvert, smoothed_cov);

        smoothed_state.ResizeTo(smoothed_cov * (fCovInvert*fSt + bCovInvert*bSt));
        smoothed_state = smoothed_cov * (fCovInvert*fSt + bCovInvert*bSt);

        return true;

}

TMatrixT< double > GFTools::getSmoothedPos	(	GFTrack *	trk,
		unsigned int	irep,
		unsigned int	ihit
	)

Get smoothed track position in plane coordinates.

Definition at line 4 of file GFTools.cxx.

References GFTrack::getHit(), GFAbsRecoHit::getHMatrix(), getSmoothedData(), GFTrack::getTrackRep(), and pos.

Referenced by GFDaf::calcWeights().

                                                                                           {

        TMatrixT<double> smoothed_state;
        TMatrixT<double> smoothed_cov;
        TMatrixT<double> pos;

        if(GFTools::getSmoothedData(trk, irep, ihit, smoothed_state, smoothed_cov)) {

                TMatrixT<double> H = trk->getHit(ihit)->getHMatrix(trk->getTrackRep(irep));
                TMatrixT<double> pos_tmp(H * smoothed_state);
                pos.ResizeTo(pos_tmp);
                pos = pos_tmp;

        }
        return pos;

}

GFDetPlane GFTools::getSmoothingPlane	(	GFTrack *	trk,
		unsigned int	irep,
		unsigned int	ihit
	)

Get smoothing plane.

Definition at line 299 of file GFTools.cxx.

References GFTrack::getBK(), and GFBookkeeping::getDetPlane().

                                                                                        {

        GFDetPlane pl;
        trk->getBK(irep)->getDetPlane("fPl",ihit,pl);
        return pl;

}

void GFTools::invertMatrix	(	const TMatrixT< double > &	mat,
		TMatrixT< double > &	inv
	)

Invert a matrix, throwing GFException when inversion fails.

Definition at line 334 of file GFTools.cxx.

References GFException::setFatal(), and status.

Referenced by GFKalman::calcGain(), GFDaf::calcWeights(), GFKalman::chi2Increment(), getBiasedSmoothedData(), GFDafHit::getHitCov(), getSmoothedChiSqu(), and getSmoothedData().

                                                                            {
        inv.ResizeTo(mat);
        bool status = 0;
        TDecompSVD invertAlgo(mat);

  // check if numerical limits are reached (i.e at least one entry < 1E-100 and/or at least one entry > 1E100)
        if (!(mat<1.E100) || !(mat>-1.E100)){
                GFException e("cannot invert matrix GFTools::invertMatrix(), entries too big (>1E100)",
                                __LINE__,__FILE__);
                e.setFatal();
                throw e;        
        }
        
        invertAlgo.SetTol(1E-50); //this is a hack because a tolerance of 1E-22 does not make any sense for doubles only for floats
        
        inv = invertAlgo.Invert(status);
        if(status == 0){
                GFException e("cannot invert matrix GFTools::invertMatrix(), status = 0",
                                __LINE__,__FILE__);
                e.setFatal();
                throw e;
        }
}