TTrackBase.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TTrackBase.cxx,v 1.19 2011/09/05 17:32:47 maoh Exp $
//-----------------------------------------------------------------------------
// Filename : TTrackBase.h
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : [email protected]
//-----------------------------------------------------------------------------
// Description : A virtual class for a track class in tracking.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#define HEP_SHORT_NAMES
#include "CLHEP/Alist/CList.h"
#include "TrkReco/TTrackBase.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TTrackHEP.h"
#include "TrkReco/TMFitter.h"
 
TTrackBase::TTrackBase(const AList<TMLink> & a)
: _links(a),
  _mc(0),
//    _superLayer(0),
//    _nSuperLayers(0),
  _fitter(0),
  _updated(false),
  _fitted(false),
  _time(0.) {
//      _layer[0] = 0;
//      _layer[1] = 0;
//      _nLayers[0] = 0;
//      _nLayers[1] = 0;
//    update();
}
 
TTrackBase::TTrackBase()
: _mc(0),
//    _superLayer(0),
//    _nSuperLayers(0),
  _fitter(0),
  _updated(true),
  _fitted(false),
  _time(0.) {
//      _layer[0] = 0;
//      _layer[1] = 0;
//      _nLayers[0] = 0;
//      _nLayers[1] = 0;
}
 
TTrackBase::~TTrackBase() {
}
 
void
TTrackBase::dump(const std::string & msg, const std::string & pre) const {
    bool mc = (msg.find("mc") != std::string::npos);
    bool pull = (msg.find("pull") != std::string::npos);
    bool flag = (msg.find("flag") != std::string::npos);
    bool detail = (msg.find("detail") != std::string::npos);
    if (detail)
    mc = pull = flag = true;
 
    if (detail || (msg.find("layer") != std::string::npos)) {
    if (! _updated) update();
    }
    if (detail || (msg.find("hits") != std::string::npos)) {
    Dump(_links, msg, pre);
    }
}
 
void
TTrackBase::update(void) const {
    _cores.removeAll();
    unsigned n = _links.length();
    for (unsigned i = 0; i < n; i++) {
    TMLink * l = _links[i];
    const TMDCWireHit & h = * l->hit();
    if (h.state() & WireHitInvalidForFit) continue;
    if (! (h.state() & WireHitFittingValid)) continue;
    _cores.append(l);
    }
    _updated = true;
}
 
double
TTrackBase::distance(const TMLink &) const {
    std::cout << "TTrackBase::distance !!! not implemented" << std::endl;
    return 0.;
}
 
int
TTrackBase::approach(TMLink &) const {
    std::cout << "TTrackBase::approach !!! not implemented" << std::endl;
    return -1;
}
 
void
TTrackBase::appendByApproach(AList<TMLink> & list, double maxSigma) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    TTrackBase::appendByApproach ... sigma=" << maxSigma << std::endl;
#endif
 
    AList<TMLink> unused;
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
    TMLink & l = * list[i];
 
    if ((_links.hasMember(l)) || (l.hit()->state() & WireHitUsed))
        continue;
 
    //...Calculate closest approach...
    int err = approach(l);
    if (err < 0) {
        unused.append(l);
        continue;
    }
 
    //...Calculate sigma...
    float distance = (l.positionOnWire() - l.positionOnTrack()).mag();
    float diff = fabs(distance - l.drift());
    float sigma = diff / l.dDrift();
 
    //...For debug...
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    sigma=" << sigma;
    std::cout << ",dist=" << distance;
    std::cout << ",diff=" << diff;
    std::cout << ",err=" << l.hit()->dDrift() << ",";
    if (sigma < maxSigma) std::cout << "ok,";
    else                  std::cout << "X,";
    l.dump("mc");
#endif
        
    //...Make sigma cut...
    if (sigma > maxSigma) {
        unused.append(l);
        continue;
    }
 
    //...OK...
    _links.append(l);
    _updated = false;
    _fitted = false;
    }
    list.remove(unused);
}
 
void
TTrackBase::appendByDistance(AList<TMLink> & list, double sigma) {
    std::cout << "TTrackBase::appendByDistance !!! not implemented" << std::endl;
    list.removeAll();
}
 
AList<TMLink> 
TTrackBase::refineMain(double sigma) {
    AList<TMLink> bad;
    unsigned n = _links.length();
    for (unsigned i = 0; i < n; i++){
//  std::cout<<" pull of link = "<<_links[i]->pull()<<std::endl;
    if (_links[i]->pull() > sigma)
        bad.append(_links[i]);
    }
    return bad;
}
 
void
TTrackBase::refine(AList<TMLink> & list, double sigma) {
    AList<TMLink> bad = refineMain(sigma);
#ifdef TRKRECO_DEBUG
    std::cout << "    refine ... sigma=" << sigma << ", # of rejected hits=";
    std::cout << bad.length() << std::endl;
#endif
#ifdef TRKRECO_DEBUG
    Dump(bad, "sort pull mc", "        ");
#endif
 
    if (bad.length()) {
    _links.remove(bad);
    list.append(bad);
    _fitted = false;
    _updated = false;
    }
}
 
void
TTrackBase::removeLinks(void){
  _links.removeAll();
}
 
void
TTrackBase::refine(double sigma) {
    AList<TMLink> bad = refineMain(sigma);
//      for (unsigned i = 0; i < bad.length(); i++) {
//      const TMDCWireHit * hit = bad[i]->hit();
//      hit->state(hit->state() | WireHitInvalidForFit);
//      }
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    refine ... sigma=" << sigma << std::endl;
    Dump(bad, "detail sort", "        ");
#endif
 
    if (bad.length()) {
    _fitted = false;
    _updated = false;
    }
}
int 
TTrackBase::DropWorst() {
  AList<TMLink> bad;
  int jbad=-1;
  double sigma=0;
  unsigned n = _links.length();
  double chi2=0;
  for (unsigned i = 0; i < n; i++){
   // if(g_sigmalyr1)g_sigmalyr1->fill(sqrt(_links[i]->pull()),_links[i]->wire()->layerId());
    chi2+=_links[i]->pull();
    if (_links[i]->pull() > sigma){
      jbad=i; 
      sigma = _links[i]->pull();
    }
  }
        
  chi2=chi2/n;
  //if(sigma>9){bad.append(_links[jbad]);if(g_sigmalyr)g_sigmalyr->fill(sqrt(sigma),_links[jbad]->wire()->layerId());}
//  if(sigma>20.25&&chi2>2.08){bad.append(_links[jbad]);}
  if(sigma>9){bad.append(_links[jbad]);}
  if (bad.length()) {
    _links.remove(bad);
    _fitted = false;
    _updated = false;
  }
  return bad.length();
}
 
unsigned
TTrackBase::testByApproach(const AList<TMLink> & list, double maxSigma) const {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    TTrackBase::testByApproach ... sigma=" << maxSigma << std::endl;
#endif
 
    unsigned nOK = 0;
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
    TMLink & l = * list[i];
    nOK += testByApproach(l, maxSigma);
    }
    return nOK;
}
 
unsigned
TTrackBase::testByApproach(const TMLink & link, double maxSigma) const {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    TTrackBase::testByApproach ... sigma=" << maxSigma << std::endl;
#endif
    TMLink l = link;
 
    //...Calculate closest approach...
    int err = approach(l);
    if (err < 0) return 0;
    //...Calculate sigma...
    float distance = l.distance();
    float diff = fabs(distance - l.hit()->drift());
    float sigma = diff / l.hit()->dDrift();
    l.pull(sigma * sigma);
 
    //...For debug...
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    sigma=" << sigma;
    std::cout << ",dist=" << distance;
    std::cout << ",diff=" << diff << ",";
    if (sigma < maxSigma) std::cout << "ok,";
    else                  std::cout << "X,";
    l.dump("mc");
#endif
        
    //...Make sigma cut...
    if (sigma < maxSigma) return 1;
 
    return 0;
}
 
//  unsigned
//  TTrackBase::nAxialHits(void) const {
//      unsigned n = _links.length();
//      unsigned a = 0;
//      for (unsigned i = 0; i < n; i++) {
//      if (_links[i]->wire()->axial()) ++a;
//      }
//      return a;
//  }
 
const AList<TMLink> &
TTrackBase::links(unsigned mask) const {
    if (mask == 0) return _links;
 
    std::cout << "TTrackBase::links !!! mask is not supportted yet" << std::endl;
    return _links;
}
 
unsigned
TTrackBase::nLinks(unsigned mask) const {
    unsigned n = _links.length();
    if (mask == 0) return n;
    unsigned nn = 0;
    for (unsigned i = 0; i < n; i++) {
    const TMDCWireHit & h = * _links[i]->hit();
    if (h.state() & mask) ++nn;
    }
    return nn;
}
 
const AList<TMLink> &
TTrackBase::cores(unsigned mask) const {
    if (mask)
    std::cout << "TTrackBase::cores !!! mask is not supported" << std::endl;
    if (! _updated) update();
    return _cores;
}
 
unsigned
TTrackBase::nCores(unsigned mask) const {
    if (mask)
    std::cout << "TTrackBase::nCores !!! mask is not supported" << std::endl;
    if (! _updated) update();
    return _cores.length();
}
 
//  TMLink *
//  TTrackBase::innerMostLink(unsigned mask) const {
//      if (mask == 0) return InnerMost(_links);
//      unsigned n = _links.length();
//      AList<TMLink> tmp;
//      for (unsigned i = 0; i < n; i++) {
//      const TMDCWireHit & h = * _links[i]->hit();
//      if (h.state() & mask) tmp.append(_links[i]);
//      }
//      return InnerMost(tmp);
//  }
 
//  TMLink *
//  TTrackBase::outerMostLink(unsigned mask) const {
//      if (mask == 0) return OuterMost(_links);
//      unsigned n = _links.length();
//      AList<TMLink> tmp;
//      for (unsigned i = 0; i < n; i++) {
//      const TMDCWireHit & h = * _links[i]->hit();
//      if (h.state() & mask) tmp.append(_links[i]);
//      }
//      return OuterMost(tmp);
//  }
 
int
TTrackBase::fit(void) {
    return _fitter->fit(* this);
}
 
void
TTrackBase::append(TMLink & a) {
    if ((a.hit()->state() & WireHitUsed)) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "TTrackBase::append !!! " << a.wire()->name()
          << " this is already used by another track!" << std::endl;
#endif
    return;
    }
    if (_links.hasMember(a)) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "TTrackBase::append !!! " << a.wire()->name()
          << " this is already included in this track!" << std::endl;
#endif
    return;
    }
    _links.append(a);
    _updated = false;
    _fitted = false;
    _fittedWithCathode = false; // added by matsu ( 1999/05/24 )
}
 
void
TTrackBase::append(const AList<TMLink> & a) {
    AList<TMLink> tmp;
    for (unsigned i = 0; i < a.length(); i++) {
    if ((_links.hasMember(a[i])) || (a[i]->hit()->state() & WireHitUsed)) {
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "    TTrackBase::append(list) !!! ";
        std::cout << a[i]->wire()->name();
        std::cout << " Hey!, this is already used! Don't mess with me.";
        std::cout << std::endl;
#endif
        continue;
    }
    else {
        tmp.append(a[i]);
    }
    }
    _links.append(tmp);
    _updated = false;
    _fitted = false;
    _fittedWithCathode = false; // added by matsu ( 1999/05/24 )
}
 
const TTrackHEP * const
TTrackBase::hep(void) const {
    unsigned n = _links.length();
    CAList<TTrackHEP> hepList;
    CList<unsigned> hepCounter;
    for (unsigned i = 0; i < n; i++) {
    const TTrackHEP * hep = _links[i]->hit()->mc()->hep();
    unsigned nH = hepList.length();
    bool found = false;
    for (unsigned j = 0; j < nH; j++) {
        if (hepList[j] == hep) {
        found = true;
        ++(* hepCounter[j]);
        }
    }
 
    if (! found) {
        hepList.append(hep);
        unsigned c = 0;
        hepCounter.append(c);
    }
    }
 
    _nHeps = hepList.length();
    _hep = 0;
    unsigned max = 0;
    for (unsigned i = 0; i < _nHeps; i++) {
    if ((* hepCounter[i]) > max) {
        max = (* hepCounter[i]);
        _hep = hepList[i];
    }
    }
 
    return _hep;
}
 
unsigned
TTrackBase::nHeps(void) const {
    hep();
    return _nHeps;
}