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    : [email protected]
//-----------------------------------------------------------------------------
// 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;
}