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    : [email protected]
//-----------------------------------------------------------------------------
// Description : A class to fit a TTrackBase object to a circle.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#include "TrkReco/TCircleFitter.h"
#include "TrkReco/TTrackBase.h"
#include "TrkReco/TCircle.h"
#include "TrkReco/TMLink.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;    
}