KalmanVertexFit Class Reference

#include <KalmanVertexFit.h>

Public Member Functions
	~KalmanVertexFit ()
void	init ()
void	initVertex (const VertexParameter vtx)
void	addTrack (const HTrackParameter)
int	filter (const int k)
void	filter ()
void	inverse (const int k)
void	smooth (const int k)
void	smooth ()
void	remove (const int k)
HepVector	pull (const int k)
double	pullmomentum (const int k)
VertexParameter	vtx () const
HepVector	x () const
HepSymMatrix	Ex () const
int	numTrack () const
int	ndof () const
double	chisq ()
int	trackID (const int k) const
std::vector< int >	trackIDList () const
double	chiF (const int k) const
double	chiS (const int k) const
HTrackParameter	hTrack (const int k) const
WTrackParameter	wTrack (const int k, const double mass) const
void	setVertexIteration (const int num)
void	setTrackIteration (const int num)
void	setChisqCut (const double chicut)
void	setTrackIterationCut (const double chicut)

Static Public Member Functions
static KalmanVertexFit *	instance ()

Detailed Description

Definition at line 16 of file KalmanVertexFit.h.

Constructor & Destructor Documentation

~KalmanVertexFit()

KalmanVertexFit::~KalmanVertexFit

(

)

Definition at line 25 of file KalmanVertexFit.cxx.

{;}

Member Function Documentation

addTrack()

void KalmanVertexFit::addTrack

(

const HTrackParameter

htrk

)

Definition at line 81 of file KalmanVertexFit.cxx.

                                                         {
 
  int ifail;
 
  HepVector p(3, 0);
  p = htrk.p();
  
  m_p.push_back(p);
  m_flag.push_back(0);
  m_chiF.push_back(999.);
 
  HepMatrix A(5, 3, 0);
  HepMatrix B(5, 3, 0);
  HepSymMatrix G(5, 0); 
 
  G = (htrk.eHel()).inverse(ifail);
 
  m_hTrkOrigin.push_back(htrk);
 
  //xum
  m_hTrkInfit.push_back(htrk); 
  //m_wTrkInfit.push_back(htrk.wTrack());
 
  m_G.push_back(G);
  m_A.push_back(A);
  m_B.push_back(B);
  
  HepSymMatrix W(3, 0);
  HepSymMatrix GB(5, 0);
  HepVector c(5, 0);
 
  m_W.push_back(W);
  m_GB.push_back(GB);
  m_c.push_back(c);
 
}

Referenced by BeamParams::execute(), and PrimaryVertex::execute().

chiF()

double KalmanVertexFit::chiF ( const int k ) const

inline

Definition at line 83 of file KalmanVertexFit.h.

{return m_chiF[k];}  // chi-square of filter for track k

Referenced by BeamParams::execute(), and PrimaryVertex::execute().

chiS()

double KalmanVertexFit::chiS

(

const int

k

)

const

Definition at line 263 of file KalmanVertexFit.cxx.

                                              {
  if(m_flag[k]==1) return 999;
  int ifail;
  HepVector x(3, 0);
  HepVector p(3, 0);
  HepSymMatrix C(3, 0);
  C = m_C - m_GB[k].similarity(m_A[k].T());
  x = C.inverse(ifail)*(m_C*m_x-m_A[k].T() * m_GB[k] * 
            (m_hTrkOrigin[k].hel()-m_c[k]));
  p = m_W[k]* m_B[k].T()*m_G[k]*
    (m_hTrkOrigin[k].hel()-m_c[k]-m_A[k]*m_x);
  HepSymMatrix chis(1, 0);
  chis = C.similarity((x-m_x).T())+
    m_G[k].similarity((m_hTrkOrigin[k].hel()-m_c[k]-
               m_A[k]*m_x-m_B[k]*p).T());
  return chis[0][0];
}

Referenced by BeamParams::execute(), PrimaryVertex::execute(), and filter().

chisq()

double KalmanVertexFit::chisq ( )

inline

Definition at line 76 of file KalmanVertexFit.h.

{return m_chisq;}

Referenced by BeamParams::execute(), PrimaryVertex::execute(), and filter().

Ex()

HepSymMatrix KalmanVertexFit::Ex

(

)

const

Definition at line 76 of file KalmanVertexFit.cxx.

                                       {
  int ifail;
  return m_C.inverse(ifail);
}

Referenced by PrimaryVertex::execute().

filter() [1/2]

void KalmanVertexFit::filter

(

)

Definition at line 221 of file KalmanVertexFit.cxx.

                             {
  
  for(int iter = 0; iter < m_maxVertexIteration; iter++) {
    //
    //  at start of each iteration, 
    //  the covariance matrix take the initial values: m_C0
    //  the vertex position take the update value by the last iteration
    //
    m_C = m_C0;   // covariance matrix of verter --> initial values
    m_chisq = 0;  // total chi-square of smooth
    m_ndof = -3;  // number of degrees of freedom
   
    for(int k = 0; k < m_hTrkOrigin.size(); k++) {
      if(m_flag[k]==1) continue;
      if(filter(k)==0) // filter track k
    m_ndof += 2;
 
    }  
    //
    //  check the chi-square of smooth, make decision to remove bad tracks
    //  user may remove more bad track by tight chi-square cut
    //  through the following interface 
    //  KalmanVertexFit::remove(k)
    //
    for(int k = 0; k < m_hTrkOrigin.size(); k++) {
      if(m_flag[k]==1) continue;
      double chi2 = chiS(k);
      if(chi2 < 0) { // remove track k from vertex
    remove(k);
    continue;
      }
      
      if(chi2 > m_chisqCutforTrack) { // remove track k from vertex
    remove(k);
    continue;
      } 
      m_chisq += chi2;
 
    }
  } // end of iteration
}

Referenced by filter().

filter() [2/2]

int KalmanVertexFit::filter

(

const int

k

)

Definition at line 135 of file KalmanVertexFit.cxx.

                                       {
  int ifail;
  double chisq[10];
  ////////////////////////////////////////////////
  // start iteration
  ///////////////////////////////////////////////
  HepVector pp(3, 0);
  HepVector xp(3, 0);
  HepSymMatrix CP(3, 0);
  xp = m_x;
  CP = m_C;
  pp = m_p[k];
 
  for(int iter = 0; iter < m_maxTrackIteration; iter++) {
    //
    // calculate derivative matrices for track k
    //
    updateMatrices(k, pp, xp);
    //
    // vertex covariance matrix updated
    //
    HepSymMatrix CK(3, 0);
    CK = CP + (m_GB[k].similarity(m_A[k].T()));
    //
    // state vector x and p updated
    //
    HepVector x(3, 0);
    HepVector p(3, 0);
    x = CK.inverse(ifail) * (CP * xp + m_A[k].T() * m_GB[k] * 
                 (m_hTrkOrigin[k].hel() - m_c[k]));
    p = m_W[k] * m_B[k].T() * m_G[k] *(m_hTrkOrigin[k].hel()- m_c[k]- 
                       m_A[k]*x);
    // chi-square of filter 
    HepSymMatrix chif(1, 0);
    chif = CP.similarity((x-xp).T()) + 
      m_G[k].similarity((m_hTrkOrigin[k].hel()- 
             m_c[k]- m_A[k]*x - m_B[k]*p).T());
    
    //
    // save current 3-momentum and vertex for track k
    //
    // 3-momentum
    pp = p;
    // current vertex and its covaiance matrix
    xp = x;
    CP = CK;
    chisq[iter] = chif[0][0];
    m_chiF[k] = chif[0][0];
    //    std::cout << "chisq, iter = " <<m_chiF[k] << ", " << iter << ", " << k <<std::endl;
    if(iter > 0) {
      double delchi = (chisq[iter]-chisq[iter-1]);
      if(fabs(delchi) < m_chisqCutforTrackIteration) break;
    }
  }
 
  if(m_chiF[k] > m_chisqCutforTrack) {
    //    do not update vertex if failed
    m_flag[k]=1;
    return 1;
  } else { 
    // update vertex if filter success
    m_x = xp;
    m_C = CP;
    m_p[k] = pp;
    return 0;
  }
}

Referenced by BeamParams::execute(), and PrimaryVertex::execute().

hTrack()

HTrackParameter KalmanVertexFit::hTrack

(

const int

k

)

const

Definition at line 381 of file KalmanVertexFit.cxx.

                                                        {
  //add xum
  //HTrackParameter htrk;
  //return htrk;
  return m_hTrkInfit[k];
}

init()

void KalmanVertexFit::init

(

)

Definition at line 29 of file KalmanVertexFit.cxx.

                           {
  m_x = HepVector(3, 0);
  m_C0 = HepSymMatrix(3, 0);
  m_C = HepSymMatrix(3, 0);
  m_chisq = 0;
  m_ndof = -3;
 
  m_flag.clear();
  m_p.clear();
  m_hTrkOrigin.clear();
  m_hTrkInfit.clear(); //xum
  //m_wTrkInfit.clear();
 
  m_G.clear();
  m_A.clear();
  m_B.clear();
  m_c.clear();
  m_W.clear();
  m_GB.clear();
  m_chiF.clear();
 
  //
  // chi-square cut, number of iteration etc
  //
  m_chisqCutforTrack = 50;
  m_maxVertexIteration = 3;
  m_maxTrackIteration = 5;
  m_chisqCutforTrackIteration = 0.1;
}

Referenced by BeamParams::execute(), and PrimaryVertex::execute().

initVertex()

void KalmanVertexFit::initVertex

(

const VertexParameter

vtx

)

Definition at line 59 of file KalmanVertexFit.cxx.

                                                          {
  int ifail;
  m_x = vtx.Vx();
  m_C0 = (vtx.Evx()).inverse(ifail);
  m_C = m_C0;
}

Referenced by BeamParams::execute(), and PrimaryVertex::execute().

instance()

KalmanVertexFit * KalmanVertexFit::instance ( )

static

Definition at line 19 of file KalmanVertexFit.cxx.

                                            {
  if(m_pointer) return m_pointer;
  m_pointer = new KalmanVertexFit();
  return m_pointer;
}

Referenced by BeamParams::execute(), and PrimaryVertex::execute().

inverse()

void KalmanVertexFit::inverse

(

const int

k

)

Definition at line 203 of file KalmanVertexFit.cxx.

                                         {
  //
  // inverse kalman filter for track k
  //
 
  int ifail;
  HepVector x(3, 0);
  HepVector p(3, 0);
  HepSymMatrix C(3, 0);
  C = m_C - m_GB[k].similarity(m_A[k].T());
  x = C.inverse(ifail)*(m_C*m_x-m_A[k].T() * m_GB[k] * 
            (m_hTrkOrigin[k].hel()-m_c[k]));
  //
  //  update vertex position and its covariance matrix
  //
  m_x = x; 
  m_C = C; 
}

Referenced by remove().

ndof()

int KalmanVertexFit::ndof ( ) const

inline

Definition at line 75 of file KalmanVertexFit.h.

{return m_ndof;}

Referenced by BeamParams::execute(), and PrimaryVertex::execute().

numTrack()

int KalmanVertexFit::numTrack

(

)

const

Definition at line 118 of file KalmanVertexFit.cxx.

                                    {
  int num = 0;
  for(int k = 0; k < m_hTrkOrigin.size(); k++) 
    if(m_flag[k] == 0) num = num +1;
  
  return num;
}

Referenced by PrimaryVertex::execute().

pull()

HepVector KalmanVertexFit::pull

(

const int

k

)

Definition at line 468 of file KalmanVertexFit.cxx.

                                           {
  HepVector n_pull(5, 0);
  n_pull[0] = n_pull[1] = n_pull[2] = n_pull[3] = n_pull[4] = 999.;
  if(m_flag[k] == 1) return n_pull; 
  for (int i = 0 ; i < 5; i++) {
    double cov = m_hTrkOrigin[k].eHel()[i][i] - m_hTrkInfit[k].eHel()[i][i];
    if (cov == 0.) continue;    
    n_pull[i] = (m_hTrkOrigin[k].hel()[i] - m_hTrkInfit[k].hel()[i]) / sqrt(cov);
  }
  return n_pull;
}

Referenced by BeamParams::execute(), PrimaryVertex::execute(), and pullmomentum().

pullmomentum()

double KalmanVertexFit::pullmomentum

(

const int

k

)

Definition at line 480 of file KalmanVertexFit.cxx.

                                                {
  double pull = 999.;
  if(m_flag[k] == 1) return pull;
 
  double kappa_origin = m_hTrkOrigin[k].kappa();
  double kappa_infit = m_hTrkInfit[k].kappa();
  double lamb_origin = m_hTrkOrigin[k].lambda();
  double lamb_infit = m_hTrkInfit[k].lambda();
 
  double Vkappa_origin = m_hTrkOrigin[k].eHel()[2][2];
  double Vkappa_infit = m_hTrkInfit[k].eHel()[2][2];
  double Vlamb_origin = m_hTrkOrigin[k].eHel()[4][4];
  double Vlamb_infit = m_hTrkInfit[k].eHel()[4][4];
  double V_kappa_lamb_origin = m_hTrkOrigin[k].eHel()[4][2];
  double V_kappa_lamb_infit = m_hTrkInfit[k].eHel()[4][2];
 
  double P_origin = calculationP(kappa_origin, lamb_origin);
  //cout << "P_origin = " << P_origin << endl;
  //P_origin = m_hTrkOrigin[k].p()[0] * m_hTrkOrigin[k].p()[0]
  //          + m_hTrkOrigin[k].p()[1] * m_hTrkOrigin[k].p()[1]
  //          + m_hTrkOrigin[k].p()[2] * m_hTrkOrigin[k].p()[2];
  //cout << "P = " << P_origin << endl;
  double P_infit = calculationP(kappa_infit, lamb_infit);
 
  double SigmaP_origin = calculationSigmaP(kappa_origin, lamb_origin, Vkappa_origin,
                        Vlamb_origin, V_kappa_lamb_origin);
  double SigmaP_infit = calculationSigmaP(kappa_infit, lamb_infit, Vkappa_infit,
                        Vlamb_infit, V_kappa_lamb_infit);
  if ((SigmaP_origin - SigmaP_infit) <= 0.)  return pull;
  pull = (P_origin - P_infit) / sqrt (SigmaP_origin - SigmaP_infit);
 
  return pull;
}

Referenced by BeamParams::execute(), and PrimaryVertex::execute().

remove()

void KalmanVertexFit::remove

(

const int

k

)

Definition at line 355 of file KalmanVertexFit.cxx.

                                        {
  //
  // remove track k from the vertex fit
  //
  /*
    m_hTrkOrigin.erase(m_hTrkOrigin.begin()+k);
    m_p.erase(m_p.begin()+k);
    m_flag.erase(m_flag.begin()+k);
    m_A.erase(m_A.begin()+k);
    m_B.erase(m_B.begin()+k);
    m_c.erase(m_c.begin()+k);
    m_G.erase(m_G.begin()+k);
    m_GB.erase(m_GB.begin()+k);
    m_W.erase(m_W.begin()+k);
    m_chiF.erase(m_chiF.begin()+k);
  */
 
  if(m_flag[k]==1) return; // track k already removed
  inverse(k);
  m_ndof -= 2;
  m_flag[k]=1;
  //xum
  //m_chisq -=chiS(k);
}

Referenced by BeamParams::execute(), PrimaryVertex::execute(), and filter().

setChisqCut()

void KalmanVertexFit::setChisqCut ( const double chicut )

inline

Definition at line 97 of file KalmanVertexFit.h.

{m_chisqCutforTrack = chicut;}

Referenced by BeamParams::execute(), and PrimaryVertex::execute().

setTrackIteration()

void KalmanVertexFit::setTrackIteration ( const int num )

inline

Definition at line 96 of file KalmanVertexFit.h.

{m_maxTrackIteration = num;}

Referenced by BeamParams::execute(), and PrimaryVertex::execute().

setTrackIterationCut()

void KalmanVertexFit::setTrackIterationCut ( const double chicut )

inline

Definition at line 98 of file KalmanVertexFit.h.

{m_chisqCutforTrackIteration = chicut;} //xum

Referenced by BeamParams::execute(), and PrimaryVertex::execute().

setVertexIteration()

void KalmanVertexFit::setVertexIteration ( const int num )

inline

Definition at line 95 of file KalmanVertexFit.h.

{m_maxVertexIteration = num;}

Referenced by BeamParams::execute(), and PrimaryVertex::execute().

smooth() [1/2]

void KalmanVertexFit::smooth

(

)

smooth() [2/2]

void KalmanVertexFit::smooth

(

const int

k

)

Definition at line 281 of file KalmanVertexFit.cxx.

                                        {
  // xum
  //
  // update htrk and wtrk parameters
  //
  if(m_flag[k] == 1) return;
  int ifail;
 
  HepVector pp(3, 0);
  HepVector xp(3, 0);
  HepSymMatrix CP(3, 0);
  xp = m_x;
  CP = m_C;
  pp = m_p[k];
 
  updateMatrices(k, pp, xp);
  pp = m_W[k] * m_B[k].T() * m_G[k] * (m_hTrkOrigin[k].hel() - m_c[k] - m_A[k] * xp);
  
  updateMatrices(k, pp, xp);
  //for htrk
  HepVector helix(5, 0);
  helix = m_c[k] + m_A[k] * xp + m_B[k] * pp;
  HepMatrix E(3, 3, 0);
  E = -CP.inverse(ifail) * m_A[k].T() * m_G[k] * m_B[k] * m_W[k];
  HepSymMatrix D(3, 0);
  D = m_W[k] + CP.similarity(E.T());
  
  HepMatrix middle(5, 5, 0);
  HepSymMatrix error(5, 0);
  middle = (CP.inverse(ifail)).similarity(m_A[k]) + m_A[k] * E * m_B[k].T() +
          (m_A[k] * E * m_B[k].T()).T() + D.similarity(m_B[k]);  
  error.assign(middle);
 
  m_hTrkInfit[k].setHel(helix);
  m_hTrkInfit[k].setEHel(error);
  
  m_p[k] = pp;
/*
  //for wtrk
  double mass = m_hTrkOrigin[k].xmass(m_hTrkOrigin[k].partID());
  double Energy = sqrt(pp[0]*pp[0] + pp[1]*pp[1] + pp[2]*pp[2] + mass * mass);
 
  HepMatrix Awtrk(7, 3, 0), Bwtrk(7, 3, 0);
  Awtrk[4][0] = Awtrk[5][1] = Awtrk[6][2] = 1;
  Bwtrk[0][0] = Bwtrk[1][1] = Bwtrk[2][2] = 1;
  Bwtrk[3][0] = pp[0] / Energy;
  Bwtrk[3][1] = pp[1] / Energy;
  Bwtrk[3][2] = pp[2] / Energy;
    
  HepVector w(7, 0);
  HepSymMatrix Ew(7, 0);
  w[0] = pp[0];   w[1] = pp[1];  w[2] = pp[2];  w[3] = Energy; 
  w[4] = xp[0];   w[5] = xp[1];  w[6] = xp[2];
  
  HepSymMatrix Gwtrk(7, 0);
  Gwtrk = m_hTrkOrigin[k].wTrack().Ew().inverse(ifail);
  HepSymMatrix Wwtrk(3, 0);
  Wwtrk = Gwtrk.similarity(Bwtrk.T()).inverse(ifail);
 
  HepMatrix Ewtrk(3, 3, 0);
  Ewtrk = -CP.inverse(ifail) * m_Awtrk[k].T() * m_Gwtrk[k] * m_Bwtrk[k] * m_Wwtrk[k];
  HepSymMatrix Dwtrk(3, 0);
  Dwtrk = m_Wwtrk[k] + CP.similarity(Ewtrk.T());
 
  HepMatrix Ewmiddle(7, 7, 0);
  Ewmiddle = (CP.inverse(ifail)).similarity(m_Awtrk[k]) + m_Awtrk[k] * Ewtrk * m_Bwtrk[k].T() +
          (m_Awtrk[k] * Ewtrk * m_Bwtrk[k].T()).T() + Dwtrk.similarity(m_Bwtrk[k]);         
  Ew.assign(Ewmiddle);
 
  m_wTrkInfit[k].setCharge(m_hTrkOrigin[k].charge());
  m_wTrkInfit[k].setW(w);
  m_wTrkInfit[k].setEw(Ew);*/
}

Referenced by BeamParams::execute(), and PrimaryVertex::execute().

trackID()

int KalmanVertexFit::trackID ( const int k ) const

inline

Definition at line 77 of file KalmanVertexFit.h.

{return m_hTrkOrigin[k].trackID();}

Referenced by trackIDList().

trackIDList()

std::vector< int > KalmanVertexFit::trackIDList

(

)

const

Definition at line 126 of file KalmanVertexFit.cxx.

                                                  {
  std::vector<int> trackid;
  for(int k = 0; k < m_hTrkOrigin.size(); k++)
    if(m_flag[k] == 0) trackid.push_back(trackID(k));
 
  return trackid;
}

Referenced by PrimaryVertex::execute().

vtx()

VertexParameter KalmanVertexFit::vtx

(

)

const

Definition at line 66 of file KalmanVertexFit.cxx.

                                           {
  int ifail;
  VertexParameter v;
  v.setVx(m_x);
 
  std::vector<int> trackid;
    v.setEvx(m_C.inverse(ifail));
  return v;
}

Referenced by initVertex().

wTrack()

WTrackParameter KalmanVertexFit::wTrack	(	const int	k,
		const double	mass ) const

Definition at line 388 of file KalmanVertexFit.cxx.

                                                                            {
  if(m_flag[k] == 1) return m_hTrkOrigin[k].wTrack();
  int ifail;
 
  HepVector pp(3, 0);
  HepVector xp(3, 0);
  HepSymMatrix CP(3, 0);
  xp = m_x;
  CP = m_C;
  pp = m_p[k];
 
  WTrackParameter wtrk;
  double Energy = sqrt(pp[0]*pp[0] + pp[1]*pp[1] + pp[2]*pp[2] + mass * mass);
  
  HepMatrix Awtrk(7, 3, 0), Bwtrk(7, 3, 0);
  Awtrk[4][0] = Awtrk[5][1] = Awtrk[6][2] = 1;
  Bwtrk[0][0] = Bwtrk[1][1] = Bwtrk[2][2] = 1;
  Bwtrk[3][0] = pp[0] / Energy;
  Bwtrk[3][1] = pp[1] / Energy;
  Bwtrk[3][2] = pp[2] / Energy;
    
  HepVector w(7, 0);
  HepSymMatrix Ew(7, 0);
  w[0] = pp[0];   w[1] = pp[1];  w[2] = pp[2];  w[3] = Energy;
  w[4] = xp[0];   w[5] = xp[1];  w[6] = xp[2];  
  
  HepSymMatrix Gwtrk(7, 0);
  Gwtrk = m_hTrkOrigin[k].wTrack().Ew().inverse(ifail);
  HepSymMatrix Wwtrk(3, 0);
  Wwtrk = Gwtrk.similarity(Bwtrk.T()).inverse(ifail);
  
  HepMatrix Ewtrk(3, 3, 0);
  Ewtrk = -CP.inverse(ifail) * Awtrk.T() * Gwtrk * Bwtrk * Wwtrk;
  HepSymMatrix Dwtrk(3, 0);
  Dwtrk = Wwtrk + CP.similarity(Ewtrk.T());
 
  HepMatrix Ewmiddle(7, 7, 0);
  Ewmiddle = (CP.inverse(ifail)).similarity(Awtrk) + Awtrk * Ewtrk * Bwtrk.T() +
          (Awtrk * Ewtrk * Bwtrk.T()).T() + Dwtrk.similarity(Bwtrk);
  Ew.assign(Ewmiddle);
 
  wtrk.setCharge(m_hTrkOrigin[k].charge());
  wtrk.setW(w);
  wtrk.setEw(Ew);
 
  return wtrk;
}

x()

HepVector KalmanVertexFit::x ( ) const

inline

Definition at line 68 of file KalmanVertexFit.h.

{return m_x;}

Referenced by chiS(), BeamParams::execute(), PrimaryVertex::execute(), filter(), and inverse().

---

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

• KalmanVertexFit.h
• KalmanVertexFit.cxx