/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Reconstruction/TrkReco/TrkReco-00-08-59-patch4-slc6tag/src/TSegment0.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TSegment0.cxx,v 1.7 2010/03/31 09:58:59 liucy Exp $
//-----------------------------------------------------------------------------
// Filename : TSegment0.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class to manage a group of TMLink's.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------

#include "TrkReco/TTrack.h"
#include "TrkReco/TSegment0.h"
#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
//#include "cdc/Range.h"
#include "TrkReco/Range.h"

TSegment0::TSegment0()
: TTrackBase(),
  _innerWidth(0),
  _outerWidth(0),
  _nLayer(0),
  _clusterType(0),
  _duality(0.),
  _nDual(0),
  _angle(0.) {
    _fitted = false;
}

TSegment0::TSegment0(const AList<TMLink> & a)
: TTrackBase(a),
  _innerWidth(0),
  _outerWidth(0),
  _nLayer(0),
  _clusterType(0),
  _duality(0.),
  _nDual(0),
  _angle(0.) {
    _links.sort(SortByWireId);
    _fitted = false;
}

TSegment0::~TSegment0() {
}

void
TSegment0::dump(const std::string & msg, const std::string & pre) const {
    if (! _fitted) update();
    bool def = false;
    if (msg == "") def = true;

    if (def || msg.find("cluster") != std::string::npos || msg.find("detail") != std::string::npos) {
        std::cout << pre;
        std::cout << "#links=" << _links.length();
        std::cout << "(" << _innerWidth << "," << _outerWidth << ":";
        std::cout << clusterType() << ")," << _nDual << "," << _duality << ",";
        std::cout << _angle << std::endl;
    }
    if (def || msg.find("vector") != std::string::npos || msg.find("detail") != std::string::npos) {
        std::cout << pre;
        std::cout << "pos" << _position << "," << "dir" << _direction;
        std::cout << std::endl;
    }
    if (! def) TTrackBase::dump(msg, pre);
}

void
TSegment0::update(void) const {
    _clusterType = 0;
    _position = ORIGIN;
    _direction = ORIGIN;
    _inners.removeAll();
    _outers.removeAll();

    if (_links.length() == 0) return;

    _innerMostLayer = _links[0]->wire()->layerId();
    _outerMostLayer = _innerMostLayer;
    for (unsigned i = 1; i < _links.length(); i++) {
        unsigned id = _links[i]->wire()->layerId();
        if (id < _innerMostLayer) _innerMostLayer = id;
        else if (id > _outerMostLayer) _outerMostLayer = id;
    }
    _nLayer = _outerMostLayer - _innerMostLayer + 1;

    double centerX = _links[0]->position().x();
    HepPoint3D inner = ORIGIN;
    HepPoint3D outer = ORIGIN;
    unsigned nInner = 0;
    unsigned nOuter = 0;
    for (unsigned i = 0; i < _links.length(); i++) {
        HepPoint3D tmp = _links[i]->position();
        double diff = tmp.x() - centerX;
        if (diff > M_PI) {
            tmp.setX(centerX - 2. * M_PI + diff);
        }
        else if (diff < - M_PI) {
            tmp.setX(centerX + 2. * M_PI + diff);
        }

        _links[i]->conf(tmp);
        _position += tmp;
        unsigned id = _links[i]->wire()->layerId();
        if (id == _innerMostLayer) {
            inner += tmp;
            ++nInner;
            _inners.append(_links[i]);
        }
        if (id == _outerMostLayer) {
            outer += tmp;
            ++nOuter;
            _outers.append(_links[i]);
        }
    }
    _innerWidth = Width(_inners);
    _outerWidth = Width(_outers);
    _position *= (1. / (float) _links.length());

    inner *= (1. / (float) nInner);
    outer *= (1. / (float) nOuter);
    _direction = (inner - outer).unit();

    _fitted = true;
}

double
TSegment0::distance(const TSegment0 & c) const {
    HepVector3D dir = c.position() - _position;
    if (dir.x() > M_PI) dir.setX(dir.x() - 2. * M_PI);
    else if (dir.x() < - M_PI) dir.setX(2. * M_PI + dir.x());

    float radial = fabs(_direction.dot(dir));
    float radial2 = radial * radial;

    return (dir.mag2() - radial2) > 0.0 ? sqrt(dir.mag2() - radial2) : 0.;
}

double
TSegment0::distance(const HepPoint3D & p, const HepVector3D & v) const {
    HepVector3D dir = _position - p;
    if (dir.x() > M_PI) dir.setX(dir.x() - 2. * M_PI);
    else if (dir.x() < - M_PI) dir.setX(2. * M_PI + dir.x());

    float radial = fabs(v.unit().dot(dir));
    float radial2 = radial * radial;

    return (dir.mag2() - radial2) > 0.0 ? sqrt(dir.mag2() - radial2) : 0.;
}

Range
TSegment0::rangeX(double min, double max) const {
#ifdef TRKRECO_DEBUG_DETAIL
    if (min > max) {
        std::cout << "TSegment0::range !!! bad arguments:min,max=";
        std::cout << min << "," << max << std::endl;
    }
#endif

    unsigned n = _links.length();
    if (n == 0) return Range(0., 0.);

    //...Search for a center...
    bool found = false;
    double center;
    for (unsigned i = 0; i < n; i++) {
        double x = _links[i]->position().x();
        if (x < min || x > max) continue;
        center = x;
        found = true;
        break;
    }
    if (! found) return Range(0., 0.);

#ifdef TRKRECO_DEBUG_DETAIL
//     std::cout << "    center=" << center << std::endl;
#endif

    double distanceR = 0.;
    double distanceL = 0.;
    double distanceMax = max - min;
    for (unsigned i = 0; i < n; i++) {
        double x = _links[i]->position().x();
        if (x < min || x > max) continue;

        double distance0, distance1;
        if (x > center) {
            distance0 = x - center;
            distance1 = distanceMax - distance0;
        }
        else {
            distance1 = center - x;
            distance0 = distanceMax - distance1;
        }

        if (distance0 < distance1) {
            if (distance0 > distanceR) distanceR = distance0;
        }
        else {
            if (distance1 > distanceL) distanceL = distance1;
        }

#ifdef TRKRECO_DEBUG_DETAIL
//      std::cout << "    ";
//      std::cout << _links[i]->wire()->layerId() << "-";
//      std::cout << _links[i]->wire()->localId() << ",";
//      std::cout << _links[i]->position().x();
//      std::cout << ",0,1=" << distance0 << "," << distance1;
//      std::cout << ",l,r=" << distanceL << "," << distanceR;
//      std::cout << std::endl;
#endif
    }

    double right = center + distanceR;
    double left = center - distanceL;

    return Range(left, right);
}

void 
TSegment0::updateType(void) const {
    if (! nLinks()) return;
    if (! _fitted) update();

    //...Parameter...
    unsigned fat = 3;
    unsigned tall = 3;

    //...Calculate
    updateDuality();

    //...Long or short...
    if ((_innerWidth < fat) && (_outerWidth < fat)) {
        _clusterType = 1;

        //...Check length across a super layer...
        if (_nLayer > tall) _clusterType = 2;
        return;
    }

    //...A...
    else if (_outerWidth < fat) {
        _clusterType = 3;
        return;
    }

    //...V...
    else if (_innerWidth < fat) {
        _clusterType = 4;
        return;
    }

    //...X or parallel...
    else {
        bool space = true;
        for (unsigned i = _innerMostLayer; i <= _outerMostLayer; i++) {
            unsigned n = 0;
            AList<TMLink> tmp;
            for (unsigned j = 0; j < _links.length(); j++)
                if (_links[j]->wire()->layerId() == i) {
                    ++n;
                    tmp.append(_links[j]);
                }

            if (n == Width(tmp)) {
                space = false;
                break;
            }
        }

        _clusterType = 5;
        if (space) _clusterType = 6;
        return;
    }
}

AList<TSegment0>
TSegment0::split(void) const {
    AList<TSegment0> list;

    //...Do not split if cluster type is 1, 2, or 7...
    unsigned t = clusterType();
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    ... splitting : type=" << t << std::endl;
#endif 
    if (t == 0) return list;
    else if (t == 2) {
        // beta 5 if (_nDual > 2 && _duality > 0.7 && _angle > 0.7)
        //          return splitDual();
        if (_nDual > 2 && _duality > 0.7) return splitDual();
        return list;
    }
    else if (t == 1) return list;
    else if (t == 7) return list;

    //...Parallel...
    else if (t == 6) return splitParallel();

    //...Avoid splitting of X or parallel...(future implementation)...
    else if (t > 4) return splitComplicated();

    //...A or V...
    return splitAV();
}

AList<TSegment0>
TSegment0::splitAV(void) const {
    unsigned t = clusterType();
    AList<TMLink> seeds[2];

    //...Calculate corner of V or A...
    const AList<TMLink> * corners;
    if (t == 3) corners = & _outers;
    else        corners = & _inners;
    HepPoint3D corner;
    for (unsigned i = 0; i < corners->length(); i++)
        corner += (* corners)[i]->wire()->xyPosition();
    corner *= 1. / (float) corners->length();
    seeds[0].append(* corners);
    seeds[1].append(* corners);

    //...Calculdate edges...
    const AList<TMLink> * links;
    if (t == 3) links = & _inners;
    else        links = & _outers;
    AList<TMLink> edge = Edges(* links);
    HepPoint3D edgePos[2];
    HepVector3D v[2];
    for (unsigned i = 0; i < 2; i++) {
        edgePos[i] = edge[i]->wire()->xyPosition();
        v[i] = (edgePos[i] - corner).unit();
    }
    seeds[0].append(edge[0]);
    seeds[1].append(edge[1]);

#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    corner:" << corner << std::endl;
    std::cout << "    edge:" << edgePos[0] << "(";
    std::cout << edge[0]->wire()->layerId() << "-";
    std::cout << edge[0]->wire()->localId() << ")";
    std::cout << v[0] << std::endl;
    std::cout << "         ";
    std::cout << edgePos[1] << "(";
    std::cout << edge[1]->wire()->layerId() << "-";
    std::cout << edge[1]->wire()->localId() << ")";
    std::cout << v[1] << std::endl;
#endif

    //...Examine each hits...
    unsigned n = _links.length();
    for (unsigned i = 0; i < n; i++) {
        TMLink * l = _links[i];
        if (edge.member(l)) continue;
        if (corners->member(l)) continue;

        HepVector3D p = l->wire()->xyPosition() - corner;
        double p2 = p.mag2();

        double dist[2];
        for (unsigned j = 0; j < 2; j++) {
            dist[j] = v[j].dot(p);
            double d2 = dist[j] * dist[j];
            dist[j] = (p2 - d2) > 0. ? sqrt(p2 - d2) : 0.;
        }
        if (dist[0] < dist[1]) seeds[0].append(l);
        else                   seeds[1].append(l);

#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "        p:" << p << std::endl;
        std::cout << "    " << l->wire()->layerId() << "-";
        std::cout << l->wire()->localId() << ":" << dist[0];
        std::cout << "," << dist[1] << std::endl;
#endif
    }

    AList<TSegment0> list;
    for (unsigned i = 0; i < 2; i++) {
        if (seeds[i].length()) {
            TSegment0 * nc = new TSegment0(seeds[i]);
            AList<TSegment0> ncx = nc->split();
            if (ncx.length() == 0) {
                nc->solveDualHits();
                list.append(nc);
            }
            else {
                list.append(ncx);
                delete nc;
            }
        }
    }
    return list;
}

AList<TSegment0>
TSegment0::splitComplicated(void) const {

    //...Select best hits...
    AList<TSegment0> newClusters;
    AList<TMLink> goodHits;
    unsigned n = _links.length();
    for (unsigned i = 0; i < n; i++) {
        const TMDCWireHit * h = _links[i]->hit();
        unsigned state = h->state();
        if (! (state & WireHitContinuous)) continue;
        if (! (state & WireHitIsolated)) continue;
        if ((! (state & WireHitPatternLeft)) &&
            (! (state & WireHitPatternRight))) continue;
        goodHits.append(_links[i]);
    }
    if (goodHits.length() == 0) return newClusters;

    //...Main loop...
    goodHits.sort(SortByWireId);
    AList<TMLink> original(_links);
    while (goodHits.length()) {

        //...Select an edge hit...
        TMLink * seed = goodHits.last();
        const TMDCWire * wire = seed->wire();
        unsigned localId = wire->localId();
        AList<TMLink> used;
        unsigned nn = _links.length();
        for (unsigned i = 0; i < nn; i++) {
            TMLink * t = _links[i];
            const TMDCWire * w = t->wire();

            //...Straight?...
            if (abs((int) w->localId() - (int) localId) < 2) used.append(t);
        }

#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "        seed is " << seed->wire()->name() << std::endl;
        std::cout << "            ";
        for (unsigned i = 0; i < used.length(); i++) {
            std::cout << used[i]->wire()->name() << ",";
        }
        std::cout << std::endl;
#endif

        //...Create new cluster...
        if (used.length() == 0) continue;
        if (used.length() == nLinks()) return newClusters;
        TSegment0 * c = new TSegment0(used);
        AList<TSegment0> cx = c->split();
        if (cx.length() == 0) {
            c->solveDualHits();
            newClusters.append(c);
        }
        else {
            newClusters.append(cx);
            delete c;
        }
        goodHits.remove(used);
        original.remove(used);
    }

    //...Remainings...
    if ((original.length()) && (original.length() < nLinks())) {
        TSegment0 * c = new TSegment0(original);
        AList<TSegment0> cx = c->split();
        if (cx.length() == 0) {
            c->solveDualHits();
            newClusters.append(c);
        }
        else {
            newClusters.append(cx);
            delete c;
        }
    }

    return newClusters;
}

AList<TSegment0>
TSegment0::splitParallel(void) const {
    AList<TMLink> seeds[2];
    AList<TSegment0> newClusters;
    for (unsigned i = _innerMostLayer; i <= _outerMostLayer; i++) {
        AList<TMLink> list = SameLayer(_links, i);
        AList<TMLink> outerList = Edges(list);

#ifdef TRKRECO_DEBUG_DETAIL
        if (outerList.length() != 2) {
            std::cout << "TSegment0::splitParallel !!! ";
            std::cout << "This is not a parallel cluster" << std::endl;
        }
#endif

        seeds[0].append(outerList[0]);
        seeds[1].append(outerList[1]);
        if (list.length() == 2) continue;

        const TMDCWire & wire0 = * outerList[0]->wire();
        const TMDCWire & wire1 = * outerList[1]->wire();
        for (unsigned k = 0; k < list.length(); k++) {
            TMLink * t = list[k];
            if (t == outerList[0]) continue;
            if (t == outerList[1]) continue;

            if (abs(wire0.localIdDifference(* t->wire())) <
                abs(wire1.localIdDifference(* t->wire())))
                seeds[0].append(t);
            else
                seeds[1].append(t);
        }
    }
    
    if ((seeds[0].length()) && (seeds[0].length() < nLinks())) {
        TSegment0 * c0 = new TSegment0(seeds[0]);
        AList<TSegment0> c0x = c0->split();
        if (c0x.length()) {
            newClusters.append(c0x);
            delete c0;
        }
        else {
            c0->solveDualHits();
            newClusters.append(c0);
        }
    }

    if ((seeds[1].length()) && (seeds[1].length() < nLinks())) {
        TSegment0 * c1 = new TSegment0(seeds[1]);
        AList<TSegment0> c1x = c1->split();
        if (c1x.length()) {
            newClusters.append(c1x);
            delete c1;
        }
        else {
            c1->solveDualHits();
            newClusters.append(c1);
        }
    }

    return newClusters;
}

void
TSegment0::updateDuality(void) const {
    _duality = 0.;
    _nDual = 0;
    HepPoint3D x[2];
    for (unsigned i = _innerMostLayer; i <= _outerMostLayer; i++) {
        AList<TMLink> list = SameLayer(_links, i);
        if (i == _innerMostLayer) {
            for (unsigned j = 0; j < list.length(); j++)
                x[0] += list[j]->hit()->xyPosition();
            x[0] *= 1. / double(list.length());
        }
        else if (i == _outerMostLayer) {
            for (unsigned j = 0; j < list.length(); j++)
                x[1] += list[j]->hit()->xyPosition();
            x[1] *= 1. / double(list.length());
        }

        if (list.length() == 2) {
            if (Width(list) != 2) continue;
            const TMDCWireHit * h0 = list[0]->hit();
            const TMDCWireHit * h1 = list[1]->hit();

            double distance = (h0->xyPosition() - h1->xyPosition()).mag();
            distance = fabs(distance - h0->drift() - h1->drift());
            _duality += distance;
            ++_nDual;
        }
    }
    if (_nDual > 0) _duality /= float(_nDual);
    _angle = (x[1] - x[0]).unit().dot(x[0].unit());
}

AList<TSegment0>
TSegment0::splitDual(void) const {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    TSegment0::splitDual called" << std::endl;
#endif
    AList<TMLink> seeds[2];
    AList<TMLink> unknown;
    for (unsigned i = _innerMostLayer; i <= _outerMostLayer; i++) {
        AList<TMLink> list = SameLayer(_links, i);

        if (list.length() == 2) {
            if (Width(list) == 2) {
                seeds[0].append(list[0]);
                seeds[1].append(list[1]);
                continue;
            }
        }
        unknown.append(list);
    }

    if (unknown.length() > 0) {
#ifdef TRKRECO_DEBUG_DETAIL
        if (seeds[0].length() == 0)
            std::cout << "TSegment0::splitDual !!! no TMLink for seed 0" << std::endl;
        if (seeds[1].length() == 0)
            std::cout << "TSegment0::splitDual !!! no TMLink for seed 1" << std::endl;
#endif  
        const HepPoint3D & p0 = seeds[0][0]->xyPosition();
        const HepPoint3D & p1 = seeds[1][0]->xyPosition();
        HepVector3D v0 = (seeds[0].last()->xyPosition() - p0).unit();
        HepVector3D v1 = (seeds[1].last()->xyPosition() - p1).unit();

        for (unsigned i = 0; i < unknown.length(); i++) {
            TMLink * t = unknown[i];
            HepVector3D x0 = t->xyPosition() - p0;
            double d0 = (x0 - (x0.dot(v0) * v0)).mag();
            HepVector3D x1 = t->xyPosition() - p1;
            double d1 = (x1 - (x1.dot(v0) * v0)).mag();

            if (d0 < d1) seeds[0].append(t);
            else         seeds[1].append(t);
        }
    }

    AList<TSegment0> newClusters;
    newClusters.append(new TSegment0(seeds[0]));
    newClusters.append(new TSegment0(seeds[1]));
    return newClusters;
}

int
TSegment0::solveDualHits(void) {
    updateType();
    if (_clusterType == 0) return 0;
    if (_nDual == 0) return 0;

    AList<TMLink> seeds;
    AList<TMLink> duals;
    for (unsigned i = _innerMostLayer; i <= _outerMostLayer; i++) {
        AList<TMLink> list = SameLayer(_links, i);

        if (list.length() == 1) {
            seeds.append(list[0]);
        }
        else if (list.length() == 2) {
            if (Width(list) > 1) {
                const TMDCWireHit * h0 = list[0]->hit();
                const TMDCWireHit * h1 = list[1]->hit();
                double distance = (h0->xyPosition() - h1->xyPosition()).mag();
                distance = fabs(distance - h0->drift() - h1->drift());
                if (distance > 0.5) duals.append(list);
#ifdef TRKRECO_DEBUG_DETAIL
//              h0->dump();
//              h1->dump();
//              std::cout << "duality distance = " << distance << std::endl;
//              std::cout << "i = " << i << std::endl;
#endif
            }
        }
        else if (list.length() == 0) {
            continue;
        }
#ifdef TRKRECO_DEBUG_DETAIL
        else {
            std::cout << "TSegment0::solveDualHits !!! this is not expected 2";
            std::cout << std::endl;
            this->dump("cluster hits mc", "    ");
        }
#endif
    }

    //...Solve them...
    if (seeds.length() < 2) return -1;
    AList<TMLink> outers = InOut(seeds);
    const HepPoint3D & p = outers[0]->xyPosition();
    HepVector3D v = (outers[1]->xyPosition() - p).unit();
    unsigned n = duals.length() / 2;
    for (unsigned i = 0; i < n; i++) {
        TMLink * t0 = duals[i * 2 + 0];
        TMLink * t1 = duals[i * 2 + 1];
        HepVector3D x0 = t0->xyPosition() - p;
        HepVector3D x1 = t1->xyPosition() - p;
        double d0 = (x0 - (x0.dot(v) * v)).mag();
        double d1 = (x1 - (x1.dot(v) * v)).mag();
        if (d0 < d1) _links.remove(t1);
        else         _links.remove(t0);
    }
    return n;
}

/*
#include "tuple/BelleTupleManager.h"

void
CheckSegments(const CAList<TSegment0> & list) {
    static bool first = true;
    static BelleHistogram * hError;
    static BelleHistogram * hNHeps[11];
    static BelleHistogram * hNHits[3];

    if (first) {
        first = false;
        extern BelleTupleManager * BASF_Histogram;
        BelleTupleManager * m = BASF_Histogram;

        hError = m->histogram("segment errors", 10, 0, 10, 20000);

        hNHeps[0] = m->histogram("segment nheps sl 0", 10, 0., 10.);
        hNHeps[1] = m->histogram("segment nheps sl 1", 10, 0., 10.);
        hNHeps[2] = m->histogram("segment nheps sl 2", 10, 0., 10.);
        hNHeps[3] = m->histogram("segment nheps sl 3", 10, 0., 10.);
        hNHeps[4] = m->histogram("segment nheps sl 4", 10, 0., 10.);
        hNHeps[5] = m->histogram("segment nheps sl 5", 10, 0., 10.);
        hNHeps[6] = m->histogram("segment nheps sl 6", 10, 0., 10.);
        hNHeps[7] = m->histogram("segment nheps sl 7", 10, 0., 10.);
        hNHeps[8] = m->histogram("segment nheps sl 8", 10, 0., 10.);
        hNHeps[9] = m->histogram("segment nheps sl 9", 10, 0., 10.);
        hNHeps[10] = m->histogram("segment nheps sl 10", 10, 0., 10.);

        hNHits[0] = m->histogram("segment nhits, nheps = 1", 20, 0., 20.);
        hNHits[1] = m->histogram("segment nhits, nheps = 2", 20, 0., 20.);
        hNHits[2] = m->histogram("segment nhits, nheps >= 3", 20, 0., 20.);

        std::cout << "CheckSegments ... initialized" << std::endl;
        return;
    }

    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
        const TSegment0 & s = * list[i];
        const AList<TMLink> & links = s.links();
        unsigned nLinks = links.length();
        if (nLinks == 0) {
            hError->accumulate(0.5);
            continue;
        }

        unsigned sl = links[0]->wire()->superLayerId();
        unsigned nHeps = s.nHeps();
        hNHeps[sl]->accumulate((float) nHeps + .5);
        if (nHeps == 1)      hNHits[0]->accumulate((float) nLinks + .5);
        else if (nHeps == 2) hNHits[1]->accumulate((float) nLinks + .5);
        else                 hNHits[2]->accumulate((float) nLinks + .5);
    }
}

void
CheckSegmentLink(const TSegment0 & base,
                 const TSegment0 & next,
                 float distance,
                 float dirAngle) {
    static bool first = true;
    static BelleHistogram * hAngle[2];
    static BelleHistogram * hAngle1[2];
    static BelleHistogram * hDistance[2];
    static BelleHistogram * hDirAngle[2];
    static BelleHistogram * hDirAngle1[2];

    if (first) {
        first = false;
        extern BelleTupleManager * BASF_Histogram;
        BelleTupleManager * m = BASF_Histogram;

        hAngle[0] = m->histogram("segment correct link, angle",
                                 50, -1., 1., 21000);
        hAngle[1] = m->histogram("segment wrong link, angle", 50, -1., 1.);
        hAngle1[0] = m->histogram("segment correct link, angle, wide",
                                  50, .8, 1.);
        hAngle1[1] = m->histogram("segment wrong link, angle, wide",
                                  50, .8, 1.);
        hDistance[0] = m->histogram("segment correct link, dist", 50, 0., 1.);
        hDistance[1] = m->histogram("segment wrong link, dist", 50, 0., 1.);
        hDirAngle[0] = m->histogram("segment correct link, dir angle",
                                    50, -1, 1.);
        hDirAngle[1] = m->histogram("segment wrong link, dir angle",
                                    50, -1, 1.);
        hDirAngle1[0] = m->histogram("segment correct link, dir angle, wide",
                                    50, .8, 1.);
        hDirAngle1[1] = m->histogram("segment wrong link, dir angle, wide",
                                    50, .8, 1.);

        std::cout << "CheckSegmentLink ... initialized" << std::endl;
        return;
    }

    const TTrackHEP * const hep0 = base.hep();
    const TTrackHEP * const hep1 = next.hep();

    float angle = base.direction().dot(next.direction());

    unsigned id = 0;
    if (hep0 != hep1) id = 1;
    hAngle[id]->accumulate(angle);
    hAngle1[id]->accumulate(angle);
    hDistance[id]->accumulate(distance);
    hDirAngle[id]->accumulate(dirAngle);
    hDirAngle1[id]->accumulate(dirAngle);
}
*/

unsigned
NCoreLinks(const CAList<TSegment0> & list) {
    unsigned n = 0;
    unsigned nList = list.length();
    for (unsigned i = 0; i < nList; i++) {
        const AList<TMLink> & links = list[i]->links();
        for (unsigned j = 0; j < links.length(); j++) {
            unsigned state = links[j]->hit()->state();
            if ((! (state & WireHitPatternLeft)) &&
                (! (state & WireHitPatternRight)))
                ++n;
        }
    }
    return n;
}

AList<TMLink>
Links(const TSegment0 & s, const TTrack & t) {
    AList<TMLink> a;

    const AList<TMLink> & links = s.links();
    const AList<TMLink> & trackLinks = t.links();
    unsigned n = links.length();
    for (unsigned i = 0; i < n; i++) {
        if (trackLinks.hasMember(links[i]))
            a.append(links[i]);
    }

    return a;
}

---