/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Reconstruction/TrkReco/TrkReco-00-08-59-patch4-slc6tag/src/TMLine.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TMLine.cxx,v 1.7 2012/05/28 05:16:29 maoh Exp $
//-----------------------------------------------------------------------------
// Filename : TMLine.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class to represent a line in tracking.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------

#include "TrkReco/TMLine.h"
#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TTrackHEP.h"

const TLineFitter
TMLine::_fitter = TLineFitter("TMLine Default Line Fitter");

TMLine::TMLine() 
: TTrackBase(),
  _a(0.),
  _b(0.),
  _det(0.),
  _fittedUpdated(false),
  _chi2(0.),
  _reducedChi2(0.) {

    //...Set a defualt fitter...
    fitter(& TMLine::_fitter);
}

TMLine::TMLine(const AList<TMLink> & a)
: TTrackBase(a),
  _a(0.),
  _b(0.),
  _det(0.),
  _fittedUpdated(false),
  _chi2(0.),
  _reducedChi2(0.) {

    //...Set a defualt fitter...
    fitter(& TMLine::_fitter);
}

TMLine::~TMLine() {
}

void
TMLine::dump(const std::string & msg, const std::string & pre) const {
    bool def = false;
    if (msg == "") def = true;

    if (def || msg.find("line") != std::string::npos || msg.find("detail") != std::string::npos) {
        std::cout << pre;
        std::cout << "#links=" << _links.length();
        std::cout << ",a=" << _a;
        std::cout << ",b=" << _b;
        std::cout << ",det=" << _det;
        std::cout << std::endl;
    }
    if (! def) TTrackBase::dump(msg, pre);
}

//  int
//  TMLine::fitx(void) {
//      if (_fitted) return 0;

//      unsigned n = _links.length();
//      double sum = double(n);
//      double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
//      for (unsigned i = 0; i < n; i++) {
//      TMLink & l = * _links[i];

//      double x = l.position().x();
//      double y = l.position().y();
//      sumX  += x;
//      sumY  += y;
//      sumX2 += x * x;
//      sumXY += x * y;
//      sumY2 += y * y;
//      }

//      _det = sum * sumX2 - sumX * sumX;
//  #ifdef TRKRECO_DEBUG_DETAIL
//      std::cout << "    TMLine::fit ... det=" << _det << std::endl;
//  #endif
//      if(_det == 0. && n != 2){
//      return -1;
//      }else if(_det == 0. && n == 2){
//      double x0 = _links[0]->position().x();
//      double y0 = _links[0]->position().y();
//      double x1 = _links[1]->position().x();
//      double y1 = _links[1]->position().y();
//      if(x0 == x1)return -1;
//      _a = (y0-y1)/(x0-x1);
//      _b = -_a*x1 + y1;
        
//      _fitted = true;
//      return 0;
//      }
//      _a = (sumXY * sum - sumX * sumY) / _det;
//      _b = (sumX2 * sumY - sumX * sumXY) / _det;    

//      _fitted = true;
//      return 0;
//  }

double
TMLine::chi2(void) const {
#ifdef TRKRECO_DEBUG
    if (! _fitted)
        std::cout << "TMLine::chi2 !!! fit not performed" << std::endl;
#endif

    if (_fittedUpdated) return _chi2;
    _chi2 = 0.;
    unsigned n = _links.length();
    for (unsigned i = 0; i < n; i++) {
        TMLink & l = * _links[i];

        double x = l.position().x();
        double y = l.position().y();
        double c = y - _a * x - _b;
        _chi2 += c * c;
    }
    _fittedUpdated = true;
    return _chi2;
}

void
TMLine::refine(AList<TMLink> & list, float maxSigma) {
    AList<TMLink> bad;
    unsigned n = _links.length();
    for (unsigned i = 0; i < n; i++) {
        TMLink & l = * _links[i];
        double dist = distance(l);
        if (dist > maxSigma) bad.append(l);
    }

#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    TMLine::refine ... rejected hits:max distance=" << maxSigma;
    std::cout << std::endl;
    bad.sort(SortByWireId);
    for (unsigned i = 0; i < bad.length(); i++) {
        TMLink & l = * _links[i];
        std::cout << "        ";
        std::cout << l.wire()->layerId() << "-";
        std::cout << l.wire()->localId();
        std::cout << "(";
        if (l.hit()->mc()) {
            if (l.hit()->mc()->hep()) std::cout << l.hit()->mc()->hep()->id();
            else std::cout << "0";
        }
        std::cout << "),";
        std::cout << l.position() << "," << distance(l);
        if (distance(l) > maxSigma) std::cout << " X";
        std::cout << std::endl;
    }
#endif

    if (bad.length()) {
        _links.remove(bad);
        list.append(bad);
        _fitted = false;
        _fittedUpdated = false;
    }
}

int
TMLine::fit2() {
  std::cout<<"   "<<__FILE__<<"   "<<__LINE__<<" ERROR : nsl"<<std::endl;
  //    if (_fitted) return 0;

    unsigned n = _links.length();
    int mask[100] = {0};
    int nsl[11] = {64,80,96,128,144,160,192,208,240,256,288};
    int npos = 0, nneg = 0;
    for (unsigned i = 0; i < n -1 ; i++) {
        TMLink & l = * _links[i];
        for (unsigned j = i+1; j < n  ; j++) {
           TMLink & s = * _links[j];
           if(l.hit()->wire()->layerId() == s.hit()->wire()->layerId()){
             //... Check 3 consective hits
             if(i > 0 && (mask[i-1] == 1 && mask[j] == 1) ){
               TMLink & t = * _links[i-1];
               if(l.hit()->wire()->layerId() == t.hit()->wire()->layerId()){
                 mask[i] = 1;
               }
             }
             int ilast = nsl[l.hit()->wire()->superLayerId()]-1;
             int ilocal = l.hit()->wire()->localId();
             int jlocal = s.hit()->wire()->localId();
             if(ilocal > 0 && ilocal < ilast){
               if(abs(jlocal-ilocal) > 1 ) {
                  mask[i] = 1;
                  mask[j] = 1;
               }
             }else if(ilocal == 0){
               if(jlocal > 1 && jlocal < ilast){
                 mask[i] = 1;
                 mask[j] = 1;
               }
             }else if(ilocal == ilast){
               if(jlocal > 0 && jlocal < ilast-1){
                 mask[i] = 1;
                 mask[j] = 1;
               }
             }
           }
        }
        //...
        if(mask[i] == 0){
          if(l.position().y() >= 0) npos += 1;
          if(l.position().y() <  0) nneg += 1;
        }
    }            

    //....
    double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
    int nused = 0;
    int lyid[2];
    for (unsigned i = 0; i < n; i++) {

        if(mask[i] == 1) continue;

        TMLink & l = * _links[i];

        double x = l.position().x();
        double y = l.position().y();
        if(abs(npos-nneg) > 3){
          if(npos > nneg && y < 0 ) continue;
          if(npos < nneg && y > 0 ) continue;
        }
        sumX  += x;
        sumY  += y;
        sumX2 += x * x;
        sumXY += x * y;
        sumY2 += y * y;
        if(nused < 2){
          lyid[nused] = l.hit()->wire()->layerId();
        }
        nused += 1;
    }

     if(nused < 2 || (nused == 2 && lyid[0]==lyid[1])) {
       return -2;
     }
     double sum = double(nused);
    _det = sum * sumX2 - sumX * sumX;
    if(_det == 0.) {
        return -1;
    }
    _a = (sumXY * sum - sumX * sumY) / _det;
    _b = (sumX2 * sumY - sumX * sumXY) / _det;    

    _fitted = true;
    return 0;
}

int
TMLine::fit2s() {
  //    if (_fitted) return 0;

    unsigned n = _links.length();
    int mask[100] = {0};
    int npos = 0, nneg = 0;
    for (unsigned i = 0; i < n -1 ; i++) {
        TMLink & l = * _links[i];
        for (unsigned j = i+1; j < n  ; j++) {
           TMLink & s = * _links[j];
           if(l.hit()->wire()->layerId() == s.hit()->wire()->layerId()){
              mask[i] = 1;
              mask[j] = 1;
           }
        }
        //...
        if(mask[i] == 0){
          if(l.position().y() >= 0) npos += 1;
          if(l.position().y() <  0) nneg += 1;
        }
    }            

    //....
    double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
    int nused = 0;
    int lyid[2];
    for (unsigned i = 0; i < n; i++) {

        if(mask[i] == 1) continue;

        TMLink & l = * _links[i];

        double x = l.position().x();
        double y = l.position().y();
        if(npos > nneg && y < 0 ) continue;
        if(npos < nneg && y > 0 ) continue;

        sumX  += x;
        sumY  += y;
        sumX2 += x * x;
        sumXY += x * y;
        sumY2 += y * y;
        nused += 1;
    }

     if(nused < 4) {
       return -2;
     }
     double sum = double(nused);
    _det = sum * sumX2 - sumX * sumX;
    if(_det == 0.) {
        return -1;
    }
    _a = (sumXY * sum - sumX * sumY) / _det;
    _b = (sumX2 * sumY - sumX * sumXY) / _det;    

    _fitted = true;
    return 0;
}
int
TMLine::fit2sp() {
  //    if (_fitted) return 0;

    unsigned n = _links.length();
    int mask[100] = {0};
    double phi_ave = 0.;
    int nphi = 0;
    double Crad = 180./3.141592;
    for (unsigned i = 0; i < n -1 ; i++) {
        TMLink & l = * _links[i];
        for (unsigned j = i+1; j < n  ; j++) {
           TMLink & s = * _links[j];
           if(l.hit()->wire()->layerId() == s.hit()->wire()->layerId()){
              mask[i] = 1;
              mask[j] = 1;
           }
        }
        //...
        if(mask[i] != 1){
          double phi = Crad*atan2(l.position().y(), l.position().x());
          phi_ave += phi;
          nphi += 1;
        }
    }            

    //...
    if(mask[n-1] != 1){
      TMLink & l = * _links[n-1];
      double phi = Crad*atan2(l.position().y(), l.position().x());
      phi_ave += phi;
      nphi += 1;
    }
    double phi_max = 0.;
    double phi_min = 0.;
    if(nphi> 0){
      phi_max = phi_ave/n + 40;
      phi_min = phi_ave/n - 40;
    }

    //....
    double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
    int nused = 0;
    int lyid[2];
    for (unsigned i = 0; i < n; i++) {

        if(mask[i] == 1) continue;

        TMLink & l = * _links[i];

        double x = l.position().x();
        double y = l.position().y();
        double phi = Crad*atan2(l.position().y(), l.position().x());
        if(phi > phi_max && phi<phi_min ) continue;

        sumX  += x;
        sumY  += y;
        sumX2 += x * x;
        sumXY += x * y;
        sumY2 += y * y;
        nused += 1;
    }

     if(nused < 4) {
       return -2;
     }
     double sum = double(nused);
    _det = sum * sumX2 - sumX * sumX;
    if(_det == 0.) {
        return -1;
    }
    _a = (sumXY * sum - sumX * sumY) / _det;
    _b = (sumX2 * sumY - sumX * sumXY) / _det;    

    _fitted = true;
    return 0;
}

int
TMLine::fit2p() {
  //    if (_fitted) return 0;

  std::cout<<"   "<<__FILE__<<"   "<<__LINE__<<" ERROR : nsl"<<std::endl;
    unsigned n = _links.length();
    int mask[100] = {0};
    int nsl[11] = {64,80,96,128,144,160,192,208,240,256,288};
    double phi_ave = 0.;
    int nphi = 0;
    double Crad = 180./3.141592;
    for (unsigned i = 0; i < n -1 ; i++) {
        TMLink & l = * _links[i];
        for (unsigned j = i+1; j < n  ; j++) {
           TMLink & s = * _links[j];
           if(l.hit()->wire()->layerId() == s.hit()->wire()->layerId()){
             //... Check 3 consective hits
             if(i > 0 && (mask[i-1] == 1 && mask[j] == 1) ){
               TMLink & t = * _links[i-1];
               if(l.hit()->wire()->layerId() == t.hit()->wire()->layerId()){
                 mask[i] = 1;
               }
             }
             int ilast = nsl[l.hit()->wire()->superLayerId()]-1;
             int ilocal = l.hit()->wire()->localId();
             int jlocal = s.hit()->wire()->localId();
             if(ilocal > 0 && ilocal < ilast){
               if(abs(jlocal-ilocal) > 1 ) {
                  mask[i] = 1;
                  mask[j] = 1;
               }
             }else if(ilocal == 0){
               if(jlocal > 1 && jlocal < ilast){
                 mask[i] = 1;
                 mask[j] = 1;
               }
             }else if(ilocal == ilast){
               if(jlocal > 0 && jlocal < ilast-1){
                 mask[i] = 1;
                 mask[j] = 1;
               }
             }
           }
        }
        //...
        //...
        if(mask[i] != 1){
          double phi = Crad*atan2(l.position().y(), l.position().x());
          phi_ave += phi;
          nphi += 1;
        }
    }            

    //...
    if(mask[n-1] != 1){
      TMLink & l = * _links[n-1];
      double phi = Crad*atan2(l.position().y(), l.position().x());
      phi_ave += phi;
      nphi += 1;
    }
    double phi_max = 0.;
    double phi_min = 0.;
    if(nphi> 0){
      phi_max = phi_ave/n + 40;
      phi_min = phi_ave/n - 40;
    }

    //....
    double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
    int nused = 0;
    int lyid[2];
    for (unsigned i = 0; i < n; i++) {

        if(mask[i] == 1) continue;

        TMLink & l = * _links[i];

        double x = l.position().x();
        double y = l.position().y();
        double phi = Crad*atan2(l.position().y(), l.position().x());
        if(phi > phi_max && phi<phi_min ) continue;

        sumX  += x;
        sumY  += y;
        sumX2 += x * x;
        sumXY += x * y;
        sumY2 += y * y;
        if(nused < 2){
          lyid[nused] = l.hit()->wire()->layerId();
        }
        nused += 1;
    }

     if(nused < 2 || (nused == 2 && lyid[0]==lyid[1])) {
       return -2;
     }
     double sum = double(nused);
    _det = sum * sumX2 - sumX * sumX;
    if(_det == 0.) {
        return -1;
    }
    _a = (sumXY * sum - sumX * sumY) / _det;
    _b = (sumX2 * sumY - sumX * sumXY) / _det;    

    _fitted = true;
    return 0;
}

void
TMLine::removeChits() {

    unsigned n = _links.length();
    int nlyr[50] = {0};
    int nneg = 0, npos = 0;
    for (unsigned i = 0; i < n -1 ; i++) {
        TMLink & l = * _links[i];
        nlyr[l.hit()->wire()->layerId()] += 1;
        if(l.position().y() < 0.){
           nneg += 1;
        }else {
           npos += 1;
        }
    }            

    //...
    AList<TMLink> bad;
    for (unsigned i = 0; i < n; i++) {

        TMLink & l = * _links[i];

        //...if # of hits in a wire layer, don't use...
        if (nlyr[l.hit()->wire()->layerId()] > 3) {
            bad.append(l);
            continue;
        }
        //...remove extremely bad poinits
        if(abs(nneg - npos) > 3) {
          if(npos > nneg && l.position().y() < 0 ) bad.append(l);
          if(npos < nneg && l.position().y() > 0 ) bad.append(l);
        }
    }
    //...
    if (bad.length() > 0 && bad.length() < n) {
        _links.remove(bad);
    }

    //... For the next fit
   _fitted = false;
   _fittedUpdated = false;
}

void
TMLine::removeSLY(AList<TMLink> & list) {
      _links.remove(list);
      _fitted = false;
      _fittedUpdated = false;
}

void
TMLine::appendSLY(AList<TMLink> & list) {
      _links.append(list);
      _fitted = false;
      _fittedUpdated = false;
}

void
TMLine::appendByszdistance(AList<TMLink> & list, unsigned isl, float maxSigma) {

    //... intialize
    unsigned nb = _links.length();
    
    //....Select good hit
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
        TMLink & l = * list[i];
        if(l.hit()->wire()->superLayerId() == isl){
          double dist = distance(l);
          if (dist < maxSigma) {
            _links.append(l);
          }
        }
    }

    unsigned na = _links.length();
    if(nb != na){
      AList<TMLink> bad;
      //... remove duplicated hits
      for(unsigned i = 0; i<na ;i++){
         TMLink & l = * _links[i];
         if(i < na - 1) {
           TMLink & lnext = * _links[i+1];
           if(l.hit()->wire()->layerId() == lnext.hit()->wire()->layerId() ){
             if(l.hit()->wire()->localId() == lnext.hit()->wire()->localId()){
               bad.append(l);
             }
           }
         }
      }
      if(bad.length() > 0) _links.remove(bad);
      _fitted = false;
      _fittedUpdated = false;
    }
}

double
TMLine::reducedChi2(void) const {
#ifdef TRKRECO_DEBUG
    if (! _fitted)
        std::cout << "TMLine::reducedChi2 !!! fit not performed" << std::endl;
#endif

    if (_fittedUpdated) return _reducedChi2;
    double chi2 = 0.;
    double scale = 20.;
    unsigned n = _links.length();
    for (unsigned i = 0; i < n; i++) {
        TMLink & l = * _links[i];

        double x = l.position().x();
        double y = l.position().y();
        double c = y - _a * x - _b;
        double err = 1.;
        if (l.hit()) err = scale * l.hit()->dDrift();
        chi2 += c * c / err / err;
    }

    _reducedChi2 = chi2/(n-2);
    _fittedUpdated = true;
    return _reducedChi2;
}

#if defined(__GNUG__)

int
SortByB(const TMLine ** a, const TMLine ** b) {
    if (fabs((* a)->b()) > fabs((* b)->b()))
        return 1;
    else if (fabs((* a)->b()) == fabs((* b)->b()))
        return 0;
    else
        return -1;
}

#else

extern "C" int
SortByB(const void * av, const void * bv) {
    const TMLine ** a((const TMLine **) av);
    const TMLine ** b((const TMLine **) bv);
    if (fabs((* a)->b()) > fabs((* b)->b()))
        return 1;
    else if (fabs((* a)->b()) == fabs((* b)->b()))
        return 0;
    else
        return -1;
}

#endif

---