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    : [email protected]
//-----------------------------------------------------------------------------
// 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;
}