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

---