T3DLine.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: T3DLine.cxx,v 1.9 2010/03/31 09:58:59 liucy Exp $
//-----------------------------------------------------------------------------
// Filename : T3DLine.cc
// Section  : Tracking
// Owner    : Kenji Inami
// Email    : inami@bmail.kek.jp
//-----------------------------------------------------------------------------
// Description : A class to represent a 3D line in tracking.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#include "TrkReco/T3DLine.h"
#include "TrkReco/T3DLineFitter.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TTrack.h"
// #include "TrkReco/TMDCWire.h"
 
const T3DLineFitter T3DLine::_fitter = T3DLineFitter( "T3DLine Default fitter" );
 
T3DLine::T3DLine()
    : TTrackBase()
    , _dr( 0 )
    , _phi0( 0 )
    , _dz( 0 )
    , _tanl( 0 )
    , _chi2( 0 )
    , _ndf( 0 )
    , _Ea( SymMatrix( 4, 0 ) )
    , _pivot( ORIGIN ) {
 
  _cos_phi0 = 1;
  _sin_phi0 = 0;
 
  //...Set a default fitter...
  fitter( &T3DLine::_fitter );
 
  _fitted            = false;
  _fittedWithCathode = false;
}
 
T3DLine::T3DLine( const TTrack& a )
    : TTrackBase( (TTrackBase&)a )
    , _dr( a.helix().dr() )
    , _phi0( a.helix().phi0() )
    , _dz( a.helix().dz() )
    , _tanl( a.helix().tanl() )
    , _chi2( 0 )
    , _ndf( 0 )
    , _Ea( SymMatrix( 4, 0 ) )
    , _pivot( a.helix().pivot() ) {
 
  _cos_phi0 = cos( _phi0 );
  _sin_phi0 = sin( _phi0 );
 
  //...Set a default fitter...
  fitter( &T3DLine::_fitter );
 
  _fitted            = false;
  _fittedWithCathode = false;
}
 
T3DLine::~T3DLine() {}
 
T3DLine::T3DLine( const T3DLine& a )
    : TTrackBase( (TTrackBase&)a )
    , _dr( a.dr() )
    , _phi0( a.phi0() )
    , _dz( a.dz() )
    , _tanl( a.tanl() )
    , _chi2( a.chi2() )
    , _ndf( a.ndf() )
    , _Ea( a.Ea() )
    , _pivot( a.pivot() ) {
 
  _cos_phi0 = cos( _phi0 );
  _sin_phi0 = sin( _phi0 );
 
  //...Set a default fitter...
  fitter( &T3DLine::_fitter );
 
  _fitted            = false;
  _fittedWithCathode = false;
}
 
double T3DLine::dr( void ) const { return _dr; }
 
double T3DLine::phi0( void ) const { return _phi0; }
 
double T3DLine::dz( void ) const { return _dz; }
 
double T3DLine::tanl( void ) const { return _tanl; }
 
double T3DLine::cosPhi0( void ) const { return _cos_phi0; }
 
double T3DLine::sinPhi0( void ) const { return _sin_phi0; }
 
const HepPoint3D& T3DLine::pivot( void ) const { return _pivot; }
 
Vector T3DLine::a( void ) const {
  Vector ta( 4 );
  ta[0] = _dr;
  ta[1] = _phi0;
  ta[2] = _dz;
  ta[3] = _tanl;
  return ( ta );
}
 
const SymMatrix& T3DLine::Ea( void ) const { return ( _Ea ); }
 
Helix T3DLine::helix( void ) const {
  HepVector a( 5 );
  a[0] = _dr;
  a[1] = _phi0;
  a[2] = 1e-10;
  a[3] = _dz;
  a[4] = _tanl;
  Helix _helix( _pivot, a );
  return _helix;
}
 
unsigned T3DLine::ndf( void ) const { return _ndf; }
 
double T3DLine::chi2( void ) const { return _chi2; }
 
double T3DLine::reducedchi2( void ) const {
  if ( _ndf == 0 )
  {
    std::cout << "error at T3DLine::reducedchi2  ndf=0" << std::endl;
    return 0;
  }
  return ( _chi2 / _ndf );
}
 
HepPoint3D T3DLine::x( double t ) const {
  double     tx = _pivot.x() + _dr * _cos_phi0 - t * _sin_phi0;
  double     ty = _pivot.y() + _dr * _sin_phi0 + t * _cos_phi0;
  double     tz = _pivot.z() + _dz + t * _tanl;
  HepPoint3D p( tx, ty, tz );
  return p;
}
 
HepPoint3D T3DLine::x0( void ) const {
  double     tx = _pivot.x() + _dr * _cos_phi0;
  double     ty = _pivot.y() + _dr * _sin_phi0;
  double     tz = _pivot.z() + _dz;
  HepPoint3D p( tx, ty, tz );
  return p;
}
HepVector3D T3DLine::k( void ) const {
  HepPoint3D p( -_sin_phi0, _cos_phi0, _tanl );
  return p;
}
 
const HepPoint3D& T3DLine::pivot( const HepPoint3D& newpivot ) {
  double dr = ( _pivot.x() - newpivot.x() ) * _cos_phi0 +
              ( _pivot.y() - newpivot.y() ) * _sin_phi0 + _dr;
  double dz = _pivot.z() - newpivot.z() + _dz +
              _tanl * ( ( _pivot.x() - newpivot.x() ) * _sin_phi0 +
                        ( newpivot.y() - _pivot.y() ) * _cos_phi0 );
  _dr    = dr;
  _dz    = dz;
  _pivot = newpivot;
  return _pivot;
}
 
void T3DLine::set( const HepPoint3D& t_pivot, double t_dr, double t_phi0, double t_dz,
                   double t_tanl ) {
 
  _pivot    = t_pivot;
  _dr       = t_dr;
  _phi0     = t_phi0;
  _dz       = t_dz;
  _tanl     = t_tanl;
  _cos_phi0 = cos( _phi0 );
  _sin_phi0 = sin( _phi0 );
}
 
Vector T3DLine::a( const Vector& ta ) {
  _dr       = ta[0];
  _phi0     = ta[1];
  _dz       = ta[2];
  _tanl     = ta[3];
  _cos_phi0 = cos( _phi0 );
  _sin_phi0 = sin( _phi0 );
  return ( ta );
}
 
const SymMatrix& T3DLine::Ea( const SymMatrix& tEa ) {
  _Ea = tEa;
  return ( _Ea );
}
 
int T3DLine::approach( TMLink& l, bool doSagCorrection ) const {
 
  const TMDCWire& w                    = *l.wire();
  HepPoint3D      xw                   = w.xyPosition();
  HepPoint3D      wireBackwardPosition = w.backwardPosition();
  HepVector3D     v                    = w.direction();
 
  HepPoint3D onWire, onTrack;
 
  if ( approach_line( wireBackwardPosition, v, onWire, onTrack ) < 0 ) return ( -1 );
 
  // onWire,onTrack filled
 
  if ( !doSagCorrection )
  {
    l.positionOnWire( onWire );
    l.positionOnTrack( onTrack );
    return ( 0 ); // no sag correction
  }
  // Sag correction
  //   loop for sag correction
  double onWire_y = onWire.y();
  double onWire_z = onWire.z();
 
  unsigned nTrial = 1;
  while ( nTrial < 100 )
  {
    w.wirePosition( onWire_z, xw, wireBackwardPosition, (HepVector3D&)v );
    if ( approach_line( wireBackwardPosition, v, onWire, onTrack ) < 0 ) return ( -1 );
    if ( fabs( onWire_y - onWire.y() ) < 0.0001 ) break; // |dy|< 1 micron
    onWire_y = onWire.y();
    onWire_z = onWire.z();
 
    nTrial++;
  }
 
  l.positionOnWire( onWire );
  l.positionOnTrack( onTrack );
  return ( nTrial );
}
 
int T3DLine::approach_line( const HepPoint3D& w0, const HepVector3D& v, HepPoint3D& onLine,
                            HepPoint3D& onTrack ) const {
  //  line = [w0] + s * [v]    -> [onLine]
  //  trk  = [x0] + t * [k]    -> [onTrack]
  //  if [v]//[k] then return(-1) error
 
  const HepVector3D k   = this->k();
  const double      v_k = v.dot( k );
  const double      v_2 = v.mag2();
  const double      k_2 = k.mag2();
  const double      tk  = v_k * v_k - v_2 * k_2;
  if ( tk == 0 ) return ( -1 );
 
  const HepPoint3D x0 = this->x0();
  const HepPoint3D dx = x0 - w0;
  const double     t  = dx.dot( v_2 * k - v_k * v ) / tk;
  const double     s  = dx.dot( v_k * k - k_2 * v ) / tk;
 
  onLine  = w0 + s * v;
  onTrack = x0 + t * k;
  return ( 0 );
}
 
int T3DLine::approach_point( const HepPoint3D& p0, HepPoint3D& onTrack ) const {
  //  trk  = [x0] + t * [k]    -> [onTrack]
  //  if [v]//[k] then return(-1) error
 
  const HepVector3D k  = this->k();
  const HepPoint3D  x0 = this->x0();
  const HepPoint3D  dx = p0 - x0;
  const double      t  = dx.dot( k ) / k.mag2();
 
  onTrack = x0 + t * k;
  return ( 0 );
}