TrkHelixRep.cxx

Go to the documentation of this file.

// File and Version Information:
//      $Id: TrkHelixRep.cxx,v 1.4 2011/09/26 00:56:56 zhangy Exp $
//
// Description:
//      Implements helix track representation
//
// Author(s): Steve Schaffner
//------------------------------------------------------------------------
#include "TrkFitter/TrkHelixRep.h"
#include "MdcGeom/BesAngle.h"
#include "MdcGeom/Constants.h"
#include "PatBField/BField.h"
#include "ProxyDict/IfdIntKey.h"
#include "TrkBase/TrkExchangePar.h"
#include "TrkBase/TrkRecoTrk.h"
#include "TrkBase/TrkSimpTraj.h"
using std::endl;
using std::ostream;
 
using MdcPatRec::BField;
 
// Construct from exchange parameters
TrkHelixRep::TrkHelixRep( const TrkExchangePar& inPar, TrkRecoTrk* trk, PdtPid::PidType hypo,
                          const TrkHotList* hots )
    : TrkSimpleRep( hots, trk, hypo ), _traj( inPar ) {}
 
// Construct from exchange parameters, no hits
TrkHelixRep::TrkHelixRep( const TrkExchangePar& inPar, TrkRecoTrk* trk, PdtPid::PidType hypo,
                          int nact, int nsv, int ndc, double chi2, double stFndRng,
                          double endFndRng )
    : TrkSimpleRep( trk, hypo, nact, nsv, ndc, chi2, stFndRng, endFndRng )
    , _traj( inPar.params(), inPar.covariance() ) {}
 
// copy ctor
TrkHelixRep::TrkHelixRep( const TrkHelixRep& old, TrkRecoTrk* trk, PdtPid::PidType hypo )
    : TrkSimpleRep( old, trk, hypo ), _traj( old._traj ) {}
 
TrkHelixRep::~TrkHelixRep() {}
 
TrkHelixRep* TrkHelixRep::clone( TrkRecoTrk* theTrack ) const {
  TrkHelixRep* newRep = new TrkHelixRep( *this, theTrack, this->particleType() );
  newRep->setCurrent( fitCurrent() );
  return newRep;
}
 
TrkHelixRep* TrkHelixRep::cloneNewHypo( PdtPid::PidType hypo ) {
  TrkHelixRep* newRep = new TrkHelixRep( *this, parentTrack(), hypo );
  newRep->setValid( fitValid() );
  newRep->setCurrent( fitCurrent() ); // true for helices, anyway
  return newRep;
}
 
TrkExchangePar TrkHelixRep::helix( double ) const {
  HepVector outPar             = _traj.parameters()->parameter();
  outPar[HelixTraj::phi0Index] = BesAngle( outPar[HelixTraj::phi0Index] ).rad();
  return TrkExchangePar( outPar, _traj.parameters()->covariance() );
}
 
void TrkHelixRep::print( ostream& ostr ) const {
  ostr << "TrkHelixRep "
       << " q: " << charge() << " phi0: " << BesAngle( _traj.phi0() ).rad()
       << " om: " << _traj.omega() << " pt: " << pt() << " p: " << momentum()
       << " d0: " << _traj.d0() << " z0: " << _traj.z0() << " ct: " << _traj.tanDip()
       << "parent track: " << parentTrack()->id() << endl;
}
 
void TrkHelixRep::printAll( ostream& ostr ) const { print( ostr ); }
 
TrkSimpTraj& TrkHelixRep::simpTraj() { return _traj; }
 
const TrkSimpTraj& TrkHelixRep::simpTraj() const { return _traj; }
 
TrkDifTraj& TrkHelixRep::traj() { return _traj; }
 
const TrkDifTraj& TrkHelixRep::traj() const { return _traj; }
 
// Returns arbitrary key, distinguishing this Rep from all other reps
const IfdKey& TrkHelixRep::myKey() const {
  static IfdIntKey _theKey( 3133 );
  return _theKey;
}
 
bool TrkHelixRep::resid( const TrkHitOnTrk* h, double& residual, double& residErr,
                         bool exclude ) const {
  if ( this != h->getParentRep() ) return false;
  static HepVector derivs;
  double           deltaChi;
  TrkErrCode       errCode = h->getFitStuff( derivs, deltaChi );
  // TODO: add caching of derivatives here...
  if ( errCode.failure() ) return false;
  if ( !exclude || !h->isActive() )
  {
    HepSymMatrix c     = _traj.parameters()->covariance().similarityT( derivs );
    double       tErr2 = c( 1, 1 );
    double       rms   = h->hitRms();
    residual           = h->residual();
    double e2          = rms * rms + ( h->isActive() ? -1 : 1 ) * tErr2;
    if ( e2 < 0 ) return false;
    residErr = sqrt( e2 );
  }
  else
  {
    assert( h->isActive() );
    // FIXME: never really checked this is really OK...
    // FIXME: should compare to HelixRefiner...
    HepSymMatrix w = _traj.parameters()->covariance();
    int          ier;
    w.invert( ier );
    if ( ier != 0 ) return false;
    HepVector p0 = _traj.parameters()->parameter();
    HepVector p1 = w * p0;
    p1 -= deltaChi * derivs;
    for ( unsigned i = 0; i < p0.num_row(); ++i )
      for ( unsigned j = 0; j < i + 1; ++j )
      { w.fast( i + 1, j + 1 ) -= derivs[i] * derivs[j]; }
    w.invert( ier );
    if ( ier != 0 ) return false;
    residual = h->residual() + h->hitRms() * dot( derivs, w * p1 - p0 );
    residErr = h->hitRms() * sqrt( w( 1, 1 ) + 1 );
  }
  return true;
}