TTrackBase.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TTrackBase.cxx,v 1.19 2011/09/05 17:32:47 maoh Exp $
//-----------------------------------------------------------------------------
// Filename : TTrackBase.h
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A virtual class for a track class in tracking.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
// #define HEP_SHORT_NAMES
#include "TrkReco/TTrackBase.h"
#include "CLHEP/Alist/CList.h"
#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TMFitter.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TTrackHEP.h"
 
TTrackBase::TTrackBase( const AList<TMLink>& a )
    : _links( a )
    , _mc( 0 )
    ,
    //    _superLayer(0),
    //    _nSuperLayers(0),
    _fitter( 0 )
    , _updated( false )
    , _fitted( false )
    , _time( 0. ) {
  //      _layer[0] = 0;
  //      _layer[1] = 0;
  //      _nLayers[0] = 0;
  //      _nLayers[1] = 0;
  //    update();
}
 
TTrackBase::TTrackBase()
    : _mc( 0 )
    ,
    //    _superLayer(0),
    //    _nSuperLayers(0),
    _fitter( 0 )
    , _updated( true )
    , _fitted( false )
    , _time( 0. ) {
  //      _layer[0] = 0;
  //      _layer[1] = 0;
  //      _nLayers[0] = 0;
  //      _nLayers[1] = 0;
}
 
TTrackBase::~TTrackBase() {}
 
void TTrackBase::dump( const std::string& msg, const std::string& pre ) const {
  bool mc     = ( msg.find( "mc" ) != std::string::npos );
  bool pull   = ( msg.find( "pull" ) != std::string::npos );
  bool flag   = ( msg.find( "flag" ) != std::string::npos );
  bool detail = ( msg.find( "detail" ) != std::string::npos );
  if ( detail ) mc = pull = flag = true;
 
  if ( detail || ( msg.find( "layer" ) != std::string::npos ) )
  {
    if ( !_updated ) update();
  }
  if ( detail || ( msg.find( "hits" ) != std::string::npos ) ) { Dump( _links, msg, pre ); }
}
 
void TTrackBase::update( void ) const {
  _cores.removeAll();
  unsigned n = _links.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink*            l = _links[i];
    const TMDCWireHit& h = *l->hit();
    if ( h.state() & WireHitInvalidForFit ) continue;
    if ( !( h.state() & WireHitFittingValid ) ) continue;
    _cores.append( l );
  }
  _updated = true;
}
 
double TTrackBase::distance( const TMLink& ) const {
  std::cout << "TTrackBase::distance !!! not implemented" << std::endl;
  return 0.;
}
 
int TTrackBase::approach( TMLink& ) const {
  std::cout << "TTrackBase::approach !!! not implemented" << std::endl;
  return -1;
}
 
void TTrackBase::appendByApproach( AList<TMLink>& list, double maxSigma ) {
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    TTrackBase::appendByApproach ... sigma=" << maxSigma << std::endl;
#endif
 
  AList<TMLink> unused;
  unsigned      n = list.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink& l = *list[i];
 
    if ( ( _links.hasMember( l ) ) || ( l.hit()->state() & WireHitUsed ) ) continue;
 
    //...Calculate closest approach...
    int err = approach( l );
    if ( err < 0 )
    {
      unused.append( l );
      continue;
    }
 
    //...Calculate sigma...
    float distance = ( l.positionOnWire() - l.positionOnTrack() ).mag();
    float diff     = fabs( distance - l.drift() );
    float sigma    = diff / l.dDrift();
 
    //...For debug...
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    sigma=" << sigma;
    std::cout << ",dist=" << distance;
    std::cout << ",diff=" << diff;
    std::cout << ",err=" << l.hit()->dDrift() << ",";
    if ( sigma < maxSigma ) std::cout << "ok,";
    else std::cout << "X,";
    l.dump( "mc" );
#endif
 
    //...Make sigma cut...
    if ( sigma > maxSigma )
    {
      unused.append( l );
      continue;
    }
 
    //...OK...
    _links.append( l );
    _updated = false;
    _fitted  = false;
  }
  list.remove( unused );
}
 
void TTrackBase::appendByDistance( AList<TMLink>& list, double sigma ) {
  std::cout << "TTrackBase::appendByDistance !!! not implemented" << std::endl;
  list.removeAll();
}
 
AList<TMLink> TTrackBase::refineMain( double sigma ) {
  AList<TMLink> bad;
  unsigned      n = _links.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    //  std::cout<<" pull of link = "<<_links[i]->pull()<<std::endl;
    if ( _links[i]->pull() > sigma ) bad.append( _links[i] );
  }
  return bad;
}
 
void TTrackBase::refine( AList<TMLink>& list, double sigma ) {
  AList<TMLink> bad = refineMain( sigma );
#ifdef TRKRECO_DEBUG
  std::cout << "    refine ... sigma=" << sigma << ", # of rejected hits=";
  std::cout << bad.length() << std::endl;
#endif
#ifdef TRKRECO_DEBUG
  Dump( bad, "sort pull mc", "        " );
#endif
 
  if ( bad.length() )
  {
    _links.remove( bad );
    list.append( bad );
    _fitted  = false;
    _updated = false;
  }
}
 
void TTrackBase::removeLinks( void ) { _links.removeAll(); }
 
void TTrackBase::refine( double sigma ) {
  AList<TMLink> bad = refineMain( sigma );
  //      for (unsigned i = 0; i < bad.length(); i++) {
  //    const TMDCWireHit * hit = bad[i]->hit();
  //    hit->state(hit->state() | WireHitInvalidForFit);
  //      }
 
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    refine ... sigma=" << sigma << std::endl;
  Dump( bad, "detail sort", "        " );
#endif
 
  if ( bad.length() )
  {
    _fitted  = false;
    _updated = false;
  }
}
int TTrackBase::DropWorst() {
  AList<TMLink> bad;
  int           jbad  = -1;
  double        sigma = 0;
  unsigned      n     = _links.length();
  double        chi2  = 0;
  for ( unsigned i = 0; i < n; i++ )
  {
    // if(g_sigmalyr1)g_sigmalyr1->fill(sqrt(_links[i]->pull()),_links[i]->wire()->layerId());
    chi2 += _links[i]->pull();
    if ( _links[i]->pull() > sigma )
    {
      jbad  = i;
      sigma = _links[i]->pull();
    }
  }
 
  chi2 = chi2 / n;
  // if(sigma>9){bad.append(_links[jbad]);if(g_sigmalyr)g_sigmalyr->fill(sqrt(sigma),_links[jbad]->wire()->layerId());}
  //  if(sigma>20.25&&chi2>2.08){bad.append(_links[jbad]);}
  if ( sigma > 9 ) { bad.append( _links[jbad] ); }
  if ( bad.length() )
  {
    _links.remove( bad );
    _fitted  = false;
    _updated = false;
  }
  return bad.length();
}
 
unsigned TTrackBase::testByApproach( const AList<TMLink>& list, double maxSigma ) const {
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    TTrackBase::testByApproach ... sigma=" << maxSigma << std::endl;
#endif
 
  unsigned nOK = 0;
  unsigned n   = list.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink& l = *list[i];
    nOK += testByApproach( l, maxSigma );
  }
  return nOK;
}
 
unsigned TTrackBase::testByApproach( const TMLink& link, double maxSigma ) const {
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    TTrackBase::testByApproach ... sigma=" << maxSigma << std::endl;
#endif
  TMLink l = link;
 
  //...Calculate closest approach...
  int err = approach( l );
  if ( err < 0 ) return 0;
  //...Calculate sigma...
  float distance = l.distance();
  float diff     = fabs( distance - l.hit()->drift() );
  float sigma    = diff / l.hit()->dDrift();
  l.pull( sigma * sigma );
 
  //...For debug...
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    sigma=" << sigma;
  std::cout << ",dist=" << distance;
  std::cout << ",diff=" << diff << ",";
  if ( sigma < maxSigma ) std::cout << "ok,";
  else std::cout << "X,";
  l.dump( "mc" );
#endif
 
  //...Make sigma cut...
  if ( sigma < maxSigma ) return 1;
 
  return 0;
}
 
//  unsigned
//  TTrackBase::nAxialHits(void) const {
//      unsigned n = _links.length();
//      unsigned a = 0;
//      for (unsigned i = 0; i < n; i++) {
//      if (_links[i]->wire()->axial()) ++a;
//      }
//      return a;
//  }
 
const AList<TMLink>& TTrackBase::links( unsigned mask ) const {
  if ( mask == 0 ) return _links;
 
  std::cout << "TTrackBase::links !!! mask is not supportted yet" << std::endl;
  return _links;
}
 
unsigned TTrackBase::nLinks( unsigned mask ) const {
  unsigned n = _links.length();
  if ( mask == 0 ) return n;
  unsigned nn = 0;
  for ( unsigned i = 0; i < n; i++ )
  {
    const TMDCWireHit& h = *_links[i]->hit();
    if ( h.state() & mask ) ++nn;
  }
  return nn;
}
 
const AList<TMLink>& TTrackBase::cores( unsigned mask ) const {
  if ( mask ) std::cout << "TTrackBase::cores !!! mask is not supported" << std::endl;
  if ( !_updated ) update();
  return _cores;
}
 
unsigned TTrackBase::nCores( unsigned mask ) const {
  if ( mask ) std::cout << "TTrackBase::nCores !!! mask is not supported" << std::endl;
  if ( !_updated ) update();
  return _cores.length();
}
 
//  TMLink *
//  TTrackBase::innerMostLink(unsigned mask) const {
//      if (mask == 0) return InnerMost(_links);
//      unsigned n = _links.length();
//      AList<TMLink> tmp;
//      for (unsigned i = 0; i < n; i++) {
//      const TMDCWireHit & h = * _links[i]->hit();
//      if (h.state() & mask) tmp.append(_links[i]);
//      }
//      return InnerMost(tmp);
//  }
 
//  TMLink *
//  TTrackBase::outerMostLink(unsigned mask) const {
//      if (mask == 0) return OuterMost(_links);
//      unsigned n = _links.length();
//      AList<TMLink> tmp;
//      for (unsigned i = 0; i < n; i++) {
//      const TMDCWireHit & h = * _links[i]->hit();
//      if (h.state() & mask) tmp.append(_links[i]);
//      }
//      return OuterMost(tmp);
//  }
 
int TTrackBase::fit( void ) { return _fitter->fit( *this ); }
 
void TTrackBase::append( TMLink& a ) {
  if ( ( a.hit()->state() & WireHitUsed ) )
  {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "TTrackBase::append !!! " << a.wire()->name()
              << " this is already used by another track!" << std::endl;
#endif
    return;
  }
  if ( _links.hasMember( a ) )
  {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "TTrackBase::append !!! " << a.wire()->name()
              << " this is already included in this track!" << std::endl;
#endif
    return;
  }
  _links.append( a );
  _updated           = false;
  _fitted            = false;
  _fittedWithCathode = false; // added by matsu ( 1999/05/24 )
}
 
void TTrackBase::append( const AList<TMLink>& a ) {
  AList<TMLink> tmp;
  for ( unsigned i = 0; i < a.length(); i++ )
  {
    if ( ( _links.hasMember( a[i] ) ) || ( a[i]->hit()->state() & WireHitUsed ) )
    {
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout << "    TTrackBase::append(list) !!! ";
      std::cout << a[i]->wire()->name();
      std::cout << " Hey!, this is already used! Don't mess with me.";
      std::cout << std::endl;
#endif
      continue;
    }
    else { tmp.append( a[i] ); }
  }
  _links.append( tmp );
  _updated           = false;
  _fitted            = false;
  _fittedWithCathode = false; // added by matsu ( 1999/05/24 )
}
 
const TTrackHEP* const TTrackBase::hep( void ) const {
  unsigned          n = _links.length();
  CAList<TTrackHEP> hepList;
  CList<unsigned>   hepCounter;
  for ( unsigned i = 0; i < n; i++ )
  {
    const TTrackHEP* hep   = _links[i]->hit()->mc()->hep();
    unsigned         nH    = hepList.length();
    bool             found = false;
    for ( unsigned j = 0; j < nH; j++ )
    {
      if ( hepList[j] == hep )
      {
        found = true;
        ++( *hepCounter[j] );
      }
    }
 
    if ( !found )
    {
      hepList.append( hep );
      unsigned c = 0;
      hepCounter.append( c );
    }
  }
 
  _nHeps       = hepList.length();
  _hep         = 0;
  unsigned max = 0;
  for ( unsigned i = 0; i < _nHeps; i++ )
  {
    if ( ( *hepCounter[i] ) > max )
    {
      max  = ( *hepCounter[i] );
      _hep = hepList[i];
    }
  }
 
  return _hep;
}
 
unsigned TTrackBase::nHeps( void ) const {
  hep();
  return _nHeps;
}