#include <VtVertex.hh>

Inheritance diagram for VERTEX::Vertex:

[legend]

Collaboration diagram for VERTEX::Vertex:

[legend]

Public Member Functions
— Constructors —
	Vertex ()

	Vertex (const Vertex &rhs)

	Vertex (const Track_v &v)

	Vertex (Track &t1, Track &t2)

	~Vertex ()

— VertexIf access functions —
float	vx () const
	$x$ of vertex More...

float	vy () const
	$y$ of vertex More...

float	vz () const
	$z$ of vertex More...

float	chi2 () const
	$\chi^2$ of vertex fit More...

float	prob () const
	upper tail $\chi^2$ probability More...

unsigned short int	ndf () const
	degrees of freedom of vertex fit More...

unsigned short int	ntracks () const
	no of tracks in vertex $n$ More...

float	vtx_cov_x () const
	$\sigma_x^2$ of vertex More...

float	vtx_cov_y () const
	$\sigma_y^2$ of vertex More...

float	vtx_cov_z () const
	$\sigma_z^2$ of vertex More...

– Other access methods —
double	angle () const

double	dist () const
	available after call to rmsDistAngle() More...

double	vxerr () const
	$\sqrt{\sigma_{vx}^2}$ vertex $x$-error More...

double	vyerr () const
	$\sqrt{\sigma_{vy}^2}$ vertex $y$-error More...

double	vzerr () const
	$\sqrt{\sigma_{vz}^2}$ vertex $z$-error More...

double	pmaxfrac () const
	$\max\|p_i/p_j\|$ max. momentum fraction More...

const MATRIX::VtSymMatrix &	CS () const
	$3\times3$ Vertex covariance matrix More...

const MATRIX::VtSymMatrix &	covn () const
	$(3\cdot n+3)\times(3\cdot n+3)$ general covariance matrix More...

bool	valid () const
	is vertex valid? More...

void	set_invalid ()
	mark vertex as invalid More...

— Mass calculation routines —
@memo $m = \sqrt{\sum_i p_i^2 + m_i^2 - \sum_i \vec{p}_i^2}$
double	mass (const bool use=false) const
	rest-mass 0 or track rest mass More...

double	massCC (const bool use=false) const
	rest-mass 0 or conjugated track rest mass More...

double	mass (double m1) const
	single rest-mass calc., More...

double	mass (double m1, double m2) const
	2-track mass with particle rest-masses More...

double	mass (double m1, double m2, double m3) const
	3-track mass with rest-masses More...

double	mass (double m1, double m2, double m3, double m4) const
	4-track mass with rest-masses More...

double	mass (const Vector_d &m02, bool use=false, bool CC=false) const
	specify $n$ rest masses More...

— Vertexing routines —
const bool	findVertex2D (void)
	unweighted 2-D min. dist. More...

const bool	findVertex3D (void)
	unweighted 3-D min. dist. More...

const bool	findVertexVt ()
	Vt Kalman-filter. More...

const bool	VtEstimateVertex ()
	weighted 2-D min. dist. More...

const bool	VtEstimateVertexMath (double &x, double &y, double &z)
	estimate Vertex without changing Vertex object More...

const bool	VtEstimateVertexMathTA (double &x, double &y, double &z)
	estimate Vertex without changing Vertex object More...

const bool	VtMass ()
	Kalman filter with mass constraints. More...

— Adding tracks —
void	push_back (Track &track)

— Misc functions —
double	rmsDistAngle () const
	calc rms dist and rms angle More...

bool	calc_mother (bool use_kalman=true)
	calc track params of mother particle More...

bool	calc_mother_cov ()
	calc cov. matrix of mother particle More...

const bool	remove_last ()
	remove last track from vertex, do inverse Kalman filter More...

void	remove_worst ()
	remove track with biggest $\chi^2$ contribution, do inverse Kalman filter More...

double	track_chi2 (int i)
	get track $\chi^2$ contribution More...

int	track_worst ()
	get track with biggest $\chi^2$ contribution More...

void	use_momentum (const bool use)
	set for all tracks whether momentum should be used or not More...

void	use_kalman (const bool use)
	use refitted track parameters or not More...

bool	use_kalman () const

double	distance (double x, double y, double z) const
	$\chi^2$ distance to space point, $ndf = 3$ More...

double	distance (Track &t) const
	$\chi^2$ distance to track, $ndf = 2$ More...

double	distance (const Vertex &rhs) const
	$\chi^2$ distance to vertex, $ndf = 3$ More...

— Expert routines —
void	clear ()
	clear all track-vertex relations, makes vertex invalid More...

const bool	VtFilter ()

const bool	VtInverseFilter () const

void	VtSmoothX ()
	smooth vertex position More...

void	VtSmoothQ ()
	smooth momenta More...

const bool	VtRemoveTrack (Relation &track)
	remove track using inverse Kalman filter More...

const bool	VtFit ()
	Kalman filter + smoother. More...

const MATRIX::VtVector &	xv () const

const MATRIX::VtVector &	xvs () const
	vector of vertex-positions More...

std::ostream &	print (std::ostream &os) const
	called by operator<<() More...

double	chi2n () const
	$\chi^2$ after last filter step More...

double	chi2l () const
	$\chi^2$ contribution of last filter step More...

virtual void	remove (Relation *r)
	unbook relation More...

virtual const iterator	erase (const iterator &pos)
	erase a relation without refit More...

const MATRIX::VtSymMatrix &	bigcov ()
	compute general $(3\cdot n+3) \times (3\cdot n+3)$ covariance matrix More...

— Mass Constraint routines —
void	MassConstr (double m)
	add a mass constraint for all tracks, $m$ = mass of mother particle More...

MassC &	addMassConstr (double m=0.)
	add a mass constraint More...

void	clearMassConstr ()
	delete the mass constraints More...

unsigned int	nMassConstr () const
	number of mass constraints More...

Public Member Functions inherited from VERTEX::Track
	Track ()

	Track (const Track *const track)

	Track (const Track &rhs)

	Track (const MATRIX::VtVector &v, const MATRIX::CMatrix &c)

virtual	~Track ()

void	set (double x, double y, double z, double tx, double ty, double p, const MATRIX::CMatrix &c)

Track &	operator= (const Track &rhs)

float	x () const
	x at z=0 (in Vt: p(3,..)) More...

float	y () const
	y at z=0 (in Vt: p(4,..)) More...

float	x (float z) const
	Track $x$ position at $z$. More...

float	y (float z) const
	Track $y$ position at $z$. More...

float	z () const
	z = 0 in Vt More...

float	tx () const
	slope (in Vt: p(1,..)) More...

float	ty () const
	slope (in Vt: p(2,..)) More...

float	p () const
	momentum (in Vt: p(5,..)) More...

float	pt () const
	transv. momentum $\sqrt{p_x^2 + p_y^2}$ More...

float	pz () const
	$p \cdot e_z$ More...

float	chi2 () const
	dummy function: always return 0 More...

unsigned short int	ndf () const
	dummy function: always return 0 More...

float	phi () const
	azimuthal angle $\phi$ More...

float	theta () const
	polar angle $\theta = \cos^{-1}(e_z)$ More...

float	eta () const
	rapidity $-\log\tan(\theta/2.)$ More...

int	charge () const
	charge More...

float	energy (double rm=0.) const
	Energy with given rest-mass (in GeV) $E = \sqrt{p^2 + m^2}$. More...

float	xf (double rm=0.) const
	Feynman Variable $x_F = \frac{E+p_z}{(E+p_z)_{max}}$. More...

float	rap (double rm=0.) const
	Rapidity $y = \frac{1}{2}\ln(\frac{E+p_z}{E-p_z})$. More...

MATRIX::VtVector	evec () const
	$\vec{v} = (e_x,e_y,e_z)$ unit vector along refitted track More...

MATRIX::VtVector	tvec () const
	$\vec{v} = (t_x,t_y,1.)$ refitted slope vector More...

MATRIX::VtVector	pvec () const
	$\vec{v} = (p_x,p_y,p_z)$ refitted mom. vector More...

float	cov_x (double dz=0.) const
	get $\sigma_x^2\|_{z+dz}$ More...

float	cov_y (double dz=0.) const
	get $\sigma_y^2\|_{z+dz}$ More...

float	cov_tx () const
	get $\sigma_{t_x}^2$ More...

float	cov_ty () const
	get $\sigma_{t_y}^2$ More...

float	cov_p () const
	get $\sigma_p^2$ More...

bool	isValid () const
	dummy function: always return true More...

void	valid ()
	dummy function: do nothing More...

void	invalid ()
	dummy function: do nothing More...

bool	propagate (const double zz)
	propagate track to $z$ More...

double	px () const
	$p \cdot e_x$ More...

double	py () const
	$p \cdot e_y$ More...

double	ex () const
	$t_x \cdot e_z$ More...

double	ey () const
	$t_y \cdot e_z$ More...

double	ez () const
	$1/\sqrt{1+t_x^2+t_y^2}$ More...

double	xerr (double dz=0) const
	$\sqrt{\sigma_x^2(z+dz)}$ More...

double	yerr (double dz=0) const
	$\sqrt{\sigma_y^2(z+dz)}$ More...

double	txerr () const
	$\sqrt{\sigma^2_{tx}}$ More...

double	tyerr () const
	$\sqrt{\sigma^2_{ty}}$ More...

double	perr () const
	$\sqrt{\sigma^2_p}$ More...

const MATRIX::CMatrix &	V () const
	covariance matrix More...

const MATRIX::CMatrix &	G () const
	inverse $5\times5$ covariance matrix More...

const MATRIX::CMatrix &	GM () const
	inverse $4\times4$ covariance matrix, without momentum More...

std::ostream &	print (std::ostream &os) const
	called by cout More...

double	rm () const
	rest-mass More...

void	rm (const double mass)
	set rest-mass (needed for mass constrained fits) More...

double	rmCC () const
	conjugated rest-mass More...

void	rmCC (const double mass)
	set conjugated rest-mass More...

void	delete_mom ()

Public Member Functions inherited from VERTEX::RelationList
	RelationList ()

virtual	~RelationList ()

void	push__back (Relation *rel)
	add a relation More...

Relation *	back () const
	return last relation More...

unsigned int	size () const
	return no of relations More...

iterator	begin ()

iterator	end ()

const_iterator	begin () const

const_iterator	end () const

reverse_iterator	rbegin ()

reverse_iterator	rend ()

const_reverse_iterator	rbegin () const

const_reverse_iterator	rend () const

Relation *const	operator[] (unsigned int i) const
	direct access More...

void	print (std::ostream &os) const

void	unbook (Relation *const rel)
	get rid of relation pointer More...

const RelationList &	operator= (const RelationList &rhs)

bool	operator== (const RelationList &rhs) const

— Operators —
bool	v_use_kalman

bool	v_valid

bool	v_Mvalid

bool	v_mother

bool	v_angdist

double	v_angle

double	v_dist

double	v_bk13

double	v_chi2

MATRIX::VtVector	kal_xv

MATRIX::VtVector	kal_xvs

MATRIX::VtSymMatrix	v_CS

MATRIX::VtSymMatrix	v_CINV

MATRIX::VtSymMatrix *	v_covn

Track_v	tracks

MassC_v	xmass

const Vertex &	operator= (const Vertex &rhs)

const bool	operator== (const Vertex &rhs) const

const double	operator- (const Vertex &rhs) const

void	add_track (Track &t)

const Track *	get_track (int i) const

const unsigned int	bigdim () const

void	invalid ()

	ClassDef (Vertex, 1)

Additional Inherited Members
Public Types inherited from VERTEX::RelationList
typedef ConstRelationIterator	const_iterator

typedef ConstReverseRelationIterator	const_reverse_iterator

typedef RelationIterator	iterator

typedef ReverseRelationIterator	reverse_iterator

Protected Member Functions inherited from VERTEX::Track
	ClassDef (Track, 1)

Protected Attributes inherited from VERTEX::Track
MATRIX::VtVector	t_p

MATRIX::CMatrix	t_V

MATRIX::CMatrix	t_G

MATRIX::CMatrix	t_GM

Detailed Description

Vertex class

Constructor & Destructor Documentation

◆ Vertex() [1/4]

Vertex::Vertex ( )

                 :
    Track       (),
    v_use_kalman(false),
    v_valid     (false),
    v_Mvalid    (false),
    v_mother    (false),
    v_angdist   (false),
    v_angle     (0.),
    v_dist      (0.),
    v_bk13      (0.),
    v_chi2      (0.),
    kal_xv      (VtVector(3)),
    kal_xvs     (VtVector(3)),
    v_CS        (VtSymMatrix(3,1.e5)),
    v_CINV      (CINV0),
    v_covn      (0),
    tracks      (),
    xmass       ()
  {}

◆ Vertex() [2/4]

Vertex::Vertex ( const Vertex & rhs )

                                 :
    Track       (rhs),
    v_use_kalman(rhs.v_use_kalman),
    v_valid     (rhs.v_valid),
    v_Mvalid    (rhs.v_Mvalid),
    v_mother    (rhs.v_mother),
    v_angdist   (rhs.v_angdist),
    v_angle     (rhs.v_angle),
    v_dist      (rhs.v_dist),
    v_bk13      (rhs.v_bk13),
    v_chi2      (rhs.v_chi2),
    kal_xv      (rhs.kal_xv),
    kal_xvs     (rhs.kal_xvs),
    v_CS        (rhs.v_CS),
    v_CINV      (v_CS.dsinv()),
    v_covn      (0)
    //    tracks  (rhs.tracks),
    //    xmass   (rhs.xmass)
  {
    // do not copy relations
  }

◆ Vertex() [3/4]

VERTEX::Vertex::Vertex ( const Track_v & v )

◆ Vertex() [4/4]

Vertex::Vertex	(	Track &	t1,
		Track &	t2
	)

                                    :
    Track       (),
    v_use_kalman(false),
    v_valid     (false),
    v_Mvalid    (false),
    v_mother    (false),
    v_angdist   (false),
    v_angle     (0.),
    v_dist      (0.),
    v_bk13      (0.),
    v_chi2      (0.),
    kal_xv      (VtVector(3)),
    kal_xvs     (VtVector(3)),
    v_CS        (VtSymMatrix(3,1.e5)),
    v_CINV      (CINV0),
    v_covn      (0),
    tracks      (),
    xmass       ()
 {
    createRelation(t1, *this);
    createRelation(t2, *this);
  }

◆ ~Vertex()

Vertex::~Vertex ( )

                  {
    // delete the mass constraints
    clearMassConstr();
    // delete tracks owned by vertex
    const Track_cit end = tracks.end();
    for(Track_cit it = tracks.begin(); it != end; ++it) {
      delete *it;
    }
    // delete big covariance matrix
    if(v_covn != 0) delete v_covn;
  }

Member Function Documentation

◆ add_track()

void Vertex::add_track ( Track & t )

                                 {
    tracks.push_back(&t);   // add track to the tracks list in a way that vertex become the owner
    push_back(t);           // set the relation
  }

◆ addMassConstr()

MassC & Vertex::addMassConstr ( double m = 0. )

add a mass constraint

                                       {
    xmass.push_back(new MassC(m));
    return *xmass.back();
  }

◆ angle()

double VERTEX::Vertex::angle ( ) const

$\sqrt{n\cdot\sum_it_{x,i}^2 - (\sum_i t_{x,i})^2 + n\cdot\sum_i t_{y,i}^2 - (\sum_i t_{y,i})^2} / n$

◆ bigcov()

const VtSymMatrix & Vertex::bigcov ( )

compute general $(3\cdot n+3) \times (3\cdot n+3)$ covariance matrix

                                    {
 
    // build up big cov. matrix (in Vt: covn)
    // covn consists of vertex + track parameter cov. matrices:
    //  vtx  tr1  tr2  tr3  ...
    //
    //  CS | ES | ES | ES | ...   CC = cor(q_k,q_l) = F_k*C_n*F_l^T
    //  ---+----+----+----+          = -ES_k * F_l
    //  ES | DS | CC | CC | ...
    //  ---+----+----+----+ 
    //  ES | CC | DS | CC | ...
    //  ---+----+----+----+
    //  ES | CC | CC | DS | ...
    //  ...
 
    if(v_covn == 0)
      v_covn = new VtSymMatrix(bigdim());
 
    VtSymMatrix& covn = *v_covn;
 
    //    VtSymMatrix covn(bigdim());
 
    covn.place_at(v_CS,0,0);
 
    unsigned int ipos = 3;
    const iterator last = end();
    for(iterator it = begin(); it != last; ++it) {
      const Kalman& kalman = it->kalman;
 
      covn.place_at(kalman.ES(), ipos, 0);
      covn.place_at(kalman.DS(), ipos, ipos);
 
      // fill cor(q_k^n,q_l^n) = F_k*C_n*F_l^T = -E_k^n * F_l^T
      unsigned int jpos = 6;
      iterator tmp=it;
      for(iterator it2=++tmp; it2 != last; ++it2) {
    const Kalman& kalman2 = it2->kalman;
    
    covn.place_at(-kalman.ES()*kalman2.F().T(), ipos,jpos);
    jpos+=3;
      }
 
      ipos+=3;
    }
 
    return covn;
  }

◆ bigdim()

const unsigned int VERTEX::Vertex::bigdim ( ) const

private

◆ calc_mother()

bool Vertex::calc_mother ( bool use_kalman = true )

calc track params of mother particle

                                          {
    //    if(v_valid == false) return v_mother = false;
 
    // add up track momenta
    t_p[5] = 0;
#if defined USE_ROOT
    SVector<double,3> psum;
#else
    VtVector psum(0,0,0);
#endif
 
    const iterator last = end();
    if(use_kalman == true) {
      for(iterator it = begin(); it != last; ++it) {
    const Kalman& t = it->kalman;
    t_p[5] += t.p();
    psum   += t.pvec();
      }
 
      if(fabs(psum[2])>1e-20) {
    t_p[3] = psum[0]/psum[2]; // tx = pxsum/pzsum
    t_p[4] = psum[1]/psum[2]; // ty = pysum/pzsum
      } else {
    return v_mother = false;
      }
 
      // x,y at z = 0
      t_p[0] = vx() - vz() * tx();
      t_p[1] = vy() - vz() * ty();
      t_p[2] = 0.;
 
    } else {
      double x, y, z;
      VtEstimateVertexMath(x,y,z);
      for(iterator it = begin(); it != last; ++it) {
    const Track& t = it->track;
    t_p[5] += t.p();
    psum   += t.pvec();
      }
 
      if(fabs(psum[2])>1e-20) {
    t_p[3] = psum[0]/psum[2]; // tx = pxsum/pzsum
    t_p[4] = psum[1]/psum[2]; // ty = pysum/pzsum
      } else {
    return v_mother = false;
      }
 
      // x,y at z = 0
      t_p[0] = x - z * tx();
      t_p[1] = y - z * ty();
      t_p[2] = 0.;
    }
 
    //    t_p[5] = sqrt(sqr(pxsum) + sqr(pysum) + sqr(pzsum));
 
    // what is the difference betwenn sum of track.p() and
    // calculation sqrt(sqr(pxsum)+sqr(pysum)+sqr(pzsum))?????
    //      cout << " reco-p: " << sqrt(sqr(pxsum) + sqr(pysum) + sqr(pzsum))
    //       << " orig reco p: " << p() << endl;
 
    return v_mother = true;
  }

◆ calc_mother_cov()

bool Vertex::calc_mother_cov ( )

calc cov. matrix of mother particle

                               {
    // no track parameters computed
    if(v_mother == false || v_valid == false) return false;
    // get big covariance matrix
    const VtSymMatrix& covn = (v_Mvalid==false)? bigcov() : *v_covn;
 
    // fill jacobian
    VtMatrix    jac(bigdim(),5);
 
    // unit vectors
    double epx = ex();
    double epy = ey();
    double epz = ez();
 
    // vertex part:
    jac(0,0) = jac(1,1) = 1.;
    jac(2,0) = -tx();
    jac(2,1) = -ty();
 
    double fact3 = 1./pz();
    //    double fact3 = 1./(p() * epz);
 
    unsigned int i = 2;
    const iterator last = end();
    for(iterator it = begin(); it != last; ++it) {
      const Kalman& t = it->kalman;
      
      // unit vectors along track:
      double ex = t.ex();
      double ey = t.ey();
      double ez = t.ez();
 
      double fact2 = t.pz();
      double fact  = fact2 * fact3;
 
      ++i;
      jac(i,2) = fact * (1. - sqr(ex) + tx() * ex * ez);
      jac(i,3) = fact * ex * (ty()*ez - ey);
      jac(i,0) = -vz() * jac(i,2);
      jac(i,1) = -vz() * jac(i,3);
 
      jac(i,4) = fact2 * ( epx*(1.-sqr(ex)) - epy*ex*ey - epz*ex*ez );
      
      ++i;
      jac(i,2) = fact*ey*(tx()*ez-ex);
      jac(i,3) = fact*(1. - sqr(ey) + ty()*ey*ez);
      jac(i,0) = -vz() * jac(i,2);
      jac(i,1) = -vz() * jac(i,3);
 
      jac(i,4) = fact2 * ( epy*(1.-sqr(ey)) - epx*ex*ey - epz*ey*ez );
 
      ++i;
      jac(i,2) = fact3 * (ex - tx()*ez);
      jac(i,3) = fact3 * (ey - ty()*ez);
      jac(i,0) = -vz() * jac(i,2);
      jac(i,1) = -vz() * jac(i,3);
      
      jac(i,4) = epx*ex + epy*ey + epz*ez;
    }
 
#ifdef VtDEBUG
    cout << " jac: " << jac << endl;
#endif
 
    // compute Track covariance matrix:
    t_V = covn.product(jac);
    // compute inverse
    t_G = t_V;
    t_G.VtDsinv();
    t_GM.copy(t_V.dsinv(4)); // inverse without momentum
 
    return true;
  }

◆ chi2()

float Vertex::chi2 ( ) const

$\chi^2$ of vertex fit

236{ return v_chi2; }

◆ chi2l()

double Vertex::chi2l ( ) const

$\chi^2$ contribution of last filter step

                             { 
    const_reverse_iterator it1;
    const_reverse_iterator it2 = it1 = rbegin();
    ++it2;
    return it1->kalman.chi2() - it2->kalman.chi2();
  }

◆ chi2n()

double Vertex::chi2n ( ) const

$\chi^2$ after last filter step

                             { 
    return back()->kalman.chi2();
  }

◆ ClassDef()

VERTEX::Vertex::ClassDef	(	Vertex	,
		1
	)

private

◆ clear()

void Vertex::clear ( void )

virtual

clear all track-vertex relations, makes vertex invalid

Reimplemented from VERTEX::RelationList.

                     {
    RelationList::clear();
    invalid();
    v_mother = v_use_kalman = false;
    // delete the mass constraints
    clearMassConstr();
    // delete tracks owned by vertex
    const Track_cit end = tracks.end();
    for(Track_cit it = tracks.begin(); it != end; ++it) {
      delete *it;
    }
    tracks.clear();
    // delete big covariance matrix
    if(v_covn != 0) { 
      delete v_covn; 
      v_covn = 0;
    }
  }

◆ clearMassConstr()

void Vertex::clearMassConstr ( )

delete the mass constraints

                               {
    // delete mass constraints owned by vertex
    const MassC_cit end = xmass.end();
    for(MassC_cit it = xmass.begin(); it != end; ++it) {
      delete *it;
    }
    xmass.clear();
    return;
  }

◆ covn()

const MATRIX::VtSymMatrix & VERTEX::Vertex::covn ( ) const

$(3\cdot n+3)\times(3\cdot n+3)$ general covariance matrix

◆ CS()

const MATRIX::VtSymMatrix & VERTEX::Vertex::CS ( ) const

$3\times3$ Vertex covariance matrix

◆ dist()

double VERTEX::Vertex::dist ( ) const

available after call to rmsDistAngle()

◆ distance() [1/3]

double Vertex::distance ( const Vertex & rhs ) const

$\chi^2$ distance to vertex, $ndf = 3$

                                                 {
    const VtSymMatrix& cov1 = back()->kalman.CINV();
    const VtSymMatrix& cov2 = rhs.back()->kalman.CINV();
    VtSymMatrix Cmat = cov1+cov2;
 
    if(Cmat.VtDsinv() == false) return 1000.;
 
    const VtVector x1  = cov1 * kal_xvs;
    const VtVector x2  = cov2 * rhs.kal_xvs;
    const VtVector xav = Cmat * (x1+x2);
    const VtVector dx1 = xav - kal_xvs;
    const VtVector dx2 = xav - rhs.kal_xvs;
 
    // get a chi2 value
    double sum1 = cov1.product(dx1);
    double sum2 = cov2.product(dx2);
    return (sum1 + sum2);
  }

◆ distance() [2/3]

double Vertex::distance	(	double	x,
		double	y,
		double	z
	)		const

$\chi^2$ distance to space point, $ndf = 3$

                                                               {
    VtVector dx(xx-vx(),yy-vy(),zz-vz());
    return back()->kalman.CINV().product(dx);
  }

◆ distance() [3/3]

double Vertex::distance ( Track & t ) const

$\chi^2$ distance to track, $ndf = 2$

                                        {
    Vertex tmp;
    const const_iterator last = end();
    for(RelationList::const_iterator it = begin(); it != last; ++it) {
        Track& tp = it->track;
    tmp.push_back(tp);
    }
    tmp.push_back(t);
    if (tmp.findVertexVt() == false) return -10.;
    return tmp.track_chi2(tmp.ntracks()-1);
  }

◆ erase()

const iterator Vertex::erase ( const iterator & pos )

virtual

erase a relation without refit

Reimplemented from VERTEX::RelationList.

                                                  {
    invalid();
    return RelationList::erase(pos);
  }

◆ findVertex2D()

const bool Vertex::findVertex2D ( void )

unweighted 2-D min. dist.

                                      {  
 
    unsigned int v_nrseg = size();
    if(v_nrseg < 2) return false;
    
    double sumtx=0, sumty=0, sumrx=0, sumry=0, sumrt=0, sumtt=0;
 
    double tx, ty, rx, ry;
    const RelationList::iterator last = end();
    for(RelationList::iterator pos = begin(); pos != last; ++pos) {
      const Track& track = pos->track;
 
      // slopes and position at z=0
      tx = track.tx();  rx = track.x() - tx*track.z();
      ty = track.ty();  ry = track.y() - ty*track.z();
      
      // update cumulants
      sumtx += tx;   sumrx += rx;
      sumty += ty;   sumry += ry;
      sumrt += rx*tx   + ry*ty;
      sumtt += sqr(tx) + sqr(ty);
    }
 
#ifdef VtDEBUG
    cout << " sumtx: " << sumtx << " sumty: " << sumty
     << " sumrx: " << sumrx << " sumry: " << sumry
     << " sumrt: " << sumrt << " sumtt: " << sumtt
     << endl;
#endif
 
    // get vertex z-position
 
    double z = -(v_nrseg*sumrt - sumrx*sumtx - sumry*sumty) / 
                (v_nrseg*sumtt - sqr<double>(sumtx) - sqr<double>(sumty));
  
    // and the transverse coordinates
    double x = (sumrx + z*sumtx)/static_cast<double>(v_nrseg);
    double y = (sumry + z*sumty)/static_cast<double>(v_nrseg);
 
    // fill vertex parameters
    kal_xvs[0]  = x;         // x-vertex
    kal_xvs[1]  = y;         // y-vertex
    kal_xvs[2]  = z;         // z-vertex
 
    v_valid = true;
    return true;
  }

◆ findVertex3D()

const bool Vertex::findVertex3D ( void )

unweighted 3-D min. dist.

                                      {
 
    if(size() < 2) return false;
 
    double c[3], v[3], tt[9];
    for(int i=0; i<9; i++) tt[i] = 0.;
    for(int i=0; i<3; i++) c[i]  = v[i] = 0.;
 
    double t[3], r[3], dt;
    const RelationList::iterator last = end();
    for(RelationList::iterator pos = begin(); pos != last; ++pos) {
      const Track& track = pos->track;
 
      t[0] = track.tx();  r[0] = track.x();
      t[1] = track.ty();  r[1] = track.y();
      t[2] = 1;           r[2] = track.z();
      
      double tsq = sqr(t[0]) + sqr(t[1]) + sqr(t[2]); 
      
      for(int i=0; i<3; i++) {
    for(int j=0; j<3; j++) {
      dt = t[i]*t[j]/tsq; if(i==j) dt -= 1; 
      c[i]      += dt*r[j];
      tt[3*i+j] += dt;
    }
      }
    }
    
    int ifail = 0, ndim = 3;
    double work[3];
    Dinv1(&ndim, tt, &ndim, work, &ifail);
    
    if(ifail) {
      cout << "findVertex: matrix inversion failed" << endl;
      return false;
    }
    
    // get the vertex position
    for(int i=0; i<3; i++) {
      for(int j=0; j<3; j++) v[i] += tt[3*i+j]*c[j];
    }
    printf("findVertex3D v: %f %f %f\n",v[0],v[1],v[2]);
 
    // fill vertex parameters
    kal_xvs[0]  = v[0];      // x-vertex
    kal_xvs[1]  = v[1];      // y-vertex
    kal_xvs[2]  = v[2];      // z-vertex
 
    // require that vertex-z is smaller than smallest hit-z of rsegs
    // to avoid vertexing the downstream end of rseg's
    //    for(RelationList::iterator it = begin(); it != last; ++it) {
    //      if(vz() > it->track.z() ) return false;
    //    }
 
    v_valid = true;
    return true;
  }

◆ findVertexVt()

const bool Vertex::findVertexVt ( )

Vt Kalman-filter.

                                  {  
    v_valid = v_Mvalid = false;
    if(ntracks() < 2) { return false; }
    
    // 1st step: find approximate Vertex from average of all given tracks.
    // Initialize Kalman-filter
    VtEstimateVertex();
    // kal_xvs[0] = vest[0];
    // kal_xvs[1] = vest[1];
    // kal_xvs[2] = vest[2];
    
    // 2nd step: initialize kalman filter
    v_CS   = VtSymMatrix(3,1.e5);
    v_CINV = CINV0;
    const RelationList::iterator last = end();
    for(RelationList::iterator tit = begin(); tit != last; ++tit) {
      tit->kalman.init();
    }
    
#ifdef VtDEBUG
    cout << "-------------------------------" << endl;
    cout << "Estimated Vertex: (" << vx() << "," << vy() << "," << vz() << ")" << endl;
    cout << "Number of Tracks: " << size();
    //copy(begin(), end(), ostream_iterator<Relation>(cout," "));
    cout << endl;
#endif
 
    double prChi2 = 0.;
    const Kalman& trackn = back()->kalman; // last element
 
    // 3rd step: iterate the kalman filter
    unsigned int iter = 0;
    do {
 
      if(VtFit() == false) return false;
 
      double chi2 = trackn.chi2();
 
#ifdef VtDEBUG
      cout << "----> chi2n: " << chi2 << " prChi2: " << prChi2 << " fabs(chi2-prChi2): " 
       << fabs(chi2 - prChi2) << endl;
#endif
 
      if(fabs(chi2 - prChi2) > 1.) {
    prChi2 = chi2;
      } else
    if(iter > 0)
      break;
    } while(++iter < 20);
 
#ifdef VtDEBUG
    cout << "-------------------------------" << endl;
    cout << "No of iterations: " << iter << endl;
    cout << *this;
#endif
 
    if(iter<20) { //  && chi2()<10.
      v_valid = true;
    } else {
      //      cout << "-----> VERTEX FIT FAILED!! <--------" << endl;
    }
    return v_valid;
  }

◆ get_track()

const Track * Vertex::get_track ( int i ) const

266{return tracks.at(i);}

◆ invalid()

void Vertex::invalid ( )

private

                       {
    v_valid = v_Mvalid = false;
 
    if(v_covn != 0) {
      delete v_covn;
      v_covn = 0;
    }
 
    return;
  }

◆ mass() [1/6]

double Vertex::mass ( const bool use = false ) const

rest-mass 0 or track rest mass

                                          {
    return mass(vector<double>(ntracks(),0), use, false);
  }

◆ mass() [2/6]

double Vertex::mass	(	const Vector_d &	m02,
		bool	use = `false`,
		bool	CC = `false`
	)		const

specify $n$ rest masses

                                                                  {
 
    double esum  = 0.;
#if defined USE_ROOT
    SVector<double,3> psum;
#else
    VtVector psum(0,0,0);
#endif
 
    unsigned int i=0;
    const const_iterator last = end();
    for(RelationList::const_iterator it = begin(); it != last; ++it) {
      const Track& t = it->track;
      double p = 0.;
      if(v_use_kalman == false && use == false) {
    p     = t.p();
    psum += t.pvec();
      } else {
    const Kalman& k = it->kalman;
    p     = k.p();
    psum += k.pvec();
      }
 
      // in case no rest-mass given, use permanent rest-mass
      double mm = m02[i++];
 
      // handle track rest masses
      if(mm == 0.) {
    double rm = t.rm();
    if(CC == true) rm = t.rmCC(); // use conjugated mass
    mm = rm * rm;
      }
      esum  += sqrt(p*p + mm);
    }
    //   cout << "unllllock mutex" << endl;
    //    pthread_mutex_unlock(&VTmutex);
 
#if defined USE_ROOT
    double mass2 = esum*esum - mag2(psum);
#else
    double mass2 = esum*esum - psum(0)*psum(0) - psum(1)*psum(1) -
      psum(2)*psum(2);
#endif
    if(mass2 < 0) return 0.;
    else return sqrt(mass2);
  }

◆ mass() [3/6]

double Vertex::mass ( double m1 ) const

single rest-mass calc.,

                                     {
    return mass(vector<double>(ntracks(),m1*m1));
  }

◆ mass() [4/6]

double Vertex::mass	(	double	m1,
		double	m2
	)		const

2-track mass with particle rest-masses

                                                {
    if(ntracks() != 2) return -1;
    vector<double> tmp(2);
    tmp[0] = m1*m1, tmp[1] = m2*m2;
    return mass(tmp);
  }

◆ mass() [5/6]

double Vertex::mass	(	double	m1,
		double	m2,
		double	m3
	)		const

3-track mass with rest-masses

                                                           {
    if(ntracks() != 3) return -1;
    vector<double> tmp(3);
    tmp[0] = m1*m1, tmp[1] = m2*m2, tmp[2] = m3*m3;
    return mass(tmp);
  }

◆ mass() [6/6]

double Vertex::mass	(	double	m1,
		double	m2,
		double	m3,
		double	m4
	)		const

4-track mass with rest-masses

                                                                      {
    if(ntracks() != 4) return -1;
    vector<double> tmp(4);
    tmp[0] = m1*m1, tmp[1] = m2*m2, tmp[2] = m3*m3, tmp[4] = m4*m4;
    return mass(tmp);
  }

◆ massCC()

double Vertex::massCC ( const bool use = false ) const

rest-mass 0 or conjugated track rest mass

                                            {
    return mass(vector<double>(ntracks(),0), use, true);
  }

◆ MassConstr()

void Vertex::MassConstr ( double m )

add a mass constraint for all tracks, $m$ = mass of mother particle

                                  {
    MassC* massc = new MassC(m);
    const RelationList::iterator last = end();
    for(RelationList::iterator it = begin(); it != last; ++it) {
      massc->push_back(*it);
    }
 
    xmass.push_back(massc);
    return;
  }

◆ ndf()

unsigned short int Vertex::ndf ( ) const

degrees of freedom of vertex fit

239 { return ((v_Mvalid == true) ? 2*ntracks()-3 + xmass.size() : 2*ntracks()-3); }

◆ nMassConstr()

unsigned int VERTEX::Vertex::nMassConstr ( ) const

number of mass constraints

◆ ntracks()

unsigned short int Vertex::ntracks ( ) const

no of tracks in vertex $n$

240{ return size(); }

◆ operator-()

const double Vertex::operator- ( const Vertex & rhs ) const

                                                        {
    return sqrt(sqr(vx()-rhs.vx()) + sqr(vy()-rhs.vy()) + sqr(vz()-rhs.vz()));
  }

◆ operator=()

const Vertex & Vertex::operator= ( const Vertex & rhs )

                                                   {
    if( this == &rhs ) { return *this; }
    
    Track::operator=(rhs);
 
    v_use_kalman   = rhs.v_use_kalman;
    v_valid        = rhs.v_valid;
    v_Mvalid       = rhs.v_Mvalid;
    v_mother       = rhs.v_mother;
    v_angdist      = rhs.v_angdist;
    v_angle        = rhs.v_angle;
    v_dist         = rhs.v_dist;
    v_bk13         = rhs.v_bk13;
    v_chi2         = rhs.v_chi2;
 
    kal_xv         = rhs.kal_xv;
    kal_xvs        = rhs.kal_xvs;
    v_CS           = rhs.v_CS;
 
    // possible if tracks & xmass contain SmartPointers to the objects
    //      tracks         = rhs.tracks;
    //      xmass          = rhs.xmass;
 
    return *this;
  }

◆ operator==()

const bool Vertex::operator== ( const Vertex & rhs ) const

                                                       {
    return RelationList::operator==(rhs);
  }

◆ pmaxfrac()

double Vertex::pmaxfrac ( ) const

$\max|p_i/p_j|$ max. momentum fraction

                                {
    double maxfrac = 0.;
 
    const RelationList::const_iterator last = end();
    for(RelationList::const_iterator tit = begin(); tit != last; ++tit) {
      double p1 = fabs(tit->track.p());
 
      RelationList::const_iterator it=tit;
      for(++it; it != last; ++it) {
    double p2 = fabs(it->track.p());
    double frac = max(p1/p2, p2/p1);
    if( frac > maxfrac ) maxfrac = frac;
      }
    }
 
    return maxfrac;
  }

◆ print()

std::ostream & Vertex::print ( std::ostream & os ) const

called by operator<<()

                                                {
    os.setf(ios::right,ios::adjustfield);
    os.setf(ios::fixed,ios::floatfield);
    if(v_mother == true) {
      os << " Mother Track: ";
      Track::print(os);
      os << endl;
    }
//       if (1) return os;
    os << " Vertex:"
       << " vx: "       << setw(6) << setprecision(3) << vx()
       << " vy: "       << setw(6) << setprecision(3) << vy()
       << " vz: "       << setw(7) << setprecision(3) << vz()
//       if (1) return os;
       << " chi2: "     << setw(6) << setprecision(2) << chi2()
       << " prob: "     << setw(6) << setprecision(2) << prob()
       << " angle: "    << setw(5) << setprecision(3) << angle() 
       << " dist: "     << setw(6) << setprecision(4) << dist()
       << " mass: "     << setw(6) << setprecision(4) << mass(v_use_kalman)
       << " no of tracks: " << ntracks() 
       << endl;
 
    if(xmass.size() > 0) {
      os << "Mass constraints: " << endl;
      const MassC_cit end = xmass.end();
      for(MassC_cit it = xmass.begin(); it != end; ++it) {
    os << "   " << **it << endl;
      }
    }
 
#ifdef VtDEBUG
    os << " Relations: ";
 
    RelationList::print(os);
#endif
 
    os << "Vertex cov. matrix: " << CS();
 
#ifdef VtDEBUG
    if(v_covn != 0) {
      cout << " bigcov: " << *v_covn << endl;
    }
#endif
 
    return os;
  }

◆ prob()

float Vertex::prob ( ) const

upper tail $\chi^2$ probability

237{ return CERN_prob(chi2(), ndf()); }

VERTEX::Vertex::ndf

unsigned short int ndf() const

degrees of freedom of vertex fit

Definition: VtVertex.C:238

◆ push_back()

void Vertex::push_back ( Track & track )

                                     {
    createRelation(track, *this);
    invalid();
  }

◆ remove()

void Vertex::remove ( Relation * r )

virtual

unbook relation

Reimplemented from VERTEX::RelationList.

                                {
    RelationList::remove(r);
    invalid();
  }

◆ remove_last()

const bool Vertex::remove_last ( )

remove last track from vertex, do inverse Kalman filter

                                 {
    return VtRemoveTrack(*back());
  }

◆ remove_worst()

void Vertex::remove_worst ( )

remove track with biggest $\chi^2$ contribution, do inverse Kalman filter

                            {
    double prchi2 = 0.;
    double maxchi2 = 0.;
    Relation* worst = 0;
 
    // a sorted list would be better here
    const iterator last = end();
    for(iterator it=begin(); it!=last; ++it) {
      const Kalman& k = it->kalman;
      double dchi2 = k.chi2() - prchi2;
      prchi2 = k.chi2();
 
      if(dchi2 > maxchi2){
    maxchi2 = dchi2;
    worst = &(*it);
      }
    }
    
    if(worst != 0) VtRemoveTrack(*worst);
    //    if(worst != 0) remove(worst);
    return;
  }

◆ rmsDistAngle()

double Vertex::rmsDistAngle ( void ) const

calc rms dist and rms angle

                                        {
 
    // and the rms distance to the vertex
    v_dist = 0.;
    double tsq, rsq, rt, t[3], r[3];
    double sumtx = 0., sumty = 0.;
    double sumtx2 = 0., sumty2 = 0.;
    const RelationList::const_iterator last = end();
    for(RelationList::const_iterator pos = begin(); pos != last; ++pos) {
      const Track& track = pos->track;
 
      t[0] = track.tx();  r[0] = track.x() - kal_xvs[0];
      t[1] = track.ty();  r[1] = track.y() - kal_xvs[1];
      t[2] = 1;           r[2] = track.z() - kal_xvs[2];
 
      sumtx += t[0]; sumtx2 += sqr(t[0]);
      sumty += t[1]; sumty2 += sqr(t[1]);
 
      tsq = sqr(t[0]) + sqr(t[1]) + sqr(t[2]);
      rsq = sqr(r[0]) + sqr(r[1]) + sqr(r[2]);
      rt  = r[0]*t[0] + r[1]*t[1] + r[2]*t[2];
 
      v_dist += rsq - sqr(rt)/tsq;
    } 
 
    unsigned int v_nrseg = size();
    v_angle = sqrt(sumtx2*v_nrseg - sqr(sumtx) + sumty2*v_nrseg -
           sqr(sumty))/static_cast<double>(v_nrseg); 
    
    v_dist /= static_cast<double>(v_nrseg-1);
    v_dist  = sqrt(v_dist);
    v_angdist = true;
    return v_dist;
  }

◆ set_invalid()

void VERTEX::Vertex::set_invalid ( )

mark vertex as invalid

◆ track_chi2()

double Vertex::track_chi2 ( int i )

get track $\chi^2$ contribution

                                   {
    double prchi2 = 0.;
    double dchi2  = 0.;
    // a sorted list would be better here
    const iterator last = end();
    int ntr = 0;
    for(iterator it=begin(); it!=last; ++it) {
      const Kalman& k = it->kalman;
      dchi2 = k.chi2() - prchi2;
      if (ntr == i) return dchi2;
      prchi2 = k.chi2();
      ntr++;
    }
    
    return dchi2;
  }

◆ track_worst()

int Vertex::track_worst ( )

get track with biggest $\chi^2$ contribution

                          {
    double prchi2 = 0.;
    double maxchi2 = 0.;
 
    // a sorted list would be better here
    const iterator last = end();
    int num = -1, ntr = 0;
    for(iterator it=begin(); it!=last; ++it) {
      const Kalman& k = it->kalman;
      double dchi2 = k.chi2() - prchi2;
      prchi2 = k.chi2();
 
      if(dchi2 > maxchi2){
    maxchi2 = dchi2;
    num = ntr;
      }
      ntr++;
    }
    
    return num;
  }

◆ use_kalman() [1/2]

bool VERTEX::Vertex::use_kalman ( ) const

◆ use_kalman() [2/2]

void VERTEX::Vertex::use_kalman ( const bool use )

use refitted track parameters or not

◆ use_momentum()

void Vertex::use_momentum ( const bool use )

set for all tracks whether momentum should be used or not

                                          {
    const RelationList::iterator last = end();
    for(RelationList::iterator tit = begin(); tit != last; ++tit) {
      tit->kalman.use_momentum(use);
    }
    return;
  }

◆ valid()

bool VERTEX::Vertex::valid ( ) const

is vertex valid?

◆ VtEstimateVertex()

const bool Vertex::VtEstimateVertex ( )

weighted 2-D min. dist.

                                      {
    //return VtEstimateVertexMath(kal_xvs[0], kal_xvs[1], kal_xvs[2]);
    return VtEstimateVertexMathTA(kal_xvs[0], kal_xvs[1], kal_xvs[2]);
  }

◆ VtEstimateVertexMath()

const bool Vertex::VtEstimateVertexMath	(	double &	x,
		double &	y,
		double &	z
	)

estimate Vertex without changing Vertex object

                                                                         {
 
    double tx_sum  = 0.;
    double x_sum   = 0.;
    double xw_sum  = 0.;
    double xtx_sum = 0.;
    double tx2_sum = 0.;
    double ty_sum  = 0.;
    double y_sum   = 0.;
    double yw_sum  = 0.;
    double yty_sum = 0.;
    double ty2_sum = 0.;
 
    // fill cumulants
    const RelationList::iterator last = end();
    for(RelationList::iterator tit = begin(); tit != last; ++tit) {
      const Track& track = tit->track;
 
      double x        = track.x();
      double y        = track.y();
      double tx       = track.tx();
      double ty       = track.ty();
      double xweight  = 1./track.V().x();
      double yweight  = 1./track.V().y();
      double xweight2 = xweight*xweight;
      double yweight2 = yweight*yweight;
 
      tx_sum  += tx * xweight;
      x_sum   += x  * xweight;
      xw_sum  += xweight;
      xtx_sum += x  * tx * xweight2;
      tx2_sum += tx * tx * xweight2;
 
      ty_sum  += ty * yweight;
      y_sum   += y  * yweight;
      yw_sum  += yweight;
      yty_sum += y  * ty * yweight2;
      ty2_sum += ty * ty * yweight2;
 
    } // for track
 
    double det = -tx2_sum - ty2_sum + tx_sum*tx_sum/xw_sum + ty_sum*ty_sum/yw_sum;
 
//#ifdef VtDEBUG    
    cout << " txsum: " << tx_sum << " xsum: " << x_sum << " xwsum: " << xw_sum
     << " xtxsum: " << xtx_sum << " tx2sum: " << tx2_sum << endl;
    cout << " tysum: " << ty_sum << " ysum: " << y_sum << " ywsum: " << yw_sum
     << " ytysum: " << yty_sum << " ty2sum: " << ty2_sum << endl;
    cout << " det: " << det << endl;
//#endif
 
    if(det == 0.) {
//#ifdef VtDEBUG
      cout << "VtVertex::VtEstimateVertex error: determinant == 0.!" << endl;
//#endif
      return false;
    }
 
//      kal_xvs[2] = ( xtx_sum + yty_sum - tx_sum*x_sum/xw_sum - ty_sum*y_sum/yw_sum ) / det;
//      kal_xvs[0] = ( x_sum + tx_sum * vz() ) / xw_sum;
//      kal_xvs[1] = ( y_sum + ty_sum * vz() ) / yw_sum;
 
    z = ( xtx_sum + yty_sum - tx_sum*x_sum/xw_sum - ty_sum*y_sum/yw_sum ) / det;
    x = ( x_sum + tx_sum * z ) / xw_sum;
    y = ( y_sum + ty_sum * z ) / yw_sum;
 
    v_valid = true;
    return true;
  }

◆ VtEstimateVertexMathTA()

const bool Vertex::VtEstimateVertexMathTA	(	double &	x,
		double &	y,
		double &	z
	)

estimate Vertex without changing Vertex object

                                                                              {
  // Tomoko Ariga version inserted by VT 1/09/2011
  // calc vertex point with given segments. just topological calculation
 
    double Det,Ax,Ay,Az,Bx,By,Bz,Cx,Cy,Cz,Dx,Dy,Dz;
    Ax=0.;Ay=0.;Az=0.;Bx=0.;By=0.;Bz=0.;Cx=0.;Cy=0.;Cz=0.;Dx=0.;Dy=0.;Dz=0.;
 
    const RelationList::iterator last = end();
    for(RelationList::iterator tit = begin(); tit != last; ++tit) {
        const Track& track = tit->track;
 
        double x        = track.x();
        double y        = track.y();
        double ax       = track.tx();
        double ay       = track.ty();
        double az = 1.0;
    double z  = track.z();
    double a = ax*ax+ay*ay+az*az;
    double c = -ax*x-ay*y-az*z;
    double b = (ax*ax+ay*ay);
//      double w = 1.0/a/a; // weight for small angle tracks.
    double w = 1.0; // no weight
 
    Ax+=2.0*w/b*( a*(ay*ay+az*az) );
    Bx+=2.0*w/b*( -a*ax*ay );
    Cx+=2.0*w/b*( -a*ax*az );
    Dx+=2.0*w/b*( -(a*x+c*ax)*(ax*ax-a)-(a*y+c*ay)*ax*ay-(a*z+c*az)*az*ax );
 
    Ay+=2.0*w/b*( -a*ay*ax );
    By+=2.0*w/b*( a*(az*az+ax*ax) );
    Cy+=2.0*w/b*( -a*ay*az );
    Dy+=2.0*w/b*( -(a*y+c*ay)*(ay*ay-a)-(a*z+c*az)*ay*az-(a*x+c*ax)*ax*ay );
 
    Az+=2.0*w/b*( -a*az*ax );
    Bz+=2.0*w/b*( -a*az*ay );
    Cz+=2.0*w/b*( a*b );
    Dz+=2.0*w/b*( -(a*z+c*az)*(az*az-a)-(a*x+c*ax)*az*ax-(a*y+c*ay)*ay*az );
 
     }
 
    Det=fabs( Ax*(By*Cz-Cy*Bz)-Bx*(Ay*Cz-Cy*Az)+Cx*(Ay*Bz-By*Az) );
    if( Det < 0.00000001 ) {
      cout << "VtVertex::VtEstimateVertexMathTA error: determinant is too small! " << Det << endl;
      return false;
    }
    xx=( (By*Cz-Cy*Bz)*Dx-(Bx*Cz-Cx*Bz)*Dy+(Bx*Cy-Cx*By)*Dz)/Det;
    yy=(-(Ay*Cz-Cy*Az)*Dx+(Ax*Cz-Cx*Az)*Dy-(Ax*Cy-Cx*Ay)*Dz)/Det;
    zz=( (Ay*Bz-By*Az)*Dx-(Ax*Bz-Bx*Az)*Dy+(Ax*By-Bx*Ay)*Dz)/Det;
    
    v_valid = true;
    return true;
  }

◆ VtFilter()

const bool Vertex::VtFilter ( )

                              {
    
    // initialize by vertex position & cov. matrix
    //    const VtVector*    prkal_xv = &kal_xv0;
    const VtVector*    prkal_xv = &kal_xvs;
    const VtSymMatrix* prCINV   = &v_CINV;
    double             prChi2   = 0.;
    v_bk13 = vz();
 
    // cout << "bk13: " << v_bk13 << endl;
 
    const RelationList::iterator last = end();
    for(RelationList::iterator tit = begin(); tit != last; ++tit) {
      Kalman& kalman = tit->kalman;
 
      //      cout << "----> FILTER TRACK " << track;
#ifdef VtDEBUG
      cout << endl;
#endif
 
      // compute contributions of track k
      if(kalman.filter(v_bk13, *prCINV, *prkal_xv) == false) {
    return false; // forget this vertex
      }
 
      kalman.filter_chi2(v_bk13, prChi2, *prCINV, *prkal_xv);
    
      //      cout << "Track-Kalman chi2: " << kalman.chi2() << endl;
 
      // store pointers to previous (k-1) values
      prkal_xv = &kalman.xv();
      prCINV   = &kalman.CINV();
      prChi2   = kalman.chi2();
    } // for tracks
 
    return true;
  }

◆ VtFit()

const bool Vertex::VtFit ( )

Kalman filter + smoother.

                           {
    if(VtFilter() == false) return false;
    else {
      VtSmoothX();
      VtSmoothQ();
    }
    return true;
  }

◆ VtInverseFilter()

const bool Vertex::VtInverseFilter ( ) const

                                           {
    
    // get last track (in Vt: x_n)
    const Kalman& trackn = back()->kalman;
 
    const VtSymMatrix* prCINV = &v_CINV;
 
    const RelationList::const_iterator last = end();
    for(RelationList::const_iterator tit = begin(); tit != last; ++tit) {
      Kalman& kalman = tit->kalman;
      
      if(kalman.inverse_filter(v_bk13, trackn.CINV(), *prCINV,
                   trackn.xv()) == false) {
    return false;
      }
      prCINV   = &kalman.CINV();
    }
 
    return true;
  }

◆ VtMass()

const bool Vertex::VtMass ( )

Kalman filter with mass constraints.

                            {
    v_Mvalid = false;
 
    // should be assured:
    use_momentum(true);
    if(valid() == false && findVertexVt() == false) return false;
 
    // check kinematic boundaries
    const MassC_cit xend = xmass.end();
    for(MassC_cit it = xmass.begin(); it != xend; ++it) {
      if((*it)->check_kinematic() == false)
    return false;
    }
 
    // fill general state vectors (in Vt: alp(1,...), alpc(1,...))
    VtVector alpha (bigdim());
    VtVector alphac(bigdim());
 
    // initialize general state vector alphac (in Vt: alpc)
    alphac.place_at(xvs(),0);
 
    unsigned int i=3;
    const iterator last = end();
    for(iterator it = begin(); it != last; ++it,i+=3) {
      it->kalman.alpc_init();
      alphac.place_at(it->kalman.alpc(),i);
    }
 
#ifdef VtDEBUG
    cout << " alphac0: " << alphac << endl;
#endif
 
    double chi2old = -100.;
    double dchi2 = 0.;
    double mc_chi2 = 0.;
    const unsigned int maxiter = 100 * nMassConstr();
 
    // declare some matrices which are needed outside do-loop
    VtSymMatrix gaVga(nMassConstr());
    VtMatrix    Vgam(bigdim(),nMassConstr());
 
    // iteration
    unsigned int iter = 0;
    do {
 
      if(iter > 0) {
    // redo filter and smoother
    kal_xvs.copy(alphac,0);
    unsigned int i=3;
    for(iterator it = begin(); it != last; ++it,i+=3) {
      Kalman& kalm = it->kalman;
      kalm.qvs_nc().copy(alphac,i);
      kalm.alpc_init();
    }
    if((v_valid = VtFit()) == false)
      return false;
      }
 
      // fill general state vector alpha (in Vt: alp)
      alpha.place_at(xvs(),0);
 
      i=3;
      for(iterator it = begin(); it != last; ++it,i+=3) {
    Kalman& kalm = it->kalman;
 
    kalm.alp_init();
    kalm.calc_ealpc(); // n in Vt (eq. 40)
    alpha.place_at(kalm.alp(),i);
      }
 
#ifdef VtDEBUG
      cout << " alpha: " << alpha << endl;
#endif
 
      // get general covariance matrix v_covn (in Vt: Valp)
      const VtSymMatrix& Valp = bigcov();
 
#ifdef VtDEBUG
      cout << "Valp: " << Valp << endl;
#endif
 
      // build up matrix gam (in Vt: Gamma^T, eq. 34, 44, 45, 46)
      VtMatrix gam(bigdim(),nMassConstr());
 
      // vector \vect{g}(\vect{\alpha}^(0))
      VtVector galp(nMassConstr());
 
      unsigned int j=0;
      for(MassC_cit mit = xmass.begin(); mit != xend; ++mit,++j) {
    MassC* mC = *mit;
    mC->calc_Sl();
    galp[j] = mC->galp();
 
    i=3;
    for(iterator it = begin(); it != last; ++it,i+=3) {
      Relation& r = *it;
      if(find(r.massc.begin(),r.massc.end(),mC) != r.massc.end()) {
        gam.place_at(mC->gam(r),i,j);
      }
    }
      } // end of for(MassC_cit)
 
#ifdef VtDEBUG
      cout << " galp: " << galp << endl;
      cout << " gam: " << gam << endl;
#endif
 
      // calc (Gamma * V_alpha * Gamma^T)^-1 (in Vt: eq. 36,37,39)
      gaVga = Valp.product(gam);
      if(gaVga.VtDsinv() == false)
    return false;
      
#ifdef VtDEBUG
      cout << " gaVga^-1: " << gaVga << endl;
#endif
 
      // calc g(alpha) + Gamma * (alpha - alpha_c)
      VtVector gnum = galp + gam.T() * (alpha - alphac);
 
      // calc \vect{\mu} (in Vt: xmu, eq. 36)
      VtVector xmu = gaVga * gnum;
 
#ifdef VtDEBUG
      cout << " xmu: " << xmu << endl;
#endif
 
      // calc V_alpha*Gamma^T
      VtMatrix Vgam = Valp * gam;
 
      // calc alpha_c (in Vt: eq. 37)
      alphac = alpha - Vgam * xmu;
      
#ifdef VtDEBUG
      cout << " alphac: " << alphac << endl;
#endif
 
      // calc chi^2 (in Vt: chi2m, eq. 39)
      // this is an incremental chi2!
      // (chi2 due to mass constraint, to be added to
      // chi2 of unconstraint vertex fit!)
      mc_chi2 = gaVga.product(gnum);
 
      double chi  = mc_chi2 + chi2n();
      dchi2 = fabs(chi - chi2old);
      chi2old = chi;
 
#ifdef VtDEBUG
      cout << "chi2m: " << setprecision(7) << mc_chi2 
       << " chi: " << chi2n()
       << " fabs(chi2m + chi2() - chi2old: " << dchi2
       << endl;
#endif
    
    } while(++iter < maxiter && dchi2 > 0.00001);
 
#ifdef VtDEBUG
    cout << " iterations: " << iter << endl;
#endif
    if(iter >= maxiter) return false;
 
    VtSymMatrix& cov = *v_covn;
 
    // v_covn = Valp => v_covn = Valp - gaVga.product(Vgam.T())
    // calc v_covn (in Vt: cov(alpha_c), eq. 38)
    cov -= gaVga.product(Vgam.T());
 
#ifdef VtDEBUG
    cout << " v_covn: " << cov << endl;
#endif
 
    // get results
    v_CS.copy(cov);
    // add mass constraint chi2 term to vertex fit term
    v_chi2 += mc_chi2;
 
    kal_xvs.copy(alphac,0);
    i=3;
    for(iterator it = begin(); it != last; ++it,i+=3) {
      it->kalman.qvs_nc().copy(alphac,i);
    }
    
    return v_Mvalid = true;
  } // end of VtMass()

◆ VtRemoveTrack()

const bool Vertex::VtRemoveTrack ( Relation & track )

remove track using inverse Kalman filter

                                                   {
    if(ntracks() <= 2 || valid() == false) return false; // exclude weird cases
 
    if(VtInverseFilter() == false) return false;
 
    // get last track (in Vt: x_n)
    const Kalman& trackn = back()->kalman;
 
    // compute chi2_n -> chi2_n - chi2_k,S (eq. 22)
    double dchi2 = trackn.chi2() - trackk.kalman.chi2s();
    
    // new smoothed vertex position
    kal_xvs = trackk.kalman.xnk();
 
    const RelationList::iterator last = end();
    for(RelationList::iterator tit = begin(); tit != last; ++tit) {
      // calculate qvs
      tit->kalman.calc_qvs(kal_xvs);
    }
 
    remove(&trackk);
 
    v_chi2 = back()->kalman.set_chi2(dchi2);
 
    // refit the vertex
    return v_valid = VtFit();
  }

◆ VtSmoothQ()

void Vertex::VtSmoothQ ( )

smooth momenta

                         {
   
#ifdef VtDEBUG
    cout << "---> VtSmoothQ called! <---" << endl;
#endif
 
    const Kalman& trackn = back()->kalman;
    const RelationList::iterator last = end();
    for(RelationList::iterator tit = begin(); tit != last; ++tit) {
      tit->kalman.smooth(v_bk13, kal_xvs, trackn.C());
    }
    
    return;
  }

◆ VtSmoothX()

void Vertex::VtSmoothX ( )

smooth vertex position

                         {
    // get last track (in Vt: x_n)
    const Kalman& kaln = back()->kalman;
    
#ifdef VtDEBUG
    cout << "--> VtSmooth called!" << endl;
#endif
 
    // update vertex position saved in xvs!!
    kal_xvs = kaln.xv(); 
 
#ifdef VtDEBUG
    cout << "kal_xvs: "    << kal_xvs  << endl;
    cout << "track->C(): " << kaln.C() << endl;
#endif
 
    // save smoothed covariance matrix of vertex
    v_CS = kaln.C();
 
    // save chi2
    v_chi2 = kaln.chi2();
 
#ifdef VtDEBUG
    cout << "v_CS: " << v_CS << endl;
#endif
 
    return;
  }

◆ vtx_cov_x()

float Vertex::vtx_cov_x ( ) const

$\sigma_x^2$ of vertex

254{ return v_CS(0,0); }

◆ vtx_cov_y()

float Vertex::vtx_cov_y ( ) const

$\sigma_y^2$ of vertex

255{ return v_CS(1,1); }

◆ vtx_cov_z()

float Vertex::vtx_cov_z ( ) const

$\sigma_z^2$ of vertex

256{ return v_CS(2,2); }

◆ vx()

float Vertex::vx ( ) const

$x$ of vertex

233{ return kal_xvs[0]; }

◆ vxerr()

double VERTEX::Vertex::vxerr ( ) const

$\sqrt{\sigma_{vx}^2}$ vertex $x$-error

◆ vy()

float Vertex::vy ( ) const

$y$ of vertex

234{ return kal_xvs[1]; }

◆ vyerr()

double VERTEX::Vertex::vyerr ( ) const

$\sqrt{\sigma_{vy}^2}$ vertex $y$-error

◆ vz()

float Vertex::vz ( ) const

$z$ of vertex

235{ return kal_xvs[2]; }

◆ vzerr()

double VERTEX::Vertex::vzerr ( ) const

$\sqrt{\sigma_{vz}^2}$ vertex $z$-error

◆ xv()

const MATRIX::VtVector & VERTEX::Vertex::xv ( ) const

◆ xvs()

const MATRIX::VtVector & VERTEX::Vertex::xvs ( ) const

vector of vertex-positions

Member Data Documentation

◆ kal_xv

MATRIX::VtVector VERTEX::Vertex::kal_xv

private

◆ kal_xvs

MATRIX::VtVector VERTEX::Vertex::kal_xvs

private

◆ tracks

Track_v VERTEX::Vertex::tracks

private

◆ v_angdist

bool VERTEX::Vertex::v_angdist

mutableprivate

◆ v_angle

double VERTEX::Vertex::v_angle

mutableprivate

◆ v_bk13

double VERTEX::Vertex::v_bk13

private

◆ v_chi2

double VERTEX::Vertex::v_chi2

private

◆ v_CINV

MATRIX::VtSymMatrix VERTEX::Vertex::v_CINV

private

◆ v_covn

MATRIX::VtSymMatrix* VERTEX::Vertex::v_covn

private

◆ v_CS

MATRIX::VtSymMatrix VERTEX::Vertex::v_CS

private

◆ v_dist

double VERTEX::Vertex::v_dist

mutableprivate

◆ v_mother

bool VERTEX::Vertex::v_mother

private

◆ v_Mvalid

bool VERTEX::Vertex::v_Mvalid

private

◆ v_use_kalman

bool VERTEX::Vertex::v_use_kalman

private

◆ v_valid

bool VERTEX::Vertex::v_valid

private

◆ xmass

MassC_v VERTEX::Vertex::xmass

private

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

/home/antonio/fedra_doxygen/src/libVt++/vt++/include/VtVertex.hh
/home/antonio/fedra_doxygen/src/libVt++/vt++/VtVertex.C

Public Member Functions

— Operators —

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ Vertex() [1/4]

◆ Vertex() [2/4]

◆ Vertex() [3/4]

◆ Vertex() [4/4]

◆ ~Vertex()

Member Function Documentation

◆ add_track()

◆ addMassConstr()

◆ angle()

◆ bigcov()

◆ bigdim()

◆ calc_mother()

◆ calc_mother_cov()

◆ chi2()

◆ chi2l()

◆ chi2n()

◆ ClassDef()

◆ clear()

◆ clearMassConstr()

◆ covn()

◆ CS()

◆ dist()

◆ distance() [1/3]

◆ distance() [2/3]

◆ distance() [3/3]

◆ erase()

◆ findVertex2D()

◆ findVertex3D()

◆ findVertexVt()

◆ get_track()

◆ invalid()

◆ mass() [1/6]

◆ mass() [2/6]

◆ mass() [3/6]

◆ mass() [4/6]

◆ mass() [5/6]

◆ mass() [6/6]

◆ massCC()

◆ MassConstr()

◆ ndf()

◆ nMassConstr()

◆ ntracks()

◆ operator-()

◆ operator=()

◆ operator==()

◆ pmaxfrac()

◆ print()

◆ prob()

◆ push_back()

◆ remove()

◆ remove_last()

◆ remove_worst()

◆ rmsDistAngle()

◆ set_invalid()

◆ track_chi2()

◆ track_worst()

◆ use_kalman() [1/2]

◆ use_kalman() [2/2]

◆ use_momentum()

◆ valid()

◆ VtEstimateVertex()

◆ VtEstimateVertexMath()

◆ VtEstimateVertexMathTA()

◆ VtFilter()

◆ VtFit()

◆ VtInverseFilter()

◆ VtMass()

◆ VtRemoveTrack()

◆ VtSmoothQ()

◆ VtSmoothX()

◆ vtx_cov_x()

◆ vtx_cov_y()

◆ vtx_cov_z()

◆ vx()

◆ vxerr()