/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Reconstruction/TrkReco/TrkReco-00-08-59-patch4-slc6tag/src/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    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// 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;
}

---