---

TrLinearExtrapolator.cpp

Go to the documentation of this file.

// Gaudi
#include "GaudiKernel/ToolFactory.h"

// TrEvent
#include "Event/TrState.h"
#include "Event/TrGeneral.h"

// local
#include "TrLinearExtrapolator.h"

// Needed for the creation of TrLinearExtrapolator objects.
static const ToolFactory<TrLinearExtrapolator> s_factory;
const IToolFactory& TrLinearExtrapolatorFactory = s_factory;


TrLinearExtrapolator::TrLinearExtrapolator(const std::string& type,
                                           const std::string& name,
                                           const IInterface* parent):
  GaudiTool(type, name, parent),
  m_F()
{
  declareInterface<ITrExtrapolator>( this );

  // create transport matrix
  m_F = HepMatrix(5, 5, 1);

}

//=============================================================================
// TrLinearExtrapolator destructor.
//=============================================================================
TrLinearExtrapolator::~TrLinearExtrapolator()
{
}

//=============================================================================
// Extrapolate the TrState to z=zNew using the transport matrix  m_F
// (i.e. it performs the mathematical calculation)
//=============================================================================
void TrLinearExtrapolator::extrapolate(TrState* state) const {
  // get reference to the TrState vector and covariance
  HepVector& tX = state -> state();
  HepSymMatrix& tC = state -> covariance();

  // calculate new state
  tX = m_F * tX; // X*F  (can this be done more efficient?)
  tC = tC.similarity(m_F); // F*C*F.T()
}

//=============================================================================
// Propagate a TrState closest to the beam-line
//=============================================================================
StatusCode TrLinearExtrapolator::propagate( TrState* state, 
                                            ParticleID partId ) {

  // check for state
  if ( !state ) {
     return Error( "propagate() should be called with a pointer to a TrState !",
                   StatusCode::FAILURE);
  }

// to a linear approximation the z-position of the closest approach
// to the beam line is given by:
//   zNew = state -> z() - ( x*tx + y*ty ) / ( tx*tx + ty*ty );
  HepVector stateVec = state -> state();
  double zNew = state -> z();
  // the zNew is unchanged if the track is parallel to the beam-line!
  if ( (stateVec[2] != 0) || (stateVec[3] != 0) ) {
    zNew -=   ( stateVec[0]*stateVec[2] + stateVec[1]*stateVec[3] )
            / ( stateVec[2]*stateVec[2] + stateVec[3]*stateVec[3] );
  }
  return propagate( state, zNew, partId );
}

//=============================================================================
// Propagate a TrState to a given z-position
//=============================================================================
StatusCode TrLinearExtrapolator::propagate( TrState* state, 
                                            double zNew,
                                            ParticleID partId ) {

  // check for state
  if ( !state ) {
     return Error( "propagate() should be called with a pointer to a TrState !",
                   StatusCode::FAILURE);
  }

  // create transport matrix
  m_F = HepMatrix(5, 5, 1);

  // check current z-position
  double dz = zNew - state -> z();
  if ( fabs(dz) > TrGeneral::hiTolerance ) {
    m_F[0][2] = dz; // tx*dz
    m_F[1][3] = dz; // ty*dz
    // extrapolate
    extrapolate(state);
    state -> setZ( zNew );
  }
  
  return StatusCode::SUCCESS;
}

//=============================================================================
// Propagate a TrState to the intersection point with a given plane
//=============================================================================
StatusCode TrLinearExtrapolator::propagate( TrState* state,
                                            HepPlane3D plane,
                                            ParticleID partId ) {

  // check for state
  if ( !state ) {
     return Error( "propagate() should be called with a pointer to a TrState !",
                   StatusCode::FAILURE);
  }

  // calculation of the z-position by linear extrapolation to the plane
  // ------------------------------------------------------------------
  // given the plane defined by (a,b,c,d) and the input state of position (x0,y0,z0)
  // and slopes (Tx,Ty), the (x,y,z) of the intersection point verify the set of equations
  //    a*x + b*y + c*z + d = 0
  //    x = x0 + (z-z0)*Tx
  //    x = y0 + (z-z0)*Ty
  HepNormal3D nVec   = plane -> normal();
  HepPoint3D  posVec = state -> position();
  HepVector3D slpVec = state -> slopes();

  double den = nVec.dot( slpVec );

  if ( den != 0 ) {
    slpVec *= ( state -> z() );
    posVec -= slpVec;

    double nom = - ( nVec.dot( posVec ) + plane -> d() );

    double zNew = nom / den;

    return propagate( state, zNew, partId );
  {
  else {
    return StatusCode::FAILURE;
  }
}

//=============================================================================
// Get reference to last used transport matrix.
//=============================================================================
const HepMatrix& TrLinearExtrapolator::transportMatrix() const {
  return m_F;
}

---