TMLink.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TMLink.cxx,v 1.15 2011/10/08 01:56:15 maoh Exp $
//-----------------------------------------------------------------------------
// Filename : TMLink.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class to relate TMDCWireHit and TTrack objects.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#include "TrkReco/TMLink.h"
#include "CLHEP/Alist/ConstAList.h"
#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TTrackHEP.h"
 
// zangsl 040518 : move inlined function here
const TMDCWire* const TMLink::wire( void ) const {
  if ( _hit ) return _hit->wire();
  return NULL;
}
 
const HepPoint3D& TMLink::xyPosition( void ) const { return _hit->wire()->xyPosition(); }
// zangsl 040518 end
 
/*TMLink::TMLink(TTrack * t, const TMDCWireHit * h, const HepPoint3D & p)
: _track(t),
  _hit(h),
  _dPhi(0),
  _leftRight(0),
  _pull(0),
  _position(p),
  _link(0),
  _fit2D(0) {
//  _usecathode(0) {
    if (h) {
        _drift[0] = h->drift(0);
        _drift[1] = h->drift(1);
        _dDrift[0] = h->dDrift(0);
        _dDrift[1] = h->dDrift(1);
    }
    else {
        _drift[0] = 0.;
        _drift[1] = 0.;
        _dDrift[0] = 0.;
        _dDrift[1] = 0.;
    }
 
    for (unsigned i = 0; i < 6; ++i)
        _neighbor[i] = NULL;
    for (unsigned i = 0; i < 4; ++i)
      _arcZ[i];
 
    if (h) {
        _onTrack = _onWire = h->xyPosition();
    }
}
*/
// for Tsf
TMLink::TMLink( TTrack* t, const TMDCWireHit* h, const HepPoint3D& p, const HepPoint3D& d,
                double dr )
    : _track( t )
    , _hit( h )
    , _dPhi( 0 )
    , _leftRight( 0 )
    , _pull( 0 )
    , _position( p )
    , _positionD( d )
    , _link( 0 )
    , _fit2D( 0 )
    , _tsfTag( 0 ) {
  if ( h )
  {
    _cDrift[0] = dr; // after conformal transformation.
    _cDrift[1] = dr;
    _drift[0]  = h->drift( 0 );
    _drift[1]  = h->drift( 1 );
    _dDrift[0] = h->dDrift( 0 );
    _dDrift[1] = h->dDrift( 1 );
  }
  else
  {
    _cDrift[0] = 0.;
    _cDrift[1] = 0.;
    _drift[0]  = 0.;
    _drift[1]  = 0.;
    _dDrift[0] = 0.;
    _dDrift[1] = 0.;
  }
 
  for ( unsigned i = 0; i < 6; ++i ) _neighbor[i] = NULL;
  for ( unsigned i = 0; i < 4; ++i ) _arcZ[i];
 
  if ( h ) { _onTrack = _onWire = h->xyPosition(); }
}
 
TMLink::TMLink( const TMLink& l )
    : _track( l._track )
    , _hit( l._hit )
    , _onTrack( l._onTrack )
    , _onWire( l._onWire )
    , _position( l._position )
    , _positionD( l._positionD )
    , _dPhi( l._dPhi )
    , _leftRight( l._leftRight )
    , _pull( l._pull )
    , _link( l._link )
    , _distance( l._distance )
    ,
    // addition by matsu ( 1999/07/05 )
    //  _mclust(l._mclust ),
    //  _usecathode(l._usecathode ),
    // end of addition
    _fit2D( l._fit2D )
    , _tsfTag( l._tsfTag ) {
  _drift[0]  = l._drift[0];
  _drift[1]  = l._drift[1];
  _dDrift[0] = l._dDrift[0];
  _dDrift[1] = l._dDrift[1];
  _cDrift[0] = l._cDrift[0];
  _cDrift[1] = l._cDrift[1];
  for ( unsigned i = 0; i < 6; ++i ) _neighbor[i] = l._neighbor[i];
  for ( unsigned i = 0; i < 4; ++i ) _arcZ[i] = l._arcZ[i];
}
 
TMLink::~TMLink() {}
 
void TMLink::dump( const std::string& msg, const std::string& pre ) const {
  std::cout << pre;
  if ( _track ) std::cout << "track#=,";
  if ( _hit ) { _hit->dump( msg ); }
}
 
unsigned NLayers( const AList<TMLink>& list ) {
  unsigned l0 = 0;
  unsigned l1 = 0;
  unsigned n  = list.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    unsigned id = list[i]->wire()->layerId();
    if ( id < 32 ) l0 |= ( 1 << id );
    else l1 |= ( 1 << ( id - 32 ) );
  }
 
  unsigned l = 0;
  for ( unsigned i = 0; i < 32; i++ )
  {
    if ( l0 & ( 1 << i ) ) ++l;
    if ( l1 & ( 1 << i ) ) ++l;
  }
  return l;
}
 
void NHits( const AList<TMLink>& links, unsigned nHits[43] ) {
  for ( unsigned i = 0; i < 43; i++ ) nHits[i] = 0;
  unsigned nLinks = links.length();
  for ( unsigned i = 0; i < nLinks; i++ ) ++nHits[links[i]->wire()->layerId()];
}
 
void NHitsSuperLayer( const AList<TMLink>& links, unsigned nHits[11] ) {
  for ( unsigned i = 0; i < 11; i++ ) nHits[i] = 0;
  unsigned nLinks = links.length();
  for ( unsigned i = 0; i < nLinks; i++ ) ++nHits[links[i]->wire()->superLayerId()];
}
 
void Dump( const CAList<TMLink>& links, const std::string& msg, const std::string& pre ) {
  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 sort   = ( msg.find( "sort" ) != std::string::npos );
  bool stereo = ( msg.find( "stereo" ) != std::string::npos );
  bool detail = ( msg.find( "detail" ) != std::string::npos );
  bool pos    = ( msg.find( "position" ) != std::string::npos );
  if ( detail ) mc = pull = flag = sort = true;
 
  CAList<TMLink> tmp = links;
  if ( sort ) tmp.sort( SortByWireId );
  unsigned n       = tmp.length();
  unsigned nForFit = 0;
#define MCC_MAX 1000
  unsigned MCC0[MCC_MAX];
  unsigned MCC1[MCC_MAX];
  for ( unsigned i = 0; i < MCC_MAX; i++ )
  {
    MCC0[i] = 0;
    MCC1[i] = 0;
  }
  bool MCCOverFlow = false;
 
  std::cout << pre;
  for ( unsigned i = 0; i < n; i++ )
  {
    const TMLink& l = *tmp[i];
    std::cout << l.wire()->name();
 
    double   a    = l.pull();
    unsigned mcId = 0;
    if ( mc )
      if ( l.hit()->mc() )
        if ( l.hit()->mc()->hep() ) mcId = l.hit()->mc()->hep()->id();
    if ( pull ) { std::cout << "[" << a << "]"; }
    if ( mc )
    {
      std::cout << "(" << mcId << ")";
      if ( mcId < MCC_MAX )
      {
        ++MCC0[mcId];
        if ( l.hit()->state() & WireHitFittingValid )
        {
          if ( !( l.hit()->state() & WireHitInvalidForFit ) ) ++MCC1[mcId];
        }
      }
      else { MCCOverFlow = true; }
    }
    if ( flag )
    {
      if ( l.hit()->state() & WireHitFindingValid ) std::cout << "o";
      if ( l.hit()->state() & WireHitFittingValid )
      {
        std::cout << "+";
        if ( !( l.hit()->state() & WireHitInvalidForFit ) ) ++nForFit;
      }
      if ( l.hit()->state() & WireHitInvalidForFit ) std::cout << "x";
    }
    if ( stereo ) { std::cout << "{" << l.leftRight() << "," << l.zStatus() << "}"; }
    if ( pos ) { std::cout << ",pos=" << l.position(); }
    std::cout << ",";
  }
  std::cout << " " << n << " l(s)";
  if ( flag ) std::cout << ", fv " << nForFit << " l(s)";
  if ( mc )
  {
    unsigned nMC = 0;
    std::cout << ", mc";
    for ( unsigned i = 0; i < MCC_MAX; i++ )
    {
      if ( MCC0[i] > 0 )
      {
        ++nMC;
        std::cout << i << ":" << MCC0[i] << ",";
      }
    }
    std::cout << " total " << nMC << " contributions";
    if ( flag )
    {
      nMC = 0;
      std::cout << ", fv mc";
      for ( unsigned i = 0; i < MCC_MAX; i++ )
      {
        if ( MCC1[i] > 0 )
        {
          ++nMC;
          std::cout << i << ":" << MCC1[i] << ",";
        }
      }
      std::cout << " total " << nMC << " contributions";
    }
 
    if ( MCCOverFlow ) std::cout << "(counter overflow)";
  }
  std::cout << std::endl;
}
 
void Dump( const TMLink& link, const std::string& msg, const std::string& pre ) {
  CAList<TMLink> tmp;
  tmp.append( link );
  Dump( tmp, msg, pre );
}
 
unsigned NStereoHits( const AList<TMLink>& links ) {
  unsigned nLinks = links.length();
  unsigned n      = 0;
  for ( unsigned i = 0; i < nLinks; i++ )
    if ( links[i]->wire()->stereo() ) ++n;
  return n;
}
 
unsigned NAxialHits( const AList<TMLink>& links ) {
  unsigned nLinks = links.length();
  unsigned n      = 0;
  for ( unsigned i = 0; i < nLinks; i++ )
    if ( links[i]->wire()->axial() ) ++n;
  return n;
}
 
AList<TMLink> AxialHits( const AList<TMLink>& links ) {
  AList<TMLink> a;
  unsigned      n = links.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    if ( links[i]->wire()->axial() ) a.append( links[i] );
  }
  return a;
}
 
AList<TMLink> StereoHits( const AList<TMLink>& links ) {
  AList<TMLink> a;
  unsigned      n = links.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    if ( !links[i]->wire()->axial() ) a.append( links[i] );
  }
  return a;
}
 
TMLink* InnerMost( const AList<TMLink>& a ) {
  unsigned n     = a.length();
  unsigned minId = 9999;
  TMLink*  t     = 0;
  for ( unsigned i = 0; i < n; i++ )
  {
    unsigned id = a[i]->wire()->id();
    if ( id < minId )
    {
      minId = id;
      t     = a[i];
    }
  }
  return t;
}
 
TMLink* OuterMost( const AList<TMLink>& a ) {
  unsigned n     = a.length();
  unsigned maxId = 0;
  TMLink*  t     = 0;
  for ( unsigned i = 0; i < n; i++ )
  {
    unsigned id = a[i]->wire()->id();
    if ( id > maxId )
    {
      maxId = id;
      t     = a[i];
    }
  }
  return t;
}
 
void SeparateCores( const AList<TMLink>& input, AList<TMLink>& cores,
                    AList<TMLink>& nonCores ) {
  unsigned n = input.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink&            t = *input[i];
    const TMDCWireHit& h = *t.hit();
    if ( h.state() & WireHitFittingValid ) cores.append( t );
    else nonCores.append( t );
  }
}
 
AList<TMLink> Cores( const AList<TMLink>& input ) {
  AList<TMLink> a;
  unsigned      n = input.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink&            t = *input[i];
    const TMDCWireHit& h = *t.hit();
    if ( h.state() & WireHitFittingValid ) a.append( t );
  }
  return a;
}
 
#if defined( __GNUG__ )
int SortByWireId( const TMLink** a, const TMLink** b ) {
  if ( ( *a )->wire()->id() > ( *b )->wire()->id() ) return 1;
  else if ( ( *a )->wire()->id() == ( *b )->wire()->id() ) return 0;
  else return -1;
}
 
int SortByX( const TMLink** a, const TMLink** b ) {
  if ( ( *a )->position().x() > ( *b )->position().x() ) return 1;
  else if ( ( *a )->position().x() == ( *b )->position().x() ) return 0;
  else return -1;
}
 
#else
extern "C" int SortByWireId( const void* av, const void* bv ) {
  const TMLink** a( (const TMLink**)av );
  const TMLink** b( (const TMLink**)bv );
  if ( ( *a )->wire()->id() > ( *b )->wire()->id() ) return 1;
  else if ( ( *a )->wire()->id() == ( *b )->wire()->id() ) return 0;
  else return -1;
}
 
extern "C" int SortByX( const void* av, const void* bv ) {
  const TMLink** a( (const TMLink**)av );
  const TMLink** b( (const TMLink**)bv );
  if ( ( *a )->position().x() > ( *b )->position().x() ) return 1;
  else if ( ( *a )->position().x() == ( *b )->position().x() ) return 0;
  else return -1;
}
 
#endif
 
unsigned Width( const AList<TMLink>& list ) {
  unsigned n = list.length();
  if ( n < 2 ) return n;
 
  const TMDCWire* w      = list[0]->wire();
  unsigned        nWires = w->layer()->nWires();
  unsigned        center = w->localId();
 
#ifdef TRKRECO_DEBUG_DETAIL
  unsigned sId = w->superLayerId();
#endif
 
  unsigned left  = 0;
  unsigned right = 0;
  for ( unsigned i = 1; i < n; i++ )
  {
    w           = list[i]->wire();
    unsigned id = w->localId();
 
    unsigned distance0, distance1;
    if ( id > center )
    {
      distance0 = id - center;
      distance1 = nWires - distance0;
    }
    else
    {
      distance1 = center - id;
      distance0 = nWires - distance1;
    }
 
    if ( distance0 < distance1 )
    {
      if ( distance0 > right ) right = distance0;
    }
    else
    {
      if ( distance1 > left ) left = distance1;
    }
 
#ifdef TRKRECO_DEBUG_DETAIL
    if ( w->superLayerId() != sId )
      std::cout << "::width !!! super layer assumption violation" << std::endl;
#endif
  }
 
  return right + left + 1;
}
 
AList<TMLink> Edges( const AList<TMLink>& list ) {
  AList<TMLink> a;
 
  unsigned n = list.length();
  if ( n < 2 ) return a;
  else if ( n == 2 ) return list;
 
  const TMDCWire* w      = list[0]->wire();
  unsigned        nWires = w->layer()->nWires();
  unsigned        center = w->localId();
 
  unsigned left   = 0;
  unsigned right  = 0;
  TMLink*  leftL  = list[0];
  TMLink*  rightL = list[0];
  for ( unsigned i = 1; i < n; i++ )
  {
    w           = list[i]->wire();
    unsigned id = w->localId();
 
    unsigned distance0, distance1;
    if ( id > center )
    {
      distance0 = id - center;
      distance1 = nWires - distance0;
    }
    else
    {
      distance1 = center - id;
      distance0 = nWires - distance1;
    }
 
    if ( distance0 < distance1 )
    {
      if ( distance0 > right )
      {
        right  = distance0;
        rightL = list[i];
      }
    }
    else
    {
      if ( distance1 > left )
      {
        left  = distance1;
        leftL = list[i];
      }
    }
  }
 
  a.append( leftL );
  a.append( rightL );
  return a;
}
 
AList<TMLink> SameLayer( const AList<TMLink>& list, const TMLink& a ) {
  AList<TMLink> same;
  unsigned      id = a.wire()->layerId();
  unsigned      n  = list.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    if ( list[i]->wire()->layerId() == id ) same.append( list[i] );
  }
  return same;
}
 
AList<TMLink> SameSuperLayer( const AList<TMLink>& list, const TMLink& a ) {
  AList<TMLink> same;
  unsigned      id = a.wire()->superLayerId();
  unsigned      n  = list.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    if ( list[i]->wire()->superLayerId() == id ) same.append( list[i] );
  }
  return same;
}
 
AList<TMLink> SameLayer( const AList<TMLink>& list, unsigned id ) {
  AList<TMLink> same;
  unsigned      n = list.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    if ( list[i]->wire()->layerId() == id ) same.append( list[i] );
  }
  return same;
}
 
AList<TMLink> SameSuperLayer( const AList<TMLink>& list, unsigned id ) {
  AList<TMLink> same;
  unsigned      n = list.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    if ( list[i]->wire()->superLayerId() == id ) same.append( list[i] );
  }
  return same;
}
 
AList<TMLink> InOut( const AList<TMLink>& list ) {
  AList<TMLink> inners;
  AList<TMLink> outers;
  unsigned      n              = list.length();
  unsigned      innerMostLayer = 999;
  unsigned      outerMostLayer = 0;
  for ( unsigned i = 0; i < n; i++ )
  {
    unsigned id = list[i]->wire()->layerId();
    if ( id < innerMostLayer ) innerMostLayer = id;
    else if ( id > outerMostLayer ) outerMostLayer = id;
  }
  for ( unsigned i = 0; i < n; i++ )
  {
    unsigned id = list[i]->wire()->layerId();
    if ( id == innerMostLayer ) inners.append( list[i] );
    else if ( id == outerMostLayer ) outers.append( list[i] );
  }
  inners.append( outers );
  return inners;
}
 
unsigned SuperLayer( const AList<TMLink>& list ) {
  unsigned sl = 0;
  unsigned n  = list.length();
  for ( unsigned i = 0; i < n; i++ ) sl |= ( 1 << ( list[i]->wire()->superLayerId() ) );
  return sl;
}
 
unsigned SuperLayer( const AList<TMLink>& links, unsigned minN ) {
  unsigned n         = links.length();
  unsigned nHits[11] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
  for ( unsigned i = 0; i < n; i++ ) ++nHits[links[i]->wire()->superLayerId()];
  unsigned sl = 0;
  for ( unsigned i = 0; i < 11; i++ )
    if ( nHits[i] >= minN ) sl |= ( 1 << i );
  return sl;
}
 
unsigned NSuperLayers( const AList<TMLink>& list ) {
  unsigned l0 = 0;
  unsigned n  = list.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    unsigned id = list[i]->wire()->superLayerId();
    l0 |= ( 1 << id );
  }
 
  unsigned l = 0;
  for ( unsigned i = 0; i < 11; i++ )
  {
    if ( l0 & ( 1 << i ) ) ++l;
  }
  return l;
}
 
unsigned NSuperLayers( const AList<TMLink>& links, unsigned minN ) {
  unsigned n         = links.length();
  unsigned nHits[11] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
  for ( unsigned i = 0; i < n; i++ ) ++nHits[links[i]->wire()->superLayerId()];
  unsigned sl = 0;
  for ( unsigned i = 0; i < 11; i++ )
    if ( nHits[i] >= minN ) ++sl;
  return sl;
}
 
unsigned NMissingAxialSuperLayers( const AList<TMLink>& links ) {
  unsigned n = links.length();
  // Liuqg, change the following to BES.    unsigned nHits[6] = {0, 0, 0, 0, 0, 0};
  unsigned nHits[5] = { 0, 0, 0, 0, 0 };
  for ( unsigned i = 0; i < n; i++ )
    if ( links[i]->wire()->axial() ) ++nHits[links[i]->wire()->axialStereoLayerId() / 4];
  unsigned j = 0;
  while ( nHits[j] == 0 ) ++j;
  unsigned nMissing = 0;
  unsigned nMax     = 0;
  for ( unsigned i = j; i < 5; i++ )
  {
    if ( nHits[i] == 0 ) ++nMissing;
    else
    {
      if ( nMax < nMissing ) nMax = nMissing;
      nMissing = 0;
    }
  }
  return nMax;
}
 
const TTrackHEP& Links2HEP( const AList<TMLink>& links ) {
  const TTrackHEP*        best = NULL;
  const AList<TTrackHEP>& list = TTrackHEP::list();
  unsigned                nHep = list.length();
 
  if ( !nHep ) return *best;
 
  unsigned* N = (unsigned*)malloc( nHep * sizeof( unsigned ) );
  for ( unsigned i = 0; i < nHep; i++ ) N[i] = 0;
 
  for ( unsigned i = 0; i < links.length(); i++ )
  {
    const TMLink&    l   = *links[i];
    const TTrackHEP& hep = *l.hit()->mc()->hep();
    for ( unsigned j = 0; j < nHep; j++ )
      if ( list[j] == &hep ) ++N[j];
  }
 
  unsigned nMax = 0;
  for ( unsigned i = 0; i < nHep; i++ )
  {
    if ( N[i] > nMax )
    {
      best = list[i];
      nMax = N[i];
    }
  }
 
  return *best;
}
 
/*
double
TMLink::DriftTime(double _tof,double z) const {
    double tprop = 0.;
    double _vprop = (_layer<8) ? Constants::vpropInner : Constants::vpropOuter;
    if (0 == _layer%2){
      tprop = (0.5*_zlen + z)/_vprop; //odd
    }else{
      tprop = (0.5*_zlen - z)/_vprop; //even
    }
    double driftT;
    driftT = fabs(_rawTime - _T0Walk -1.e9*tof - tprop);
 
    //std::cout<< "lay "<<_layer<<" cell "<<_wire<<" zhit "<<z<<" tprop "<<tprop << std::endl;
    return driftT;
}
*/