TMLink.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TMLink.cxx,v 1.15 2011/10/08 01:56:15 maoh Exp $
//-----------------------------------------------------------------------------
// Filename : TMLink.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : [email protected]
//-----------------------------------------------------------------------------
// Description : A class to relate TMDCWireHit and TTrack objects.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#include "TrkReco/TMLink.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TTrackHEP.h"
#include "CLHEP/Alist/ConstAList.h"
 
//zangsl 040518 : move inlined function here
const TMDCWire * const
TMLink::wire(void) const {
    if (_hit) return _hit->wire();
    return NULL;
}
 
const HepPoint3D &
TMLink::xyPosition(void) const {
    return _hit->wire()->xyPosition();
}
//zangsl 040518 end
 
/*TMLink::TMLink(TTrack * t, const TMDCWireHit * h, const HepPoint3D & p)
: _track(t),
  _hit(h),
  _dPhi(0),
  _leftRight(0),
  _pull(0),
  _position(p),
  _link(0),
  _fit2D(0) {
//  _usecathode(0) {
    if (h) {
    _drift[0] = h->drift(0);
    _drift[1] = h->drift(1);
    _dDrift[0] = h->dDrift(0);
    _dDrift[1] = h->dDrift(1);
    }
    else {
    _drift[0] = 0.;
    _drift[1] = 0.;
    _dDrift[0] = 0.;
    _dDrift[1] = 0.;
    }
 
    for (unsigned i = 0; i < 6; ++i)
    _neighbor[i] = NULL;
    for (unsigned i = 0; i < 4; ++i)
      _arcZ[i];
 
    if (h) {
    _onTrack = _onWire = h->xyPosition();
    }
}
*/
//for Tsf
TMLink::TMLink(TTrack * t, const TMDCWireHit * h, const HepPoint3D & p, const HepPoint3D & d, double dr)
: _track(t),
  _hit(h),
  _dPhi(0),
  _leftRight(0),
  _pull(0),
  _position(p),
  _positionD(d),
  _link(0),
  _fit2D(0),
  _tsfTag(0) {
    if (h) {
        _cDrift[0] = dr;   //after conformal transformation.
        _cDrift[1] = dr;
        _drift[0] = h->drift(0);
        _drift[1] = h->drift(1);
        _dDrift[0] = h->dDrift(0);
        _dDrift[1] = h->dDrift(1);
    }
    else {
        _cDrift[0] = 0.;
        _cDrift[1] = 0.;
        _drift[0] = 0.;
        _drift[1] = 0.;
        _dDrift[0] = 0.;
        _dDrift[1] = 0.;
    }
                                                                                                                                                                                                                                                           
    for (unsigned i = 0; i < 6; ++i)
        _neighbor[i] = NULL;
    for (unsigned i = 0; i < 4; ++i)
      _arcZ[i];
                                                                                                                             
                                                                                                                             
    if (h) {
        _onTrack = _onWire = h->xyPosition();
    }
}
 
TMLink::TMLink(const TMLink & l)
: _track(l._track),
  _hit(l._hit),
  _onTrack(l._onTrack),
  _onWire(l._onWire),
  _position(l._position),
  _positionD(l._positionD),
  _dPhi(l._dPhi),
  _leftRight(l._leftRight),
  _pull(l._pull),
  _link(l._link),
  _distance(l._distance),   
  // addition by matsu ( 1999/07/05 )
//  _mclust(l._mclust ),
//  _usecathode(l._usecathode ),
  // end of addition
  _fit2D(l._fit2D),
  _tsfTag(l._tsfTag) {
    _drift[0] = l._drift[0];
    _drift[1] = l._drift[1];
    _dDrift[0] = l._dDrift[0];
    _dDrift[1] = l._dDrift[1];
    _cDrift[0] = l._cDrift[0];
    _cDrift[1] = l._cDrift[1];
    for (unsigned i = 0; i < 6; ++i)
    _neighbor[i] = l._neighbor[i];
    for (unsigned i = 0; i < 4; ++i)
      _arcZ[i] = l._arcZ[i];
 
}
 
TMLink::~TMLink() {
}
 
void
TMLink::dump(const std::string & msg, const std::string & pre) const {
    std::cout << pre;
    if (_track) std::cout << "track#=,";
    if (_hit) {
    _hit->dump(msg);
    }
}
 
unsigned
NLayers(const AList<TMLink> & list) {
    unsigned l0 = 0;
    unsigned l1 = 0;
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
    unsigned id = list[i]->wire()->layerId();
    if (id < 32) l0 |= (1 << id);
    else         l1 |= (1 << (id - 32));
    }
 
    unsigned l = 0;
    for (unsigned i = 0; i < 32; i++) {
    if (l0 & (1 << i)) ++l;
    if (l1 & (1 << i)) ++l;
    }
    return l;
}
 
void
NHits(const AList<TMLink> & links, unsigned nHits[43]) {
    for (unsigned i = 0; i < 43; i++) nHits[i] = 0;
    unsigned nLinks = links.length();
    for (unsigned i = 0; i < nLinks; i++)
    ++nHits[links[i]->wire()->layerId()];
}
 
void
NHitsSuperLayer(const AList<TMLink> & links, unsigned nHits[11]) {
    for (unsigned i = 0; i < 11; i++) nHits[i] = 0;
    unsigned nLinks = links.length();
    for (unsigned i = 0; i < nLinks; i++)
    ++nHits[links[i]->wire()->superLayerId()];
}
 
void
Dump(const CAList<TMLink> & links, const std::string & msg, const std::string & pre) {
    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 sort = (msg.find("sort") != std::string::npos);
    bool stereo = (msg.find("stereo") != std::string::npos);
    bool detail = (msg.find("detail") != std::string::npos);
    bool pos = (msg.find("position") != std::string::npos);
    if (detail)
    mc = pull = flag = sort = true;
 
    CAList<TMLink> tmp = links;
    if (sort)
    tmp.sort(SortByWireId);
    unsigned n = tmp.length();
    unsigned nForFit = 0;
#define MCC_MAX 1000
    unsigned MCC0[MCC_MAX];
    unsigned MCC1[MCC_MAX];
    for (unsigned i = 0; i < MCC_MAX; i++) {
    MCC0[i] = 0;
    MCC1[i] = 0;
    }
    bool MCCOverFlow = false;
 
    std::cout << pre;
    for (unsigned i = 0; i < n; i++) {
    const TMLink & l = * tmp[i];
    std::cout << l.wire()->name();
 
    double a = l.pull();
    unsigned mcId = 0;
    if (mc)
        if (l.hit()->mc())
        if (l.hit()->mc()->hep())
            mcId = l.hit()->mc()->hep()->id();
    if (pull) {
        std::cout << "[" << a << "]";
    }
    if (mc) {
        std::cout << "(" << mcId << ")";
        if (mcId < MCC_MAX) {
        ++MCC0[mcId];
        if (l.hit()->state() & WireHitFittingValid) {
            if (! (l.hit()->state() & WireHitInvalidForFit))
            ++MCC1[mcId];
        }
        }
        else {
        MCCOverFlow = true;
        }
    }
    if (flag) {
        if (l.hit()->state() & WireHitFindingValid)
        std::cout << "o";
        if (l.hit()->state() & WireHitFittingValid) {
        std::cout << "+";
        if (! (l.hit()->state() & WireHitInvalidForFit))
            ++nForFit;
        }
        if (l.hit()->state() & WireHitInvalidForFit)
        std::cout << "x";
    }
    if (stereo) {
        std::cout << "{" << l.leftRight() << "," << l.zStatus() << "}";
    }
    if (pos) {
        std::cout << ",pos=" << l.position();
    }
    std::cout << ",";
    }
    std::cout << " " << n << " l(s)";
    if (flag) std::cout << ", fv " << nForFit << " l(s)";
    if (mc) {
    unsigned nMC = 0;
    std::cout << ", mc";
    for (unsigned i = 0; i < MCC_MAX; i++) {
        if (MCC0[i] > 0) {
        ++nMC;
        std::cout << i << ":" << MCC0[i] << ",";
        }
    }
    std::cout << " total " << nMC << " contributions";
    if (flag) {
        nMC = 0;
        std::cout << ", fv mc";
        for (unsigned i = 0; i < MCC_MAX; i++) {
        if (MCC1[i] > 0) {
            ++nMC;
            std::cout << i << ":" << MCC1[i] << ",";
        }
        }
        std::cout << " total " << nMC << " contributions";
    }
 
    if (MCCOverFlow)
        std::cout << "(counter overflow)";
    }
    std::cout << std::endl;
}
 
void
Dump(const TMLink & link, const std::string & msg, const std::string & pre) {
    CAList<TMLink> tmp;
    tmp.append(link);
    Dump(tmp, msg, pre);
}
 
unsigned
NStereoHits(const AList<TMLink> & links) {
    unsigned nLinks = links.length();
    unsigned n = 0;
    for (unsigned i = 0; i < nLinks; i++)
    if (links[i]->wire()->stereo())
        ++n;
    return n;
}
 
unsigned
NAxialHits(const AList<TMLink> & links) {
    unsigned nLinks = links.length();
    unsigned n = 0;
    for (unsigned i = 0; i < nLinks; i++)
    if (links[i]->wire()->axial())
        ++n;
    return n;
}
 
AList<TMLink>
AxialHits(const AList<TMLink> & links) {
    AList<TMLink> a;
    unsigned n = links.length();
    for (unsigned i = 0; i < n; i++) {
    if (links[i]->wire()->axial())
        a.append(links[i]);
    }
    return a;
}
 
AList<TMLink>
StereoHits(const AList<TMLink> & links) {
    AList<TMLink> a;
    unsigned n = links.length();
    for (unsigned i = 0; i < n; i++) {
    if (! links[i]->wire()->axial())
        a.append(links[i]);
    }
    return a;
}
 
TMLink *
InnerMost(const AList<TMLink> & a) {
    unsigned n = a.length();
    unsigned minId = 9999;
    TMLink * t = 0;
    for (unsigned i = 0; i < n; i++) {
    unsigned id = a[i]->wire()->id();
    if (id < minId) {
        minId = id;
        t = a[i];
    }
    }
    return t;
}
 
TMLink *
OuterMost(const AList<TMLink> & a) {
    unsigned n = a.length();
    unsigned maxId = 0;
    TMLink * t = 0;
    for (unsigned i = 0; i < n; i++) {
    unsigned id = a[i]->wire()->id();
    if (id > maxId) {
        maxId = id;
        t = a[i];
    }
    }
    return t;
}
 
void
SeparateCores(const AList<TMLink> & input,
          AList<TMLink> & cores,
          AList<TMLink> & nonCores) {
    unsigned n = input.length();
    for (unsigned i = 0; i < n; i++) {
    TMLink & t = * input[i];
    const TMDCWireHit & h = * t.hit();
    if (h.state() & WireHitFittingValid)
        cores.append(t);
    else
        nonCores.append(t);
    }
}
 
AList<TMLink>
Cores(const AList<TMLink> & input) {
    AList<TMLink> a;
    unsigned n = input.length();
    for (unsigned i = 0; i < n; i++) {
    TMLink & t = * input[i];
    const TMDCWireHit & h = * t.hit();
    if (h.state() & WireHitFittingValid)
        a.append(t);
    }
    return a;
}
 
#if defined(__GNUG__)
int
SortByWireId(const TMLink ** a, const TMLink ** b) {
    if ((* a)->wire()->id() > (* b)->wire()->id()) return 1;
    else if
    ((* a)->wire()->id() == (* b)->wire()->id()) return 0;
    else return -1;
}
 
int
SortByX(const TMLink ** a, const TMLink ** b) {
    if ((* a)->position().x() > (* b)->position().x()) return 1;
    else if ((* a)->position().x() == (* b)->position().x()) return 0;
    else return -1;
}
 
#else
extern "C" int
SortByWireId(const void * av, const void * bv) {
  const TMLink ** a((const TMLink**)av);
  const TMLink ** b((const TMLink**)bv);
    if ((* a)->wire()->id() > (* b)->wire()->id()) return 1;
    else if
    ((* a)->wire()->id() == (* b)->wire()->id()) return 0;
    else return -1;
}
 
extern "C" int
SortByX(const void* av, const void* bv) {
  const TMLink ** a((const TMLink**)av);
  const TMLink ** b((const TMLink**)bv);
    if ((* a)->position().x() > (* b)->position().x()) return 1;
    else if ((* a)->position().x() == (* b)->position().x()) return 0;
    else return -1;
}
 
#endif
 
unsigned
Width(const AList<TMLink> & list) {
    unsigned n = list.length();
    if (n < 2) return n;
 
    const TMDCWire * w = list[0]->wire();
    unsigned nWires = w->layer()->nWires();
    unsigned center = w->localId();
 
#ifdef TRKRECO_DEBUG_DETAIL
    unsigned sId = w->superLayerId();
#endif
 
    unsigned left = 0;
    unsigned right = 0;
    for (unsigned i = 1; i < n; i++) {
    w = list[i]->wire();
    unsigned id = w->localId();
 
    unsigned distance0, distance1;
    if (id > center) {
        distance0 = id - center;
        distance1 = nWires - distance0;
    }
    else {
        distance1 = center - id;
        distance0 = nWires - distance1;
    }
 
    if (distance0 < distance1) {
        if (distance0 > right) right = distance0;
    }
    else {
        if (distance1 > left) left = distance1;
    }
 
#ifdef TRKRECO_DEBUG_DETAIL
    if (w->superLayerId() != sId)
        std::cout << "::width !!! super layer assumption violation" << std::endl;
#endif
    }
 
    return right + left + 1;
}
 
AList<TMLink>
Edges(const AList<TMLink> & list) {
    AList<TMLink> a;
 
    unsigned n = list.length();
    if (n < 2) return a;
    else if (n == 2) return list;
 
    const TMDCWire * w = list[0]->wire();
    unsigned nWires = w->layer()->nWires();
    unsigned center = w->localId();
 
    unsigned left = 0;
    unsigned right = 0;
    TMLink * leftL = list[0];
    TMLink * rightL = list[0];
    for (unsigned i = 1; i < n; i++) {
    w = list[i]->wire();
    unsigned id = w->localId();
 
    unsigned distance0, distance1;
    if (id > center) {
        distance0 = id - center;
        distance1 = nWires - distance0;
    }
    else {
        distance1 = center - id;
        distance0 = nWires - distance1;
    }
 
    if (distance0 < distance1) {
        if (distance0 > right) {
        right = distance0;
        rightL = list[i];
        }
    }
    else {
        if (distance1 > left) {
        left = distance1;
        leftL = list[i];
        }
    }
    }
 
    a.append(leftL);
    a.append(rightL);
    return a;
}
 
AList<TMLink>
SameLayer(const AList<TMLink> & list, const TMLink & a) {
    AList<TMLink> same;
    unsigned id = a.wire()->layerId();
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
    if (list[i]->wire()->layerId() == id) same.append(list[i]);
    }
    return same;
}
 
AList<TMLink>
SameSuperLayer(const AList<TMLink> & list, const TMLink & a) {
    AList<TMLink> same;
    unsigned id = a.wire()->superLayerId();
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
    if (list[i]->wire()->superLayerId() == id) same.append(list[i]);
    }
    return same;
}
 
AList<TMLink>
SameLayer(const AList<TMLink> & list, unsigned id) {
    AList<TMLink> same;
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
    if (list[i]->wire()->layerId() == id) same.append(list[i]);
    }
    return same;
}
 
AList<TMLink>
SameSuperLayer(const AList<TMLink> & list, unsigned id) {
    AList<TMLink> same;
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
    if (list[i]->wire()->superLayerId() == id) same.append(list[i]);
    }
    return same;
}
 
AList<TMLink>
InOut(const AList<TMLink> & list) {
    AList<TMLink> inners;
    AList<TMLink> outers;
    unsigned n = list.length();
    unsigned innerMostLayer = 999;
    unsigned outerMostLayer = 0;
    for (unsigned i = 0; i < n; i++) {
    unsigned id = list[i]->wire()->layerId();
    if (id < innerMostLayer) innerMostLayer = id;
    else if (id > outerMostLayer) outerMostLayer = id;
    }
    for (unsigned i = 0; i < n; i++) {
    unsigned id = list[i]->wire()->layerId();
    if (id == innerMostLayer) inners.append(list[i]);
    else if (id == outerMostLayer) outers.append(list[i]);
    }
    inners.append(outers);
    return inners;
}
 
unsigned
SuperLayer(const AList<TMLink> & list) {
    unsigned sl = 0;
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++)
    sl |= (1 << (list[i]->wire()->superLayerId()));
    return sl;
}
 
unsigned
SuperLayer(const AList<TMLink> & links, unsigned minN) {
    unsigned n = links.length();
    unsigned nHits[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
    for (unsigned i = 0; i < n; i++)
    ++nHits[links[i]->wire()->superLayerId()];
    unsigned sl = 0;
    for (unsigned i = 0; i < 11; i++)
    if (nHits[i] >= minN)
        sl |= (1 << i);
    return sl;
}
 
unsigned
NSuperLayers(const AList<TMLink> & list) {
    unsigned l0 = 0;
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
    unsigned id = list[i]->wire()->superLayerId();
    l0 |= (1 << id);
    }
 
    unsigned l = 0;
    for (unsigned i = 0; i < 11; i++) {
    if (l0 & (1 << i)) ++l;
    }
    return l;
}
 
unsigned
NSuperLayers(const AList<TMLink> & links, unsigned minN) {
    unsigned n = links.length();
    unsigned nHits[11] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
    for (unsigned i = 0; i < n; i++)
    ++nHits[links[i]->wire()->superLayerId()];
    unsigned sl = 0;
    for (unsigned i = 0; i < 11; i++)
    if (nHits[i] >= minN)
        ++sl;
    return sl;
}
 
unsigned
NMissingAxialSuperLayers(const AList<TMLink> & links) {
    unsigned n = links.length();
//Liuqg, change the following to BES.    unsigned nHits[6] = {0, 0, 0, 0, 0, 0};
    unsigned nHits[5] = {0, 0, 0, 0, 0};
    for (unsigned i = 0; i < n; i++)
    if (links[i]->wire()->axial())
        ++nHits[links[i]->wire()->axialStereoLayerId() / 4];
    unsigned j = 0;
    while (nHits[j] == 0) ++j;
    unsigned nMissing = 0;
    unsigned nMax = 0;
    for (unsigned i = j; i < 5; i++) {
    if (nHits[i] == 0) ++nMissing;
    else {
        if (nMax < nMissing) nMax = nMissing;
        nMissing = 0;
    }
    }
    return nMax;
}
 
const TTrackHEP &
Links2HEP(const AList<TMLink> & links) {
    const TTrackHEP * best = NULL;
    const AList<TTrackHEP> & list = TTrackHEP::list();
    unsigned nHep = list.length();
 
    if (! nHep) return * best;
 
    unsigned * N = (unsigned *) malloc(nHep * sizeof(unsigned));
    for (unsigned i = 0; i < nHep; i++) N[i] = 0;
 
    for (unsigned i = 0; i < links.length(); i++) {
    const TMLink & l = * links[i];
    const TTrackHEP & hep = * l.hit()->mc()->hep();
    for (unsigned j = 0; j < nHep; j++)
        if (list[j] == & hep)
        ++N[j];
    }
 
    unsigned nMax = 0;
    for (unsigned i = 0; i < nHep; i++) {
    if (N[i] > nMax) {
        best = list[i];
        nMax = N[i];
    }
    }
 
    return * best;
}
 
 
/*
double
TMLink::DriftTime(double _tof,double z) const {
    double tprop = 0.;
    double _vprop = (_layer<8) ? Constants::vpropInner : Constants::vpropOuter;
    if (0 == _layer%2){
      tprop = (0.5*_zlen + z)/_vprop; //odd
    }else{
      tprop = (0.5*_zlen - z)/_vprop; //even
    }
    double driftT;
    driftT = fabs(_rawTime - _T0Walk -1.e9*tof - tprop);
 
    //std::cout<< "lay "<<_layer<<" cell "<<_wire<<" zhit "<<z<<" tprop "<<tprop << std::endl;
    return driftT;
}
*/