/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Reconstruction/TrkReco/TrkReco-00-08-59-patch4-slc6tag/src/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/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;    
}

---