TCircle.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TCircle.cxx,v 1.10 2010/03/31 09:58:59 liucy Exp $
//-----------------------------------------------------------------------------
// Filename : TCircle.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class to represent a circle in tracking.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#include "TrkReco/TCircle.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TMLink.h"
 
const TCircleFitter TCircle::_fitter = TCircleFitter( "TCircle Default Circle Fitter" );
 
TCircle::TCircle( const AList<TMLink>& a ) : TTrackBase( a ), _radius( 0. ), _charge( 0. ) {
 
  //...Set a defualt fitter...
  fitter( &TCircle::_fitter );
}
 
TCircle::~TCircle() {}
 
void TCircle::dump( const std::string& msg, const std::string& pre ) const {
  bool def = false;
  if ( msg == "" ) def = true;
 
  if ( def || msg.find( "circle" ) != std::string::npos ||
       msg.find( "detail" ) != std::string::npos )
  {
    std::cout << pre;
    std::cout << "#links=" << _links.length();
    if ( _fitted )
    {
      std::cout << ",charge=" << _charge;
      std::cout << ",center=" << _center;
      std::cout << ",radius=" << _radius;
      std::cout << std::endl << pre;
      std::cout << "pt=" << pt();
      std::cout << ",impact=" << impact();
      std::cout << std::endl;
    }
    else { std::cout << ", not fitted yet" << std::endl; }
  }
  if ( !def ) TTrackBase::dump( msg, pre );
}
 
//  int
//  TCircle::fitx(void) {
//  #ifdef TRKRECO_DEBUG_DETAIL
//      std::cout << "    TCircle::fit ..." << std::endl;
//  #endif
//      if (_fitted) return 0;
//      unsigned n = _links.length();
//      if (n < 3) return -1;
 
//      for (unsigned i = 0; i < n; i++) {
//      TMLink * l = _links[i];
//      const TMDCWireHit * h = l->hit();
 
//      //...Check next hit...
//      Point3D point;
//      if (h->state() & WireHitPatternLeft)
//          point = h->position(WireHitLeft);
//      else if (h->state() & WireHitPatternRight)
//          point = h->position(WireHitRight);
//      else
//          point = h->xyPosition();
//      // float weight = 1. / (h->distance() * h->distance());
//      // float weight = 1. / h->distance();
 
//      _circle.add_point(point.x(), point.y()); //, weight);
//      }
 
//      if(_circle.fit()<0.0 || _circle.kappa()==0.0) return -1;
//      HepVector v(_circle.center());
//      _center.setX(v(1));
//      _center.setY(v(2));
//      _radius = _circle.radius();
 
//      //...Determine charge...Better way???
//      int qSum = 0;
//      for (unsigned i = 0; i < n; i++) {
//      TMLink * l = _links[i];
//      if (l == 0) continue;
 
//      const TMDCWireHit * h = l->hit();
//      if (h == 0) continue;
 
//      float q = (_center.cross(h->xyPosition())).z();
//      if (q > 0.) qSum += 1;
//      else        qSum -= 1;
//      }
//      if (qSum >= 0) _charge = +1.;
//      else           _charge = -1.;
//      _radius *= _charge;
 
//      _fitted = true;
//      return 0;
//  }
 
double TCircle::weight( const TMLink& l ) const {
 
  //...Axial Wires
  int maxLink = 0;
  int localID[6];
  int layerID[6];
  int LayerID = l.hit()->wire()->layerId();
  for ( int i = 0; i < 6; i++ )
  {
    if ( l.neighbor( i ) )
    {
      maxLink    = i;
      localID[i] = l.neighbor( i )->hit()->wire()->localId();
      layerID[i] = l.neighbor( i )->hit()->wire()->layerId();
    }
    else { break; }
  }
  if ( maxLink != 1 ) return 1.0;
  if ( layerID[0] == LayerID && layerID[1] == LayerID ) return 0.5;
  if ( layerID[0] + layerID[1] != LayerID * 2 ) return 1.0;
  // Liuqg  if(localID[0] != localID[1])return 1.5;//1.0 or 2.0 ??
  return 1.0;
}
 
int TCircle::fitForCurl( int ipConst ) {
  unsigned n = _links.length();
  if ( n < 3 ) return -1;
 
  // IP check
  unsigned flagIP  = 1;
  unsigned layerID = _links[0]->hit()->wire()->layerId();
  for ( unsigned i = 0; i < n; i++ )
  {
    if ( layerID != _links[i]->hit()->wire()->layerId() )
    {
      flagIP = 0;
      break;
    }
  }
  if ( ipConst != 0 ) flagIP = 1;
  if ( flagIP == 1 )
  {
    _circle.add_point( 0., 0., 0.5 );
    //++_nPointsForFit;
  }
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink* l = _links[i];
    if ( l == 0 ) continue;
 
    const TMDCWireHit* h = l->hit();
    if ( h == 0 ) continue;
 
    //...Check next hit...
    HepPoint3D point;
    point = h->xyPosition();
 
    double weight = 1.0;
    weight        = this->weight( *l );
 
    _circle.add_point( point.x(), point.y(), weight );
    //++_nPointsForFit;
  }
 
  // if (_nPointsForFit < 3) return -1;
  if ( _circle.fit() < 0.0 || _circle.kappa() == 0.0 ) return -1;
  HepVector v( _circle.center() );
  _center.setX( v( 1 ) );
  _center.setY( v( 2 ) );
  _radius = _circle.radius();
 
  //...Determine charge...Better way???
  int qSum = 0;
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink* l = _links[i];
    if ( l == 0 ) continue;
 
    const TMDCWireHit* h = l->hit();
    if ( h == 0 ) continue;
 
    float q = ( _center.cross( h->xyPosition() ) ).z();
    if ( q > 0. ) qSum += 1;
    else qSum -= 1;
  }
  if ( qSum >= 0 ) _charge = +1.;
  else _charge = -1.;
  _radius *= _charge;
 
  _fitted = true;
  return 0;
}