TMLine.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TMLine.cxx,v 1.7 2012/05/28 05:16:29 maoh Exp $
//-----------------------------------------------------------------------------
// Filename : TMLine.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class to represent a line in tracking.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#include "TrkReco/TMLine.h"
#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TTrackHEP.h"
 
const TLineFitter TMLine::_fitter = TLineFitter( "TMLine Default Line Fitter" );
 
TMLine::TMLine()
    : TTrackBase()
    , _a( 0. )
    , _b( 0. )
    , _det( 0. )
    , _fittedUpdated( false )
    , _chi2( 0. )
    , _reducedChi2( 0. ) {
 
  //...Set a defualt fitter...
  fitter( &TMLine::_fitter );
}
 
TMLine::TMLine( const AList<TMLink>& a )
    : TTrackBase( a )
    , _a( 0. )
    , _b( 0. )
    , _det( 0. )
    , _fittedUpdated( false )
    , _chi2( 0. )
    , _reducedChi2( 0. ) {
 
  //...Set a defualt fitter...
  fitter( &TMLine::_fitter );
}
 
TMLine::~TMLine() {}
 
void TMLine::dump( const std::string& msg, const std::string& pre ) const {
  bool def = false;
  if ( msg == "" ) def = true;
 
  if ( def || msg.find( "line" ) != std::string::npos ||
       msg.find( "detail" ) != std::string::npos )
  {
    std::cout << pre;
    std::cout << "#links=" << _links.length();
    std::cout << ",a=" << _a;
    std::cout << ",b=" << _b;
    std::cout << ",det=" << _det;
    std::cout << std::endl;
  }
  if ( !def ) TTrackBase::dump( msg, pre );
}
 
//  int
//  TMLine::fitx(void) {
//      if (_fitted) return 0;
 
//      unsigned n = _links.length();
//      double sum = double(n);
//      double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
//      for (unsigned i = 0; i < n; i++) {
//      TMLink & l = * _links[i];
 
//      double x = l.position().x();
//      double y = l.position().y();
//      sumX  += x;
//      sumY  += y;
//      sumX2 += x * x;
//      sumXY += x * y;
//      sumY2 += y * y;
//      }
 
//      _det = sum * sumX2 - sumX * sumX;
//  #ifdef TRKRECO_DEBUG_DETAIL
//      std::cout << "    TMLine::fit ... det=" << _det << std::endl;
//  #endif
//      if(_det == 0. && n != 2){
//      return -1;
//      }else if(_det == 0. && n == 2){
//      double x0 = _links[0]->position().x();
//      double y0 = _links[0]->position().y();
//      double x1 = _links[1]->position().x();
//      double y1 = _links[1]->position().y();
//      if(x0 == x1)return -1;
//      _a = (y0-y1)/(x0-x1);
//      _b = -_a*x1 + y1;
 
//      _fitted = true;
//      return 0;
//      }
//      _a = (sumXY * sum - sumX * sumY) / _det;
//      _b = (sumX2 * sumY - sumX * sumXY) / _det;
 
//      _fitted = true;
//      return 0;
//  }
 
double TMLine::chi2( void ) const {
#ifdef TRKRECO_DEBUG
  if ( !_fitted ) std::cout << "TMLine::chi2 !!! fit not performed" << std::endl;
#endif
 
  if ( _fittedUpdated ) return _chi2;
  _chi2      = 0.;
  unsigned n = _links.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink& l = *_links[i];
 
    double x = l.position().x();
    double y = l.position().y();
    double c = y - _a * x - _b;
    _chi2 += c * c;
  }
  _fittedUpdated = true;
  return _chi2;
}
 
void TMLine::refine( AList<TMLink>& list, float maxSigma ) {
  AList<TMLink> bad;
  unsigned      n = _links.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink& l    = *_links[i];
    double  dist = distance( l );
    if ( dist > maxSigma ) bad.append( l );
  }
 
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    TMLine::refine ... rejected hits:max distance=" << maxSigma;
  std::cout << std::endl;
  bad.sort( SortByWireId );
  for ( unsigned i = 0; i < bad.length(); i++ )
  {
    TMLink& l = *_links[i];
    std::cout << "        ";
    std::cout << l.wire()->layerId() << "-";
    std::cout << l.wire()->localId();
    std::cout << "(";
    if ( l.hit()->mc() )
    {
      if ( l.hit()->mc()->hep() ) std::cout << l.hit()->mc()->hep()->id();
      else std::cout << "0";
    }
    std::cout << "),";
    std::cout << l.position() << "," << distance( l );
    if ( distance( l ) > maxSigma ) std::cout << " X";
    std::cout << std::endl;
  }
#endif
 
  if ( bad.length() )
  {
    _links.remove( bad );
    list.append( bad );
    _fitted        = false;
    _fittedUpdated = false;
  }
}
 
int TMLine::fit2() {
  std::cout << "   " << __FILE__ << "   " << __LINE__ << " ERROR : nsl" << std::endl;
  //    if (_fitted) return 0;
 
  unsigned n         = _links.length();
  int      mask[100] = { 0 };
  int      nsl[11]   = { 64, 80, 96, 128, 144, 160, 192, 208, 240, 256, 288 };
  int      npos = 0, nneg = 0;
  for ( unsigned i = 0; i < n - 1; i++ )
  {
    TMLink& l = *_links[i];
    for ( unsigned j = i + 1; j < n; j++ )
    {
      TMLink& s = *_links[j];
      if ( l.hit()->wire()->layerId() == s.hit()->wire()->layerId() )
      {
        //... Check 3 consective hits
        if ( i > 0 && ( mask[i - 1] == 1 && mask[j] == 1 ) )
        {
          TMLink& t = *_links[i - 1];
          if ( l.hit()->wire()->layerId() == t.hit()->wire()->layerId() ) { mask[i] = 1; }
        }
        int ilast  = nsl[l.hit()->wire()->superLayerId()] - 1;
        int ilocal = l.hit()->wire()->localId();
        int jlocal = s.hit()->wire()->localId();
        if ( ilocal > 0 && ilocal < ilast )
        {
          if ( abs( jlocal - ilocal ) > 1 )
          {
            mask[i] = 1;
            mask[j] = 1;
          }
        }
        else if ( ilocal == 0 )
        {
          if ( jlocal > 1 && jlocal < ilast )
          {
            mask[i] = 1;
            mask[j] = 1;
          }
        }
        else if ( ilocal == ilast )
        {
          if ( jlocal > 0 && jlocal < ilast - 1 )
          {
            mask[i] = 1;
            mask[j] = 1;
          }
        }
      }
    }
    //...
    if ( mask[i] == 0 )
    {
      if ( l.position().y() >= 0 ) npos += 1;
      if ( l.position().y() < 0 ) nneg += 1;
    }
  }
 
  //....
  double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
  int    nused = 0;
  int    lyid[2];
  for ( unsigned i = 0; i < n; i++ )
  {
 
    if ( mask[i] == 1 ) continue;
 
    TMLink& l = *_links[i];
 
    double x = l.position().x();
    double y = l.position().y();
    if ( abs( npos - nneg ) > 3 )
    {
      if ( npos > nneg && y < 0 ) continue;
      if ( npos < nneg && y > 0 ) continue;
    }
    sumX += x;
    sumY += y;
    sumX2 += x * x;
    sumXY += x * y;
    sumY2 += y * y;
    if ( nused < 2 ) { lyid[nused] = l.hit()->wire()->layerId(); }
    nused += 1;
  }
 
  if ( nused < 2 || ( nused == 2 && lyid[0] == lyid[1] ) ) { return -2; }
  double sum = double( nused );
  _det       = sum * sumX2 - sumX * sumX;
  if ( _det == 0. ) { return -1; }
  _a = ( sumXY * sum - sumX * sumY ) / _det;
  _b = ( sumX2 * sumY - sumX * sumXY ) / _det;
 
  _fitted = true;
  return 0;
}
 
int TMLine::fit2s() {
  //    if (_fitted) return 0;
 
  unsigned n         = _links.length();
  int      mask[100] = { 0 };
  int      npos = 0, nneg = 0;
  for ( unsigned i = 0; i < n - 1; i++ )
  {
    TMLink& l = *_links[i];
    for ( unsigned j = i + 1; j < n; j++ )
    {
      TMLink& s = *_links[j];
      if ( l.hit()->wire()->layerId() == s.hit()->wire()->layerId() )
      {
        mask[i] = 1;
        mask[j] = 1;
      }
    }
    //...
    if ( mask[i] == 0 )
    {
      if ( l.position().y() >= 0 ) npos += 1;
      if ( l.position().y() < 0 ) nneg += 1;
    }
  }
 
  //....
  double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
  int    nused = 0;
  int    lyid[2];
  for ( unsigned i = 0; i < n; i++ )
  {
 
    if ( mask[i] == 1 ) continue;
 
    TMLink& l = *_links[i];
 
    double x = l.position().x();
    double y = l.position().y();
    if ( npos > nneg && y < 0 ) continue;
    if ( npos < nneg && y > 0 ) continue;
 
    sumX += x;
    sumY += y;
    sumX2 += x * x;
    sumXY += x * y;
    sumY2 += y * y;
    nused += 1;
  }
 
  if ( nused < 4 ) { return -2; }
  double sum = double( nused );
  _det       = sum * sumX2 - sumX * sumX;
  if ( _det == 0. ) { return -1; }
  _a = ( sumXY * sum - sumX * sumY ) / _det;
  _b = ( sumX2 * sumY - sumX * sumXY ) / _det;
 
  _fitted = true;
  return 0;
}
int TMLine::fit2sp() {
  //    if (_fitted) return 0;
 
  unsigned n         = _links.length();
  int      mask[100] = { 0 };
  double   phi_ave   = 0.;
  int      nphi      = 0;
  double   Crad      = 180. / 3.141592;
  for ( unsigned i = 0; i < n - 1; i++ )
  {
    TMLink& l = *_links[i];
    for ( unsigned j = i + 1; j < n; j++ )
    {
      TMLink& s = *_links[j];
      if ( l.hit()->wire()->layerId() == s.hit()->wire()->layerId() )
      {
        mask[i] = 1;
        mask[j] = 1;
      }
    }
    //...
    if ( mask[i] != 1 )
    {
      double phi = Crad * atan2( l.position().y(), l.position().x() );
      phi_ave += phi;
      nphi += 1;
    }
  }
 
  //...
  if ( mask[n - 1] != 1 )
  {
    TMLink& l   = *_links[n - 1];
    double  phi = Crad * atan2( l.position().y(), l.position().x() );
    phi_ave += phi;
    nphi += 1;
  }
  double phi_max = 0.;
  double phi_min = 0.;
  if ( nphi > 0 )
  {
    phi_max = phi_ave / n + 40;
    phi_min = phi_ave / n - 40;
  }
 
  //....
  double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
  int    nused = 0;
  int    lyid[2];
  for ( unsigned i = 0; i < n; i++ )
  {
 
    if ( mask[i] == 1 ) continue;
 
    TMLink& l = *_links[i];
 
    double x   = l.position().x();
    double y   = l.position().y();
    double phi = Crad * atan2( l.position().y(), l.position().x() );
    if ( phi > phi_max && phi < phi_min ) continue;
 
    sumX += x;
    sumY += y;
    sumX2 += x * x;
    sumXY += x * y;
    sumY2 += y * y;
    nused += 1;
  }
 
  if ( nused < 4 ) { return -2; }
  double sum = double( nused );
  _det       = sum * sumX2 - sumX * sumX;
  if ( _det == 0. ) { return -1; }
  _a = ( sumXY * sum - sumX * sumY ) / _det;
  _b = ( sumX2 * sumY - sumX * sumXY ) / _det;
 
  _fitted = true;
  return 0;
}
 
int TMLine::fit2p() {
  //    if (_fitted) return 0;
 
  std::cout << "   " << __FILE__ << "   " << __LINE__ << " ERROR : nsl" << std::endl;
  unsigned n         = _links.length();
  int      mask[100] = { 0 };
  int      nsl[11]   = { 64, 80, 96, 128, 144, 160, 192, 208, 240, 256, 288 };
  double   phi_ave   = 0.;
  int      nphi      = 0;
  double   Crad      = 180. / 3.141592;
  for ( unsigned i = 0; i < n - 1; i++ )
  {
    TMLink& l = *_links[i];
    for ( unsigned j = i + 1; j < n; j++ )
    {
      TMLink& s = *_links[j];
      if ( l.hit()->wire()->layerId() == s.hit()->wire()->layerId() )
      {
        //... Check 3 consective hits
        if ( i > 0 && ( mask[i - 1] == 1 && mask[j] == 1 ) )
        {
          TMLink& t = *_links[i - 1];
          if ( l.hit()->wire()->layerId() == t.hit()->wire()->layerId() ) { mask[i] = 1; }
        }
        int ilast  = nsl[l.hit()->wire()->superLayerId()] - 1;
        int ilocal = l.hit()->wire()->localId();
        int jlocal = s.hit()->wire()->localId();
        if ( ilocal > 0 && ilocal < ilast )
        {
          if ( abs( jlocal - ilocal ) > 1 )
          {
            mask[i] = 1;
            mask[j] = 1;
          }
        }
        else if ( ilocal == 0 )
        {
          if ( jlocal > 1 && jlocal < ilast )
          {
            mask[i] = 1;
            mask[j] = 1;
          }
        }
        else if ( ilocal == ilast )
        {
          if ( jlocal > 0 && jlocal < ilast - 1 )
          {
            mask[i] = 1;
            mask[j] = 1;
          }
        }
      }
    }
    //...
    //...
    if ( mask[i] != 1 )
    {
      double phi = Crad * atan2( l.position().y(), l.position().x() );
      phi_ave += phi;
      nphi += 1;
    }
  }
 
  //...
  if ( mask[n - 1] != 1 )
  {
    TMLink& l   = *_links[n - 1];
    double  phi = Crad * atan2( l.position().y(), l.position().x() );
    phi_ave += phi;
    nphi += 1;
  }
  double phi_max = 0.;
  double phi_min = 0.;
  if ( nphi > 0 )
  {
    phi_max = phi_ave / n + 40;
    phi_min = phi_ave / n - 40;
  }
 
  //....
  double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
  int    nused = 0;
  int    lyid[2];
  for ( unsigned i = 0; i < n; i++ )
  {
 
    if ( mask[i] == 1 ) continue;
 
    TMLink& l = *_links[i];
 
    double x   = l.position().x();
    double y   = l.position().y();
    double phi = Crad * atan2( l.position().y(), l.position().x() );
    if ( phi > phi_max && phi < phi_min ) continue;
 
    sumX += x;
    sumY += y;
    sumX2 += x * x;
    sumXY += x * y;
    sumY2 += y * y;
    if ( nused < 2 ) { lyid[nused] = l.hit()->wire()->layerId(); }
    nused += 1;
  }
 
  if ( nused < 2 || ( nused == 2 && lyid[0] == lyid[1] ) ) { return -2; }
  double sum = double( nused );
  _det       = sum * sumX2 - sumX * sumX;
  if ( _det == 0. ) { return -1; }
  _a = ( sumXY * sum - sumX * sumY ) / _det;
  _b = ( sumX2 * sumY - sumX * sumXY ) / _det;
 
  _fitted = true;
  return 0;
}
 
void TMLine::removeChits() {
 
  unsigned n        = _links.length();
  int      nlyr[50] = { 0 };
  int      nneg = 0, npos = 0;
  for ( unsigned i = 0; i < n - 1; i++ )
  {
    TMLink& l = *_links[i];
    nlyr[l.hit()->wire()->layerId()] += 1;
    if ( l.position().y() < 0. ) { nneg += 1; }
    else { npos += 1; }
  }
 
  //...
  AList<TMLink> bad;
  for ( unsigned i = 0; i < n; i++ )
  {
 
    TMLink& l = *_links[i];
 
    //...if # of hits in a wire layer, don't use...
    if ( nlyr[l.hit()->wire()->layerId()] > 3 )
    {
      bad.append( l );
      continue;
    }
    //...remove extremely bad poinits
    if ( abs( nneg - npos ) > 3 )
    {
      if ( npos > nneg && l.position().y() < 0 ) bad.append( l );
      if ( npos < nneg && l.position().y() > 0 ) bad.append( l );
    }
  }
  //...
  if ( bad.length() > 0 && bad.length() < n ) { _links.remove( bad ); }
 
  //... For the next fit
  _fitted        = false;
  _fittedUpdated = false;
}
 
void TMLine::removeSLY( AList<TMLink>& list ) {
  _links.remove( list );
  _fitted        = false;
  _fittedUpdated = false;
}
 
void TMLine::appendSLY( AList<TMLink>& list ) {
  _links.append( list );
  _fitted        = false;
  _fittedUpdated = false;
}
 
void TMLine::appendByszdistance( AList<TMLink>& list, unsigned isl, float maxSigma ) {
 
  //... intialize
  unsigned nb = _links.length();
 
  //....Select good hit
  unsigned n = list.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink& l = *list[i];
    if ( l.hit()->wire()->superLayerId() == isl )
    {
      double dist = distance( l );
      if ( dist < maxSigma ) { _links.append( l ); }
    }
  }
 
  unsigned na = _links.length();
  if ( nb != na )
  {
    AList<TMLink> bad;
    //... remove duplicated hits
    for ( unsigned i = 0; i < na; i++ )
    {
      TMLink& l = *_links[i];
      if ( i < na - 1 )
      {
        TMLink& lnext = *_links[i + 1];
        if ( l.hit()->wire()->layerId() == lnext.hit()->wire()->layerId() )
        {
          if ( l.hit()->wire()->localId() == lnext.hit()->wire()->localId() )
          { bad.append( l ); }
        }
      }
    }
    if ( bad.length() > 0 ) _links.remove( bad );
    _fitted        = false;
    _fittedUpdated = false;
  }
}
 
double TMLine::reducedChi2( void ) const {
#ifdef TRKRECO_DEBUG
  if ( !_fitted ) std::cout << "TMLine::reducedChi2 !!! fit not performed" << std::endl;
#endif
 
  if ( _fittedUpdated ) return _reducedChi2;
  double   chi2  = 0.;
  double   scale = 20.;
  unsigned n     = _links.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink& l = *_links[i];
 
    double x   = l.position().x();
    double y   = l.position().y();
    double c   = y - _a * x - _b;
    double err = 1.;
    if ( l.hit() ) err = scale * l.hit()->dDrift();
    chi2 += c * c / err / err;
  }
 
  _reducedChi2   = chi2 / ( n - 2 );
  _fittedUpdated = true;
  return _reducedChi2;
}
 
#if defined( __GNUG__ )
 
int SortByB( const TMLine** a, const TMLine** b ) {
  if ( fabs( ( *a )->b() ) > fabs( ( *b )->b() ) ) return 1;
  else if ( fabs( ( *a )->b() ) == fabs( ( *b )->b() ) ) return 0;
  else return -1;
}
 
#else
 
extern "C" int SortByB( const void* av, const void* bv ) {
  const TMLine** a( (const TMLine**)av );
  const TMLine** b( (const TMLine**)bv );
  if ( fabs( ( *a )->b() ) > fabs( ( *b )->b() ) ) return 1;
  else if ( fabs( ( *a )->b() ) == fabs( ( *b )->b() ) ) return 0;
  else return -1;
}
 
#endif