Trk::KalmanVertexOnJetAxisUpdator Class Reference

ATLAS Reconstruction Software - (C) 2005 - 2007. More...

#include <KalmanVertexOnJetAxisUpdator.h>

Inheritance diagram for Trk::KalmanVertexOnJetAxisUpdator:

Collaboration diagram for Trk::KalmanVertexOnJetAxisUpdator:

Public Member Functions
virtual StatusCode	initialize () override
	KalmanVertexOnJetAxisUpdator (const std::string &t, const std::string &n, const IInterface *p)
	Constructor.
void	add (VxTrackAtVertex *trackToAdd, const VxVertexOnJetAxis *vertexToUpdate, VxJetCandidate *candidateToUpdate) const
	Method updating the fit along the flight axis (VxJetCandidate) with a track (VxTrackOnJetAxis) which is added to the vertex along the jet axis provided (VxVertexOnJetAxis).
void	addWithFastUpdate (VxTrackAtVertex *trackToAdd, const VxVertexOnJetAxis *vertexToUpdate, VxJetCandidate *candidateToUpdate) const
void	remove (VxTrackAtVertex *trackToRemove, const VxVertexOnJetAxis *vertexToUpdate, VxJetCandidate *candidateToUpdate) const
	Method removing already added track from the fit along the flight axis, where you specify from what vertex along the jet axis the track should be removed.
void	removeWithFastUpdate (VxTrackAtVertex *trackToRemove, const VxVertexOnJetAxis *vertexToUpdate, VxJetCandidate *candidateToUpdate) const
void	updateChi2NdfInfo (VxTrackAtVertex *trackToUpdate, const VxVertexOnJetAxis *vertexToUpdate, VxJetCandidate *vertexCandidate) const
	Obsolete (will be removed - not in use, however migrated to eigen before finding out -WL) Updates only the chi2 + ndf information, in case the weight matrix is not inverted at every step, so that there is no info about the fitted position at intermediate steps of the fit.
void	updateVertexChi2 (VxJetCandidate *vertexCandidate) const
	Obsolete (will be removed - not in use, however moved to eigen before finding out -WL).
Trk::RecVertexPositions	positionUpdate (VxJetCandidate &vertexCandidate, const LinearizedTrack *trk, double trackWeight, int sign, int numVertex, bool isPrimary, bool doFastUpdate) const
	Method to do the Kalman Update of the positions of the fit along the flight axis with a single new track (which belongs to the VxVertexOnJetAxis number numVertex).
std::pair< Amg::Vector3D, Eigen::Matrix3Xd >	createTransformJacobian (const Amg::VectorX &initialjetdir, int numVertex, bool isPrimary, bool truncate) const
	Method to create the Jacobian needed to go to the space of single vertices to the space of the variables representing the full fit along the flight axis (primary vertex, phi and theta of jet direction, distance to primary vertex of the remaining vertices along the jetaxis).
void	smartInvert (Amg::MatrixX &new_vrt_cov) const
ServiceHandle< StoreGateSvc > &	evtStore ()
	The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
const ServiceHandle< StoreGateSvc > &	detStore () const
	The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
virtual StatusCode	sysInitialize () override
	Perform system initialization for an algorithm.
virtual StatusCode	sysStart () override
	Handle START transition.
virtual std::vector< Gaudi::DataHandle * >	inputHandles () const override
	Return this algorithm's input handles.
virtual std::vector< Gaudi::DataHandle * >	outputHandles () const override
	Return this algorithm's output handles.
Gaudi::Details::PropertyBase &	declareProperty (Gaudi::Property< T, V, H > &t)
void	updateVHKA (Gaudi::Details::PropertyBase &)
MsgStream &	msg () const
bool	msgLvl (const MSG::Level lvl) const

Static Public Member Functions
static const InterfaceID &	interfaceID ()

Protected Member Functions
void	renounceArray (SG::VarHandleKeyArray &handlesArray)
	remove all handles from I/O resolution
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void >	renounce (T &h)
void	extraDeps_update_handler (Gaudi::Details::PropertyBase &ExtraDeps)
	Add StoreName to extra input/output deps as needed.

Private Types
typedef ServiceHandle< StoreGateSvc >	StoreGateSvc_t

Private Member Functions
void	update (VxTrackAtVertex *trackToUpdate, const VxVertexOnJetAxis *vertexToUpdate, VxJetCandidate *candidateToUpdate, int sign, bool doFastUpdate) const
	Internal method to do the Kalman Update.
double	calculateTrackChi2 (const VxJetCandidate &vertexCandidate, const LinearizedTrack *trk, double trackWeight, int numVertex, bool isPrimary=false) const
	Internal method where the chi2/ndf update is actually done.
Gaudi::Details::PropertyBase &	declareGaudiProperty (Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
	specialization for handling Gaudi::Property<SG::VarHandleKey>

Private Attributes
double	m_initialMomentumError
StoreGateSvc_t	m_evtStore
	Pointer to StoreGate (event store by default)
StoreGateSvc_t	m_detStore
	Pointer to StoreGate (detector store by default)
std::vector< SG::VarHandleKeyArray * >	m_vhka
bool	m_varHandleArraysDeclared

Detailed Description

ATLAS Reconstruction Software - (C) 2005 - 2007.

Author

Giacinto Piacquadio, (giaci.nosp@m.nto..nosp@m.piacq.nosp@m.uadi.nosp@m.o@phy.nosp@m.sik..nosp@m.uni-f.nosp@m.reib.nosp@m.urg.d.nosp@m.e)

Christian Weiser, (chris.nosp@m.tian.nosp@m..weis.nosp@m.er@p.nosp@m.hysik.nosp@m..uni.nosp@m.-frei.nosp@m.burg.nosp@m..de)

A concrete implementation of the Vertex Updator using the Kalman filter algorithm.

This is a special implementation (part of the JetFitter algorithm) which is able to fit a vertex correcting the jet axis parameters to find the best overall agreement.

Basic approach is based on R. Fruhwirth et al. Comp. Phys. Comm 96(1996) 189

The first version implemented in TrkVertexFitterUtils, for general vertexing purposes, was created by Kirill Prokofiev and G.P.

First version: December 2006 Updated version for athena: Februar 2007

Alberts-Ludwig-Universita"t Freiburg

Definition at line 54 of file KalmanVertexOnJetAxisUpdator.h.

Member Typedef Documentation

StoreGateSvc_t

typedef ServiceHandle<StoreGateSvc> AthCommonDataStore< AthCommonMsg< AlgTool > >::StoreGateSvc_t

privateinherited

Definition at line 388 of file AthCommonDataStore.h.

Constructor & Destructor Documentation

KalmanVertexOnJetAxisUpdator()

Trk::KalmanVertexOnJetAxisUpdator::KalmanVertexOnJetAxisUpdator	(	const std::string &	t,
		const std::string &	n,
		const IInterface *	p )

Constructor.

Definition at line 25 of file KalmanVertexOnJetAxisUpdator.cxx.

                                                                                                                        :
    AthAlgTool(t,n,p),
    m_initialMomentumError(1000.)
  {
    declareInterface<KalmanVertexOnJetAxisUpdator>(this);
    declareProperty("initialMomentumError",m_initialMomentumError);
  }

Member Function Documentation

add()

void Trk::KalmanVertexOnJetAxisUpdator::add	(	VxTrackAtVertex *	trackToAdd,
		const VxVertexOnJetAxis *	vertexToUpdate,
		VxJetCandidate *	candidateToUpdate ) const

Method updating the fit along the flight axis (VxJetCandidate) with a track (VxTrackOnJetAxis) which is added to the vertex along the jet axis provided (VxVertexOnJetAxis).

To update the fit with a new vertex (VxVertexOnJetAxis) you have to iteratively add all its tracks to the fit using this method. In the RecVertexPositions datamember of the VxJetCandidate the result will appear as the value of the distance of this fitted vertex from the primary vertex. All the other vertices and the jet direction will of course be also dynamically adjusted while adding tracks to the fit.

Definition at line 45 of file KalmanVertexOnJetAxisUpdator.cxx.

                                                                                  {
    update(trackToAdd,vertexToUpdate,candidateToUpdate,1,false);
  }

addWithFastUpdate()

void Trk::KalmanVertexOnJetAxisUpdator::addWithFastUpdate	(	VxTrackAtVertex *	trackToAdd,
		const VxVertexOnJetAxis *	vertexToUpdate,
		VxJetCandidate *	candidateToUpdate ) const

Definition at line 52 of file KalmanVertexOnJetAxisUpdator.cxx.

                                                                                                {
    update(trackToAdd,vertexToUpdate,candidateToUpdate,1,true);
  }

calculateTrackChi2()

double Trk::KalmanVertexOnJetAxisUpdator::calculateTrackChi2 ( const VxJetCandidate & vertexCandidate, const LinearizedTrack * trk, double trackWeight, int numVertex, bool isPrimary = false ) const

private

Internal method where the chi2/ndf update is actually done.

Definition at line 429 of file KalmanVertexOnJetAxisUpdator.cxx.

                                                                   {
    
    
    //linearized track information 
    const VertexPositions & linearizationPositions = vertexCandidate.getLinearizationVertexPositions();
    const Amg::VectorX & initialjetdir = linearizationPositions.position();
 
    //calculate 3-dim position on jet axis from linearizationPositionVector and calculate Jacobian of transformation
    std::pair<Amg::Vector3D,Eigen::Matrix3Xd> PosOnJetAxisAndTransformMatrix = 
      createTransformJacobian(initialjetdir, numVertex, isPrimary, true);
 
    //only position jacobian changed from A ->oldA
    const AmgMatrix(5,3)& oldA = trk->positionJacobian();
 
    //now create the new jacobian which you should use
    Eigen::Matrix<double,5,Eigen::Dynamic> A=oldA*PosOnJetAxisAndTransformMatrix.second;
   
    const AmgMatrix(5,3)& B = trk->momentumJacobian();
    const AmgVector(5) & trackParameters = trk->expectedParametersAtPCA();
 
    AmgVector(5) constantTerm = trk->constantTerm() + oldA*PosOnJetAxisAndTransformMatrix.first
      -A*initialjetdir;
 
    const AmgSymMatrix(5) & trackParametersWeight  = trk->expectedWeightAtPCA();
    
    const RecVertexPositions & myPosition = vertexCandidate.getRecVertexPositions();
    const Amg::VectorX & new_vrt_position = myPosition.position();
 
    AmgSymMatrix(3) S = B.transpose()*trackParametersWeight*B;
    if (S.determinant() == 0.0) {
      ATH_MSG_WARNING ("The matrix S is not invertible, return chi2 0");
      return -0.;
    } 
    S = S.inverse().eval();
 
    //refitted track momentum  
    Amg::Vector3D newTrackMomentum = S*B.transpose()*trackParametersWeight*(trackParameters - constantTerm - A*new_vrt_position);
 
    //refitted track parameters
    AmgVector(5) refTrackParameters = constantTerm + A * new_vrt_position + B * newTrackMomentum;
 
    //parameters difference 
    AmgVector(5) paramDifference = trackParameters - refTrackParameters;
   
    return (paramDifference.transpose()*trackParametersWeight*paramDifference)(0,0)*trackWeight;
   
  }

createTransformJacobian()

std::pair< Amg::Vector3D, Eigen::Matrix3Xd > Trk::KalmanVertexOnJetAxisUpdator::createTransformJacobian	(	const Amg::VectorX &	initialjetdir,
		int	numVertex,
		bool	isPrimary,
		bool	truncate ) const

Method to create the Jacobian needed to go to the space of single vertices to the space of the variables representing the full fit along the flight axis (primary vertex, phi and theta of jet direction, distance to primary vertex of the remaining vertices along the jetaxis).

As a result the vector of position where the linearization was done and the jacobian are provided.

If the parameter truncate is true, the method also cuts the size of the matrix down to the only needed elements (if you work on vertex number 5, all the ones which come after are cut out from the jacobian to make the fit faster).

Definition at line 515 of file KalmanVertexOnJetAxisUpdator.cxx.

                                                                                         {
 
    //now modify the position jacobian in the way you need!
    const unsigned int numrow_toupdate = isPrimary ? 4: numRow(numVertex);
    const unsigned int matrixdim = truncateDimensions ?
                                   numrow_toupdate+1  :
                                   initialjetdir.rows();
 
    //store values of RecVertexOnJetAxis
    double xv = initialjetdir[Trk::jet_xv];
    double yv = initialjetdir[Trk::jet_yv];
    double zv = initialjetdir[Trk::jet_zv];
    double phi = initialjetdir[Trk::jet_phi];
    double theta = initialjetdir[Trk::jet_theta];
    double dist = 0;
    if (!isPrimary) {
      dist = initialjetdir[numrow_toupdate];
    }
    
    ATH_MSG_VERBOSE ("actual initialjetdir values : xv " << xv << " yv " << yv << " zv " << zv << " phi " << phi << " theta " << theta << " dist " << dist);
    ATH_MSG_VERBOSE ("numrow_toupdate " << numrow_toupdate << " when primary it's -5");
    
    //now create matrix transformation to go from vertex position jacobian to vertex in jet jacobian
    //(essentially Jacobian of spherical coordinate system transformation)
    Eigen::Matrix3Xd transform(3,matrixdim);
    transform.setZero();
    transform(0,Trk::jet_xv) = 1;
    transform(1,Trk::jet_yv) = 1;
    transform(2,Trk::jet_zv) = 1;
    if (!isPrimary) {
      transform(0,Trk::jet_phi)   = -dist*sin(theta)*sin(phi);
      transform(0,Trk::jet_theta) =  dist*cos(theta)*cos(phi);
      transform(0,numrow_toupdate)=  sin(theta)*cos(phi);
      transform(1,Trk::jet_phi)   =  dist*sin(theta)*cos(phi);
      transform(1,Trk::jet_theta) =  dist*cos(theta)*sin(phi);
      transform(1,numrow_toupdate)=  sin(theta)*sin(phi);
      transform(2,Trk::jet_theta) = -dist*sin(theta);
      transform(2,numrow_toupdate)=  cos(theta);
    }
    
    ATH_MSG_DEBUG ("The transform matrix xyzphitheta is: " << transform);
    
    Amg::Vector3D posOnJetAxis;
    posOnJetAxis(0) = xv+dist*cos(phi)*sin(theta);
    posOnJetAxis(1) = yv+dist*sin(phi)*sin(theta);
    posOnJetAxis(2) = zv+dist*cos(theta);
 
    ATH_MSG_DEBUG("the lin position on jet axis x " << xv+dist*cos(phi)*sin(theta) << " y " << yv+dist*sin(phi)*sin(theta) << " z " << zv+dist*cos(theta));
 
    return std::pair<Amg::Vector3D,Eigen::Matrix3Xd>(posOnJetAxis,transform);
  }

declareGaudiProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty ( Gaudi::Property< T, V, H > & hndl, const SG::VarHandleKeyType &  )

inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleKey>

Definition at line 156 of file AthCommonDataStore.h.

  {
    return *AthCommonDataStore<PBASE>::declareProperty(hndl.name(),
                                                       hndl.value(), 
                                                       hndl.documentation());
 
  }

declareProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( Gaudi::Property< T, V, H > & t )

inlineinherited

Definition at line 145 of file AthCommonDataStore.h.

                                                                     {
    typedef typename SG::HandleClassifier<T>::type htype;
    return AthCommonDataStore<PBASE>::declareGaudiProperty(t, htype());
  }

detStore()

const ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< AlgTool > >::detStore ( ) const

inlineinherited

The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 95 of file AthCommonDataStore.h.

{ return m_detStore; }

evtStore()

ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< AlgTool > >::evtStore ( )

inlineinherited

The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 85 of file AthCommonDataStore.h.

{ return m_evtStore; }

extraDeps_update_handler()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::extraDeps_update_handler ( Gaudi::Details::PropertyBase & ExtraDeps )

protectedinherited

Add StoreName to extra input/output deps as needed.

use the logic of the VarHandleKey to parse the DataObjID keys supplied via the ExtraInputs and ExtraOuputs Properties to add the StoreName if it's not explicitly given

initialize()

StatusCode Trk::KalmanVertexOnJetAxisUpdator::initialize ( )

overridevirtual

Definition at line 34 of file KalmanVertexOnJetAxisUpdator.cxx.

  {
    StatusCode sc = AthAlgTool::initialize();
    if(sc.isFailure()){
      ATH_MSG_ERROR (" Unable to initialize the AlgTool");
      return StatusCode::FAILURE;
    }
    return StatusCode::SUCCESS;
  }

inputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< AlgTool > >::inputHandles ( ) const

overridevirtualinherited

Return this algorithm's input handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

interfaceID()

const InterfaceID & Trk::KalmanVertexOnJetAxisUpdator::interfaceID ( )

inlinestatic

Definition at line 58 of file KalmanVertexOnJetAxisUpdator.h.

     {
       return IID_KalmanVertexOnJetAxisUpdator;
     };

msg()

MsgStream & AthCommonMsg< AlgTool >::msg ( ) const

inlineinherited

Definition at line 24 of file AthCommonMsg.h.

                                {
    return this->msgStream();
  }

msgLvl()

bool AthCommonMsg< AlgTool >::msgLvl ( const MSG::Level lvl ) const

inlineinherited

Definition at line 30 of file AthCommonMsg.h.

                                               {
    return this->msgLevel(lvl);
  }

outputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< AlgTool > >::outputHandles ( ) const

overridevirtualinherited

Return this algorithm's output handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

positionUpdate()

Trk::RecVertexPositions Trk::KalmanVertexOnJetAxisUpdator::positionUpdate	(	VxJetCandidate &	vertexCandidate,
		const LinearizedTrack *	trk,
		double	trackWeight,
		int	sign,
		int	numVertex,
		bool	isPrimary,
		bool	doFastUpdate ) const

Method to do the Kalman Update of the positions of the fit along the flight axis with a single new track (which belongs to the VxVertexOnJetAxis number numVertex).

Definition at line 140 of file KalmanVertexOnJetAxisUpdator.cxx.

                                                                          {
 
    //linearized track information
 
    const int numrow_toupdate = isPrimary ? 4 : numRow(numVertex);
    
    const VertexPositions & linearizationPositions = candidateToUpdate.getLinearizationVertexPositions();
 
    const Amg::VectorX & initialjetdir = linearizationPositions.position();
    
    //calculate 3-dim position on jet axis from linearizationPositionVector and calculate Jacobian of transformation
    std::pair<Amg::Vector3D,Eigen::Matrix3Xd> PosOnJetAxisAndTransformMatrix = 
      createTransformJacobian(initialjetdir, numVertex, isPrimary,
                              /*truncate=*/ !doFastUpdate); // original fast updator didn't have feature
    
    //only position jacobian changed from A ->oldA
    const AmgMatrix(5,3)& oldA = trk->positionJacobian();
 
    ATH_MSG_DEBUG("the old jacobian xyz d0z0phithetaqoverp vs xyz " << oldA);
 
    //now create the new jacobian which you should use
    Eigen::Matrix<double,5,Eigen::Dynamic> A = oldA*PosOnJetAxisAndTransformMatrix.second;
 
    ATH_MSG_DEBUG("the new jacobian " << A);
    
    const AmgMatrix(5,3)& B = trk->momentumJacobian();
    const AmgVector(5)& trackParameters = trk->expectedParametersAtPCA();
 
    //   the constant term has to be recalculated because of the use of the new parameters
 
    const AmgVector(5) constantTerm = trk->constantTerm()
      + oldA*PosOnJetAxisAndTransformMatrix.first
      - A*initialjetdir.segment(0,numrow_toupdate+1);
 
    //vertex to be updated, needs to be copied
    RecVertexPositions myPosition = candidateToUpdate.getRecVertexPositions();
 
    if (doFastUpdate) {
 
      const AmgSymMatrix(5)& trackParametersCovariance = trk->expectedCovarianceAtPCA();
 
      // old jetvertex in covariance formalism
      const Amg::VectorX & old_vrt_pos = myPosition.position();
      const Amg::MatrixX & old_vrt_cov = myPosition.covariancePosition();
 
      //now create three additional variables for the momenta error
      // (which is nearly infinite... use a big number here... then pay attention to numerical instabilities...)
      AmgSymMatrix(3) old_vrt_cov_momentum; old_vrt_cov_momentum.setZero();
      old_vrt_cov_momentum(0,0) = m_initialMomentumError*m_initialMomentumError;
      old_vrt_cov_momentum(1,1) = m_initialMomentumError*m_initialMomentumError;
      old_vrt_cov_momentum(2,2) = m_initialMomentumError*m_initialMomentumError/10000.;
 
      //R_k_k-1 = V_k + A C_k-1 A_T + B D_k-1 B_T
      AmgSymMatrix(5) old_residual_cov =
        trackParametersCovariance + A * old_vrt_cov * A.transpose() +
        B * old_vrt_cov_momentum * B.transpose();
 
      //the nice thing of the method is that the matrix to invert (old_residual) is just 5x5,
      //which is much better than n.track+5 x n.track+5 !!!
      if (old_residual_cov.determinant() == 0. ) {
        ATH_MSG_WARNING ("The old_residual matrix inversion failed");
        ATH_MSG_WARNING ("same vertex as before is returned");
        return Trk::RecVertexPositions(myPosition);
      }
 
      AmgSymMatrix(5) old_residual_cov_inv = old_residual_cov.inverse().eval();
 
      Eigen::Matrix<double,Eigen::Dynamic,5> Kk1 = old_vrt_cov*A.transpose()*old_residual_cov_inv;
      AmgVector(5) residual_vector=trackParameters-constantTerm-A*old_vrt_pos;
 
      //obtain new position
      Amg::VectorX new_vrt_pos=old_vrt_pos+Kk1*residual_vector;
      Amg::MatrixX new_vrt_cov=old_vrt_cov+Kk1*old_residual_cov*Kk1.transpose()-2.*Kk1*A*old_vrt_cov;
 
      double chi2 = myPosition.fitQuality().chiSquared() +
                    residual_vector.transpose()*old_residual_cov_inv*residual_vector;
 
      //NOT SO NICE: BUT: if there was already a track at this vertex,
      //then add 2, otherwise add 1 (additional parameter in the fit decreases +2 ndf)
      double ndf=myPosition.fitQuality().numberDoF()+sign*2.;
      return Trk::RecVertexPositions(new_vrt_pos,new_vrt_cov,ndf,chi2);
    } 
    
    const AmgSymMatrix(5) & trackParametersWeight  = trk->expectedWeightAtPCA();
    
    if (trackParametersWeight.determinant()<=0) {
      ATH_MSG_WARNING(" The determinant of the inverse of the track covariance matrix is negative: " << trackParametersWeight.determinant());
      if(trk->expectedCovarianceAtPCA().determinant()<=0) {
    ATH_MSG_WARNING(" As well as the determinant of the track covariance matrix: " << trk->expectedCovarianceAtPCA().determinant());
      }
    }
 
    //vertex to be updated, needs to be copied
    myPosition = candidateToUpdate.getRecVertexPositions();
    const Amg::MatrixX & old_full_vrt_cov = myPosition.covariancePosition();
    Eigen::FullPivLU<Amg::MatrixX> lu_decomp(old_full_vrt_cov);
    if(!lu_decomp.isInvertible()){
      ATH_MSG_DEBUG ("The vertex-positions covariance matrix is not invertible");
      ATH_MSG_DEBUG ("The copy of initial vertex returned");
      return Trk::RecVertexPositions(myPosition);
    }
 
    const Amg::VectorX & old_full_vrt_pos = myPosition.position();
    Amg::VectorX old_vrt_pos = old_full_vrt_pos.segment(0,numrow_toupdate+1);
 
    const Amg::MatrixX & old_full_vrt_weight = old_full_vrt_cov.inverse();
 
    if (trackParametersWeight.determinant()<=0) {
      ATH_MSG_WARNING(" The determinant of the track covariance matrix is zero or negative: " << trackParametersWeight.determinant());
    }
 
    Amg::MatrixX old_vrt_weight(numrow_toupdate+1,numrow_toupdate+1);
    old_vrt_weight = old_full_vrt_weight.block(0,0,numrow_toupdate+1,numrow_toupdate+1);
 
    //making the intermediate quantities:
    //W_k = (B^T*G*B)^(-1)
    AmgSymMatrix(3) S = B.transpose()*(trackParametersWeight*B);
 
    if (S.determinant() == 0.0) {
      ATH_MSG_WARNING ("The S matrix is not invertible");
      ATH_MSG_WARNING ("A copy of initial vertex returned");
      return Trk::RecVertexPositions(myPosition);
    }
    S = S.inverse().eval();
 
    //G_b = G_k - G_k*B_k*W_k*B_k^(T)*G_k
    AmgSymMatrix(5) gB = trackParametersWeight - trackParametersWeight*(B*(S*B.transpose()))*trackParametersWeight.transpose();
 
    ATH_MSG_DEBUG("Gain factor obtained: "<<trackParametersWeight*(B*(S*B.transpose()))*trackParametersWeight.transpose());
    ATH_MSG_DEBUG("Resulting Gain Matrix: "<<gB);
   
    //new vertex weight matrix, called "cov" for later inversion
    Amg::MatrixX new_vrt_cov = old_vrt_weight + trackWeight * sign * A.transpose() * ( gB * A );
 
    if (sign<0) {
      ATH_MSG_WARNING(" ATTENTION! Sign is " << sign);
    }
 
    if (new_vrt_cov.determinant()<=0) {
      ATH_MSG_WARNING(std::scientific << "The new vtx weight  matrix determinant is negative: "<< new_vrt_cov.determinant());
    }
 
    // now invert back to covariance
    try {
      // Temporarily remove smart inversion as it seems to cause negative errors from time to time
      // Symmetrize the matix before inversion, to give ride of the precision problem of Eigen.
      new_vrt_cov = (new_vrt_cov+new_vrt_cov.transpose().eval())/2.0;
      smartInvert(new_vrt_cov);
      if (new_vrt_cov.determinant() == 0.0) {
    ATH_MSG_WARNING ("The reduced weight matrix is not invertible, returning copy of initial vertex.");
    const Trk::RecVertexPositions& r_vtx(myPosition);
    return r_vtx;
      }
    }
 
    catch (std::string a) {
      ATH_MSG_WARNING( a << " Previous vertex returned " );
      const Trk::RecVertexPositions& r_vtx(myPosition);
      return r_vtx;
    }
 
    ATH_MSG_DEBUG(" new vertex covariance " << new_vrt_cov);
 
    if (new_vrt_cov.determinant()<=0) {
      ATH_MSG_DEBUG("The new vtx cov. matrix determinant is negative: "<< new_vrt_cov.determinant());
    }
 
    //new vertex position
    Amg::VectorX new_vrt_position =
      new_vrt_cov*(old_vrt_weight * old_vrt_pos + trackWeight * sign * (A.transpose() * (gB *(trackParameters - constantTerm))));
    ATH_MSG_DEBUG(" new position " << new_vrt_position);
 
    //refitted track momentum
    Amg::Vector3D newTrackMomentum = S*B.transpose()*trackParametersWeight*
      (trackParameters - constantTerm - A*new_vrt_position);
    ATH_MSG_DEBUG(" new momentum : " << newTrackMomentum);
 
    //refitted track parameters
    AmgVector(5) refTrackParameters = constantTerm + A * new_vrt_position + B * newTrackMomentum;
 
    //parameters difference
    AmgVector(5) paramDifference = trackParameters - refTrackParameters;
    Amg::VectorX posDifference = new_vrt_position - old_vrt_pos;
 
    double chi2 = myPosition.fitQuality().chiSquared()+
      (paramDifference.transpose()*trackParametersWeight*paramDifference)(0,0)*trackWeight*sign+
      (posDifference.transpose()*old_vrt_weight*posDifference)(0,0); // matrices are 1x1 but make scalar through (0,0)
 
    ATH_MSG_DEBUG(" new chi2 : " << chi2);
 
    //NOT SO NICE: BUT: if there was already a track at this vertex,
    //then add 2, otherwise add 1 (additional parameter in the fit decreases +2 ndf)
    double ndf=myPosition.fitQuality().numberDoF()+sign*2.;
 
    for (int i=0; i<new_vrt_cov.rows(); i++){
      bool negative(false);
      if (new_vrt_cov(i,i)<=0.) {
    ATH_MSG_WARNING(" Diagonal element ("<<i<<","<<i<<") of covariance matrix after update negative: "<<new_vrt_cov(i,i) <<". Giving back previous vertex.");
    negative=true;
      }
      if (negative) {
    const Trk::RecVertexPositions& r_vtx(myPosition);
    return r_vtx;
      }
    }
 
    Amg::VectorX new_full_vrt_pos(old_full_vrt_pos);
    new_full_vrt_pos.segment(0,new_vrt_position.rows()) = new_vrt_position;
    Amg::MatrixX new_full_vrt_cov(myPosition.covariancePosition());
    new_full_vrt_cov.block(0,0,new_vrt_cov.rows(),new_vrt_cov.cols()) = new_vrt_cov;
 
    Trk::RecVertexPositions r_vtx(new_full_vrt_pos,new_full_vrt_cov,ndf, chi2);
    return r_vtx;
    
    //method which avoids inverting huge covariance matrix still needs to be implemented
 
  } //end of position update method

remove()

void Trk::KalmanVertexOnJetAxisUpdator::remove	(	VxTrackAtVertex *	trackToRemove,
		const VxVertexOnJetAxis *	vertexToUpdate,
		VxJetCandidate *	candidateToUpdate ) const

Method removing already added track from the fit along the flight axis, where you specify from what vertex along the jet axis the track should be removed.

Usefull for smoothing and for estimating the compatibility of single vertices with the rest of the fit (the smoothed chi2 of the vertex can be obtained by iteratively removing all its tracks and storing the intermediate result at every time).

Just adding the chi2 of every track involved (as in Fruehwirth's reference) would lead here to misleading results.

Definition at line 59 of file KalmanVertexOnJetAxisUpdator.cxx.

                                                                                      {
    update(trackToRemove,vertexToUpdate,candidateToUpdate,-1,false);
  }

removeWithFastUpdate()

void Trk::KalmanVertexOnJetAxisUpdator::removeWithFastUpdate	(	VxTrackAtVertex *	trackToRemove,
		const VxVertexOnJetAxis *	vertexToUpdate,
		VxJetCandidate *	candidateToUpdate ) const

Definition at line 66 of file KalmanVertexOnJetAxisUpdator.cxx.

                                                                                                    {
    update(trackToRemove,vertexToUpdate,candidateToUpdate,-1,true);
  }

renounce()

std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > AthCommonDataStore< AthCommonMsg< AlgTool > >::renounce ( T & h )

inlineprotectedinherited

Definition at line 380 of file AthCommonDataStore.h.

  {
    h.renounce();
    PBASE::renounce (h);
  }

renounceArray()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::renounceArray ( SG::VarHandleKeyArray & handlesArray )

inlineprotectedinherited

remove all handles from I/O resolution

Definition at line 364 of file AthCommonDataStore.h.

                                                          {
    handlesArray.renounce();
  }

smartInvert()

void Trk::KalmanVertexOnJetAxisUpdator::smartInvert

(

Amg::MatrixX &

new_vrt_cov

)

const

Definition at line 571 of file KalmanVertexOnJetAxisUpdator.cxx.

  {
 
    int numRows = new_vrt_weight.rows();
    if (numRows<=6) {
      if (new_vrt_weight.determinant() ==0) {
        throw std::string("The reduced weight matrix is not invertible. Previous vertex returned ");
      }
      new_vrt_weight = new_vrt_weight.inverse().eval();
      return;
    }
 
    AmgSymMatrix(5) A = new_vrt_weight.block<5,5>(0,0);
    //Eigen::Dynamic we are not sure about the number of Row minus 5, so a dynamic one is given here.
    Eigen::Matrix<double,5,Eigen::Dynamic> B(5,numRows-5);
    B.setZero();
    for (int i=0; i<5; ++i) {
      for (int j=5; j<numRows; ++j) {
        B(i,j-5)=new_vrt_weight(i,j);
      }
    }
 
    // will be of dim (numrows-5)
    Amg::MatrixX D = new_vrt_weight.block(5,5,numRows-5,numRows-5);
 
    // Eigen's diagonal matrices don't have the full MatrixBase members like determinant(),
    // similarity() etc needed later -> therefore do not use them and invert by hand.
    // The inverse of D is a diagonal matrix of the inverted diagonal elements -WL
    // Eigen::DiagonalMatrix<double,Eigen::Dynamic> DdiagINV(D); do something...
    for (unsigned int i=0; i<D.rows(); i++) {
      if ( D(i,i) == 0.) {
        throw std::string("Cannot invert diagonal matrix...");
        break;
      }
 
      //Set the non-diagonal elements of D to be zero, which is as expected; the numerical proplem will change them into non-zero.
      for(unsigned int j=0; j<i; j++){
    D(i,j) = 0;
    D(j,i) = 0;
      }
      D(i,i) = 1./D(i,i);
    }
 
    Amg::MatrixX E = A - B*(D*B.transpose());
    if (E.determinant() == 0.) {
      throw std::string("Cannot invert E matrix...");
    }
    E = E.inverse().eval();
 
    Amg::MatrixX finalWeight(numRows,numRows); 
    finalWeight.setZero();
    finalWeight.block<5,5>(0,0) = E;
    Eigen::Matrix<double,5,Eigen::Dynamic> F = -E*(B*D);
    finalWeight.block(0,5,5,D.rows()) = F;
    finalWeight.block(5,0,D.rows(),5) = F.transpose();
    finalWeight.block(5,5,D.rows(),D.rows()) = 
        D+(D*((B.transpose()*(E*B))*D.transpose()));
 
    new_vrt_weight = finalWeight;
 
    if (new_vrt_weight.determinant()<=0) {
      ATH_MSG_DEBUG("smartInvert() new_vrt_weight FINAL det. is: " << new_vrt_weight.determinant());
    }
 
  }

sysInitialize()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysInitialize ( )

overridevirtualinherited

Perform system initialization for an algorithm.

We override this to declare all the elements of handle key arrays at the end of initialization. See comments on updateVHKA.

Reimplemented in asg::AsgMetadataTool, AthCheckedComponent< AthAlgTool >, AthCheckedComponent<::AthAlgTool >, and DerivationFramework::CfAthAlgTool.

sysStart()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysStart ( )

overridevirtualinherited

Handle START transition.

We override this in order to make sure that conditions handle keys can cache a pointer to the conditions container.

update()

void Trk::KalmanVertexOnJetAxisUpdator::update ( VxTrackAtVertex * trackToUpdate, const VxVertexOnJetAxis * vertexToUpdate, VxJetCandidate * candidateToUpdate, int sign, bool doFastUpdate ) const

private

Internal method to do the Kalman Update.

This is where the real fit iteration is done.

Definition at line 73 of file KalmanVertexOnJetAxisUpdator.cxx.

                                                                              {
    
    //check that all pointers are there...
    if (trackToUpdate==nullptr || vertexToUpdate==nullptr ||
    candidateToUpdate==nullptr || std::abs(sign)!=1) {
      ATH_MSG_WARNING (" Empty pointers then calling fit method update. No fit will be done...");
      return;
    }
 
    //getting tracks at vertex and verify if your trackToUpdate is there
    bool found = false;
    const std::vector<Trk::VxTrackAtVertex*> & tracksAtVertex(vertexToUpdate->getTracksAtVertex());
    for (const auto& track : tracksAtVertex){
      if (track==nullptr) {
    ATH_MSG_WARNING (" Empty pointer in vector of VxTrackAtVertex of the VxVertexOnJetAxis which is being fitted. No fit will be done...");
    return;
      } 
      if (track==trackToUpdate) {
    found=true;
    break;
      }
    }  
    
    if (!found) {
      ATH_MSG_WARNING ("Track was not found in the VxVertexOnJetAxis's list. No fit will be done...");
      return;
    }
    
    double trkWeight = trackToUpdate->weight();
    ATH_MSG_VERBOSE ("Weight is: " << trkWeight);
    
    const Trk::RecVertexPositions & old_vrt = candidateToUpdate->getRecVertexPositions();
    
    if (trackToUpdate->linState()==nullptr) {
      ATH_MSG_WARNING ("Linearized state not associated to track. Aborting fit... " << old_vrt);
      return;
    }
    
    RecVertexPositions fit_vrt;
    if (vertexToUpdate!=candidateToUpdate->getPrimaryVertex()) {
      //if not in primary vertex, then special update
      fit_vrt = positionUpdate(*candidateToUpdate, trackToUpdate->linState(),
                   trkWeight, sign,
                   vertexToUpdate->getNumVertex(), false,
                   doFastUpdate);
    } else {
      //if in primary vertex, use dist=0
      fit_vrt = positionUpdate(*candidateToUpdate, trackToUpdate->linState(),
                   trkWeight, sign,
                   vertexToUpdate->getNumVertex(), true,
                   doFastUpdate);
    }
    
    //forming a final candidate  
    candidateToUpdate->setRecVertexPositions(fit_vrt);
 
    trackToUpdate->setTrackQuality(FitQuality(fit_vrt.fitQuality().chiSquared()-old_vrt.fitQuality().chiSquared(), 2));
    
    
  }//end of update method

updateChi2NdfInfo()

void Trk::KalmanVertexOnJetAxisUpdator::updateChi2NdfInfo	(	VxTrackAtVertex *	trackToUpdate,
		const VxVertexOnJetAxis *	vertexToUpdate,
		VxJetCandidate *	vertexCandidate ) const

Obsolete (will be removed - not in use, however migrated to eigen before finding out -WL) Updates only the chi2 + ndf information, in case the weight matrix is not inverted at every step, so that there is no info about the fitted position at intermediate steps of the fit.

Dropped because the chi2 obtained is not reliable.

Definition at line 362 of file KalmanVertexOnJetAxisUpdator.cxx.

                                                                      {
 
    //check that all pointers are there...
    if (vertexCandidate==nullptr || trackToUpdate==nullptr) {
      ATH_MSG_WARNING( " Empty pointers then calling fit method updateChi2NdfInfo. No update will be done..." );
      return;
    }
 
    //getting tracks at vertex and verify if your trackToUpdate is there
    bool found = false;
    const std::vector<Trk::VxTrackAtVertex*> & tracksAtVertex(vertexToUpdate->getTracksAtVertex());
 
    for (const auto& track : tracksAtVertex) {
      if (track==nullptr) {
    ATH_MSG_WARNING( " Empty pointer in vector of VxTrackAtVertex of the VxVertexOnJetAxis whose chi2 is being updated. No update will be done..." );
    return;
      } 
      if (track==trackToUpdate) {
    found=true;
    break;
      }
    }
    
    if (!found) {
      ATH_MSG_WARNING( " Track was not found in the VxVertexOnJetAxis's list. No update will be done..." );
      return;
    }
 
    double trkWeight = trackToUpdate->weight();
    ATH_MSG_VERBOSE ("Weight is: " << trkWeight);
    
    Trk::RecVertexPositions vertexPositions = vertexCandidate->getRecVertexPositions();
    
    if (trackToUpdate->linState()==nullptr) {
      ATH_MSG_WARNING( "Linearized state not associated to track. Aborting chi2 update... " << vertexPositions );
      return;
    }
 
    double chi2(vertexPositions.fitQuality().chiSquared());
    ATH_MSG_VERBOSE ("old chi2: " << chi2);
 
    if (vertexToUpdate!=vertexCandidate->getPrimaryVertex()) {
      //if not in primary vertex, then special update
      chi2+=calculateTrackChi2(*vertexCandidate, trackToUpdate->linState(),
                   trkWeight,
                   vertexToUpdate->getNumVertex(), false);
    } else {
      chi2+=calculateTrackChi2(*vertexCandidate, trackToUpdate->linState(),
                   trkWeight,
                   vertexToUpdate->getNumVertex(), true);
    }
    
    vertexPositions.setFitQuality(FitQuality(chi2,vertexPositions.fitQuality().numberDoF()));
    
    ATH_MSG_VERBOSE (" new chi2: " << chi2);
 
 
    trackToUpdate->setTrackQuality(FitQuality(chi2-vertexPositions.fitQuality().chiSquared(),
                          trackToUpdate->trackQuality().numberDoF()));
    
    vertexCandidate->setRecVertexPositions(vertexPositions);
    
  }

updateVertexChi2()

void Trk::KalmanVertexOnJetAxisUpdator::updateVertexChi2

(

VxJetCandidate *

vertexCandidate

)

const

Obsolete (will be removed - not in use, however moved to eigen before finding out -WL).

Updates the chi2/ndf of the entire fit with the final information about the position.

Definition at line 481 of file KalmanVertexOnJetAxisUpdator.cxx.

                                                                                           {
 
    //check that all pointers are there...
    if (vertexCandidate==nullptr) { 
      ATH_MSG_WARNING( " Empty pointers then calling chi2 update method updateVertexChi2. No update will be done..." );
      return;
    }
 
    const Trk::RecVertexPositions & OldPosition = vertexCandidate->getConstraintVertexPositions();
    const Amg::VectorX & old_vrt_position = OldPosition.position();
    const Amg::MatrixX & old_vrt_weight = OldPosition.covariancePosition().inverse();
 
    const RecVertexPositions & myPosition = vertexCandidate->getRecVertexPositions();
    const Amg::VectorX & new_vrt_position = myPosition.position();
 
    AmgVector(5) posDifference = (new_vrt_position - old_vrt_position).segment(0,5);
 
    double chi2 = (posDifference.transpose()*old_vrt_weight.block<5,5>(0,0)*posDifference)(0,0);
 
#ifdef Updator_DEBUG
    std::cout << " vertex diff chi2 : " << chi2 << std::endl;
#endif
 
    Trk::RecVertexPositions vertexPositions = vertexCandidate->getRecVertexPositions();
    
    vertexPositions.setFitQuality(FitQuality(vertexPositions.fitQuality().chiSquared()+chi2,
                         vertexPositions.fitQuality().numberDoF()));
 
    vertexCandidate->setRecVertexPositions(vertexPositions);
 
  }

updateVHKA()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::updateVHKA ( Gaudi::Details::PropertyBase & )

inlineinherited

Definition at line 308 of file AthCommonDataStore.h.

                                               {
    // debug() << "updateVHKA for property " << p.name() << " " << p.toString() 
    //         << "  size: " << m_vhka.size() << endmsg;
    for (auto &a : m_vhka) {
      std::vector<SG::VarHandleKey*> keys = a->keys();
      for (auto k : keys) {
        k->setOwner(this);
      }
    }
  }

Member Data Documentation

m_detStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_detStore

privateinherited

Pointer to StoreGate (detector store by default)

Definition at line 393 of file AthCommonDataStore.h.

m_evtStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_evtStore

privateinherited

Pointer to StoreGate (event store by default)

Definition at line 390 of file AthCommonDataStore.h.

m_initialMomentumError

double Trk::KalmanVertexOnJetAxisUpdator::m_initialMomentumError

private

Definition at line 174 of file KalmanVertexOnJetAxisUpdator.h.

m_varHandleArraysDeclared

bool AthCommonDataStore< AthCommonMsg< AlgTool > >::m_varHandleArraysDeclared

privateinherited

Definition at line 399 of file AthCommonDataStore.h.

m_vhka

std::vector<SG::VarHandleKeyArray*> AthCommonDataStore< AthCommonMsg< AlgTool > >::m_vhka

privateinherited

Definition at line 398 of file AthCommonDataStore.h.

---

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

• KalmanVertexOnJetAxisUpdator.h
• KalmanVertexOnJetAxisUpdator.cxx