TCircleFitter.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TCircleFitter.cxx,v 1.8 2010/03/31 09:58:59 liucy Exp $
//-----------------------------------------------------------------------------
// Filename : TCircleFitter.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class to fit a TTrackBase object to a circle.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#include "TrkReco/TCircleFitter.h"
#include "TrkReco/TCircle.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TTrackBase.h"
// #include "lpav/Lpav.h"
// #include "TrkReco/Lpav.h"
#include "TrackUtil/Lpav.h"
 
TCircleFitter::TCircleFitter( const std::string& name )
    : TMFitter( name ), _charge( 0. ), _radius( 0. ), _center( HepPoint3D( 0., 0., 0. ) ) {}
 
TCircleFitter::~TCircleFitter() {}
 
int TCircleFitter::fit( TTrackBase& t ) const {
 
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    TCircleFitter::fit ..." << std::endl;
#endif
 
  //...Already fitted ?...
  if ( t.fitted() ) return TFitAlreadyFitted;
 
  //...Check # of hits...
  if ( t.links().length() < 3 ) return TFitErrorFewHits;
 
  //...Hit loop...
  Lpav circle;
  /*
      circle.add_point(-173.011, 94.3169);
      circle.add_point(-187.6, 101.528);
      circle.add_point(-202.018, 108.601);
      circle.add_point(-216.639, 115.691);
      circle.add_point(-231.152, 122.49);
      circle.add_point(-246.207, 129.496);
      circle.add_point(-261.179, 136.407);
      circle.add_point(-275.198, 142.764);
      circle.add_point(-290.519, 149.572);
      circle.add_point(-304.42, 155.7);
      circle.add_point(-319.947, 162.285);
      circle.add_point(-335.495, 168.731);
      circle.add_point(-608.394, 264.577);
      circle.add_point(-612.842, 265.869);
      circle.add_point(-627.175, 270.011);
      circle.add_point(-643.654, 274.719);
      circle.add_point(-657.809, 278.683);
      circle.add_point(-674.464, 283.209);
      circle.add_point(-690.7, 287.541);
      circle.add_point(-704.925, 291.091);
    */
  unsigned n = t.links().length();
 
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink*            l = t.links()[i];
    const TMDCWireHit* h = l->hit();
 
    //...Check next hit...
    HepPoint3D 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();
 
    //  if (h->state() & WireHitPatternLeft) cout<<"h: "<<h->wire()->layerId()<<" local: "
    //                      <<h->wire()->localId()<<" TCirclefitter:
    // L"<<endl;    else cout<<"h: "<<h->wire()->layerId()<<" local: "
    //      <<h->wire()->localId()<<" TCirclefitter: R"<<endl;
 
    circle.add_point( point.x(), point.y() ); //, weight);
                                              // yuany
 
    //  cout<<"TCirclefitter("<<sqrt(point.x()*point.x()+point.y()*point.y())<<","
    //  <<point.y()/point.x()<<") dis:"
    //  <<(h->drift(0))*10000/40<<" point "<<point<<" L "<<(h->state()&WireHitPatternLeft)
    //  <<" R "<<(h->state()&WireHitPatternRight)<<endl;
 
    //  cout<<"TCirclefitter "<<point.x()<<"    "<<point.y()<<"    L
    //"<<(h->state()&WireHitPatternLeft)<<" R "<<(h->state()&WireHitPatternRight)<<endl;
  }
 
  if ( circle.fit() < 0.0 || circle.kappa() == 0.0 ) return TFitFailed;
  // yuany
  //    cout<<"circle chi2:"<<circle.fit()<<" kappa:"<<circle.kappa()<<"
  // center:"<<circle.center()<<" radius:"<<circle.radius()<<endl;
 
  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 = t.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;
 
  if ( t.objectType() == Circle ) ( (TCircle&)t ).property( _charge, _radius, _center );
  fitDone( t );
  return 0;
}