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