Ext_errmx.cxx

Go to the documentation of this file.

// File: Ext_errmx.cc
//
// Description: Handle error matrix( x, y, z, px, py, pz ).
//              The used coordinate system is the cartesian 
//              BESIII coordinate system.The error matrix is 
//              now calculated in the Track Coordinate System( TCS ) 
//              instead of the BESIII Coordinate System( BCS )
//              in the previous version.
//              
//              Modified from BELLE:
//                                   Creation: 08-Jan-1998
//                                   Version: 05-Mar-1999
//              by Wang Liangliang
// 
// Date: 2005.3.30                          
//
 
#include <iostream>
 
#include "TrkExtAlg/Ext_errmx.h"
 
extern bool Ext_err_valid( bool msg, const HepSymMatrix &error, 
              const int dimension );
static const double Eps(1.0e-12);   // Infinitesimal number.
static const double Infinite(1.0e+12);  // Infinite number.
 
// Default constructor
 
Ext_errmx::Ext_errmx() : m_err(Ndim_err,0), m_err3(3,0), m_R(3,3,0), 
    m_err2(2), m_valid(0) {}
 
// Copy constructor
 
Ext_errmx::Ext_errmx( const Ext_errmx &err ) : m_err(err.m_err),
    m_err3(err.m_err3), m_R(err.m_R), m_err2(err.m_err2), 
    m_valid(err.m_valid) {}
 
Ext_errmx::Ext_errmx( const HepSymMatrix &err ) : m_err(err), 
    m_err3(3,0), m_R(3,3,0), m_err2(2), m_valid(1) {}
 
/*
  Put the error matrix.
*/
 
void Ext_errmx::put_err( const double error[] )
{
  int ne = 0;
 
  for( int i = 1; i <= m_err.num_col(); i++ ){
    for( int j = 1; j <= i; j++ ){
      m_err.fast( i, j ) = error[ ne++ ];
    }
  }
  m_valid = 1;
}
 
/*
  Setup m_err3 and m_R for get_plane_err(s).
*/
void Ext_errmx::set_plane_errs( const Hep3Vector &nx, const Hep3Vector &ny, 
                   const Hep3Vector &nz ) const
{
// Setup the rotation matrix.
 
  m_R( 1, 1 ) = nx.x();
  m_R( 1, 2 ) = nx.y();
  m_R( 1, 3 ) = nx.z();
  m_R( 2, 1 ) = ny.x();
  m_R( 2, 2 ) = ny.y();
  m_R( 2, 3 ) = ny.z();
  m_R( 3, 1 ) = nz.x();
  m_R( 3, 2 ) = nz.y();
  m_R( 3, 3 ) = nz.z();
 
// Get the 3x3 sub matrix and Rotate the error matrix.
 
  m_err3 = (m_err.sub( 1, 3 )).similarity( m_R );
 
}
 
/*
  Get the projection of error on the plane along the readout direction. 
*/
 
double Ext_errmx::get_plane_err( const Hep3Vector &np, const Hep3Vector &nr ) 
    const
// ( Inputs )
// np -- Unit vector to the direction of the track.
// nr -- Readout direction unit vector on the plane.
//
// ( Outputs )
// HepVector[0] = Error( sigma ) along the readout direction.
// HepVector[1] = Error( sigma ) along the direction perp. to 
//          the readout direction.
{
 
// Check the validity of the error matrix.
 
  if( !(*this).valid( 0 ) ){
    //std::cout << "%WARNING at Ext_get_plane_err: You are trying to calculate"
//  << " error \n           using an invalid error matrix." << std::endl;
  }
 
// Construct 3 TCS axes.
 
// x: nx --- unit vector which is perp to np and on the plane of nr and np.
// y: nz x nx
// z: nz( = np ) --- direction of the track.
 
  double nr_np( nr*np );
  double denom( 1.0 - nr_np*nr_np );
  double error( 0.0 );
 
  if( denom > Eps ){    // Track is not parallel to the readout direction.
 
    double fac( 1.0 / sqrt( denom ) );
    Hep3Vector nx( ( nr - nr_np * np ) * fac );
    Hep3Vector ny( np.cross( nx ) );
 
    (*this).set_plane_errs( nx, ny, np );
 
    double sigma2( m_err3( 1, 1 ) );    // Error along nx.
    if( sigma2 > 0 ){
      error = sqrt( sigma2 ) * fac; // Error projection.
    }
 
  } else {  // Track is parallel to the readout direction.
 
    error = Infinite;
 
  }
  return( error );
}
 
/*
  Get a projection of error on the plane along the readout direction
  and its perpendicular direction on the readout plane. 
*/
 
const HepVector &Ext_errmx::get_plane_errs( const Hep3Vector &np, 
        const Hep3Vector &nr, const Hep3Vector &nt ) const
// ( Inputs )
// np -- Unit vector to the direction of the track.
// nr -- Readout direction unit vector on the plane.
// nt -- Unit vector that is perp. to nr on the readout plane.
//
// ( Outputs )
// HepVector[0] = Error( sigma ) along the readout direction.
// HepVector[1] = Error( sigma ) along the direction perp. to 
//          the readout direction.
{
 
// Check the validity of the error matrix.
 
  if( !(*this).valid( 1 ) ){
//    std::cout << "%WARNING at Ext_get_plane_errs: You are trying to calculate"
//  << " error \n           using an invalid error matrix." << std::endl;
  }
 
  double nr_np( nr*np );
  double denom_r( 1.0 - nr_np*nr_np );
 
  if( denom_r > Eps ){  // nr is not parallel to the track direction: np.
 
    double nt_np( nt*np );
    double denom_t( 1.0 - nt_np*nt_np );
    double fac_r( 1.0 / sqrt( denom_r ) );
    Hep3Vector nx( ( nr - nr_np * np ) * fac_r );
    Hep3Vector ny( np.cross( nx ) );
 
    (*this).set_plane_errs( nx, ny, np );
 
    double sigma2( m_err3( 1, 1 ) );    // Error along nx.
    if( sigma2 > 0 ){
      m_err2( 1 ) = sqrt( sigma2 ) * fac_r; // Error projection.
    } else {
      m_err2( 1 ) = 0.0;
    }
 
    if( denom_t > Eps ){   // nt is not parallel to the track direction: np.
      double fac_t( 1.0 / sqrt( denom_t ) );
      sigma2 = m_err3( 2, 2 );
      if( sigma2 > 0 ){
    m_err2( 2 ) = sqrt( sigma2 ) * fac_t;   // Error projection.
      } else {
    m_err2( 2 ) = 0.0;
      }
    } else {    // nt is parallel to the track direction: np.
      m_err2( 2 ) = (*this).get_plane_err( np, nt );
    }
  } else {  // nr is parallel to the track direction: np.
    m_err2( 1 ) = (*this).get_plane_err( np, nr );
    m_err2( 2 ) = (*this).get_plane_err( np, nt );
  }
  return( m_err2 );
}
 
/*
  Get the 2D projected track error perpendicular to a given vector at the 
  current point. pv: momentum vector, for example.
  view=1. Hep3Vector(1)= error(sigma) along the direction 
        which is perpendicular to pv on the xy plane.
          Hep3Vector(2)= error(sigma) along the direction 
        which is (pv) x (Hep3Vector(1) direction).
  view=2. Hep3Vector(1)= error(sigma) along the direction 
        which is perpendicular to pv on the zy plane.
          Hep3Vector(2)= error(sigma) along the direction 
        which is (pv) x (Hep3Vector(1) direction).
  view=3. Hep3Vector(1)= error(sigma) along the direction 
        which is perpendicular to pv on the zx plane.
          Hep3Vector(2)= error(sigma) along the direction 
        which is (pv) x (Hep3Vector(1) direction).
*/
 
const Hep3Vector * Ext_errmx::get_tvs( const int view, 
                     const Hep3Vector &pv ) const
// ( Inputs )
// view  -- 2D projection view. = 1 xy, =2 zy, =3 zx.
// pv    -- A vector for which the perpendicular error will be calculated,
//      for example, momentum of the track.
{
 
  Hep3Vector    np( pv.unit() );
  Hep3Vector    nr;
 
  switch( view ){
 
  case 1:   // xy view
 
    if( np.x() != 0 || np.y() != 0 ){
      nr.setX(  np.y() );
      nr.setY( -np.x() );
      nr = nr.unit();
    } else {        // Pointing to z-direction.
      nr.setX( 1 );
    }
    break;
 
  case 2:   // zy view
 
    if( np.y() != 0 || np.z() != 0 ){
      nr.setY( -np.z() );
      nr.setZ(  np.y() );
      nr = nr.unit();
    } else {        // Pointing to x-direction.
      nr.setZ( 1 );
    }
    break;
 
  case 3:   // zx view
 
    if( np.z() != 0 || np.x() != 0 ){
      nr.setX( -np.z() );
      nr.setZ(  np.x() );
      nr = nr.unit();
    } else {        // Pointing to z-direction.
      nr.setZ( 1 );
    }
    break;
  } /* End of switch */
 
  Hep3Vector nt( np.cross( nr ) );
  const HepVector & err_v = (*this).get_plane_errs( np, nr, nt );
  *m_nv     = err_v[0]*nr;
  *(m_nv+1) = err_v[1]*nt;
  return( m_nv );
}
 
/*
  Get the 2D projected track error perpendicular to a given vector at the 
  current point. pv: momentum vector, for example.
  Hep3Vector(1)= error(sigma) along the direction which is
        perpendicular to pv on the plane formed by pv and z-axis.
  Hep3Vector(2)= error(sigma) along the direction which is
        (pv) x (Hep3Vector(1) direction).
*/
 
const Hep3Vector * Ext_errmx::get_tvs( const Hep3Vector &pv ) const
// ( Inputs )
// pv    -- A vector for which the perpendicular error will be calculated,
//      for example, momentum of the track.
{
 
  Hep3Vector    np( pv.unit() );
  Hep3Vector    nz( 0.0, 0.0, 1.0 );
  Hep3Vector    nt( (nz.cross(np)).unit() );
  Hep3Vector    nr( nt.cross(np) );
 
  const HepVector & err_v = (*this).get_plane_errs( np, nr, nt );
  *m_nv     = err_v[0]*nr;
  *(m_nv+1) = err_v[1]*nt;
  return( m_nv );
 
}
 
/*
  valid( msg ). Check the validity of the diagonal elements.
*/
 
bool Ext_errmx::valid( bool msg ) const 
{
  if( m_valid ){
    if( Ext_err_valid( msg, m_err, 6 ) ){
      return 1;
    } else {
      m_valid = 0;
      return 0;
    }
  } else {
    return 0;
  }
}
 
// Assign "=" operator
 
Ext_errmx &Ext_errmx::operator=( const Ext_errmx &err )
{
  if( this != &err ){
    m_err   = err.m_err;
    m_err3  = err.m_err3;
    m_err2  = err.m_err2;
    m_R     = err.m_R;
    m_valid = err.m_valid;
    *m_nv   = *err.m_nv;
    *(m_nv+1)   = *(err.m_nv+1);
  }
  return *this;
}
 
/*
 ostream "<<" operator
*/
std::ostream &operator<<( std::ostream &s, const Ext_errmx &err )
{
  s << " m_valid: " << err.m_valid << '\n'
    << "m_err: " << err.m_err << " m_err3: " << err.m_err3 
    << " m_R: " << err.m_R << " m_err2: " << err.m_err2
    << " *m_nv: " << *err.m_nv 
    << " *(m_nv+1): " << *(err.m_nv+1)
    << std::endl;
  return s;
}