TBuilder.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TBuilder.cxx,v 1.26 2012/05/28 05:16:29 maoh Exp $
//-----------------------------------------------------------------------------
// Filename : TBuilder.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : [email protected]
//-----------------------------------------------------------------------------
// Description : A class to build a track.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#include "CLHEP/String/Strings.h"
#include "TrkReco/TBuilder.h"
#include "TrkReco/TMDC.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TCircle.h"
#include "TrkReco/TLine0.h"
#include "TrkReco/TMLine.h"
#include "TrkReco/TTrack.h"
#include "TrkReco/TSegment.h"
#include "TrkReco/TPoint2D.h"
#include "TrkReco/TLine2D.h"
#include "TrkReco/TRobustLineFitter.h"
#ifdef TRKRECO_DEBUG
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TTrackHEP.h"
#endif
#ifdef TRKRECO_WINDOW
#include "TrkReco/TWindow.h"
TWindow sz("sz");
#endif
 
TBuilder::TBuilder(const std::string & a,
           float maxSigma,
           float maxSigmaStereo,
           float salvageLevel,
           float szSegmentDistance,
           float szLinkDistance,
           unsigned fittingFlag)
: _name(a),
  _fitter("TBuilder Fitter"),
  _maxSigma(maxSigma),
  _maxSigmaStereo(maxSigmaStereo),
  _salvageLevel(sqrt(salvageLevel)),
//  _minNCores(3), //20060307
  _minNCores(2), //liuqg20070312
  _szSegmentDistance(szSegmentDistance),
  _szLinkDistance(szLinkDistance) {
    if (fittingFlag & 1) _fitter.sag(true);
    if (fittingFlag & 2) _fitter.propagation(true);
    if (fittingFlag & 4) _fitter.tof(true);
    if (fittingFlag & 8) _fitter.freeT0(true);
}
 
TBuilder::~TBuilder() {
}
 
void
TBuilder::dump(const std::string & msg, const std::string & pre) const {
}
 
TTrack *
TBuilder::buildRphi(const AList<TSegment> & list) const {
 
#ifdef TRKRECO_DEBUG
    std::cout << "... building rphi by segments : # of segments = ";
    std::cout << list.length() << std::endl;
    for (unsigned i = 0; i < list.length(); i++)
    list[i]->dump("hits sort flag", "    ");
#endif
    //yuany
    //    for (unsigned i = 0; i < list.length(); i++)
    //    list[i]->dump("", "    ");
 
    //...Pick up links...
    AList<TMLink> links = Links(list);
    //yuany
//    for(unsigned i=0;i<links.length();i++){
//      cout<<"wire name "<<links[i]->wire()->name()<<endl;
//    }
 
    //...Main funtion...
    TTrack * t = buildRphi(links);
 
    //...Check used segments...
    if (t) {
    const AList<TMLink> & usedLinks = t->links();
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
        TSegment & segment = * list[i];
        AList<TMLink> used = Links(segment, * t);
        if (used.length()) {
        t->segments().append(segment);
        segment.tracks().append(t);
        }
    }
    }
 
    return t;
 
}
 
TTrack *
TBuilder::buildRphi(const AList<TMLink> & list) const {
#ifdef TRKRECO_DEBUG
    std::cout << "... building rphi by links : # of links = ";
    std::cout << list.length() << std::endl;
#endif
 
    //...Classify TMLink's...
    AList<TMLink> cores;
    AList<TMLink> nonCores;
    SeparateCores(list, cores, nonCores);
 
#ifdef TRKRECO_DEBUG
    cout<<"  ... cores..."<<endl;
    for(int ii = 0; ii < cores.length(); ++ii) {
    cout<<"layer: "<<cores[ii]->wire()->layerId()
        <<"  local: "<<cores[ii]->wire()->localId()<<endl;
    }
    cout<<"  ...noncores..."<<endl;
    for(int ii = 0; ii < nonCores.length(); ++ii) {
        cout<<"layer: "<<nonCores[ii]->wire()->layerId()
            <<"  local: "<<nonCores[ii]->wire()->localId()<<endl;
    }
#endif
    //...Check # of links...
    unsigned nCores = cores.length();
    if (nCores < _minNCores) {
#ifdef TRKRECO_DEBUG
    std::cout << "... building rphi failure : # of cores(=" << nCores;
    std::cout << ") is less then " << _minNCores << std::endl;
#endif
 
    return NULL;
    }
 
    //...Make a circle...
#ifdef TRKRECO_DEBUG
    std::cout <<"links in list = "<<list.length()<<std::endl;
    std::cout << "... making a circle : # cores =" << cores.length() << std::endl;
#endif
    TCircle c(cores);
    int err = c.fit();
    if (err < 0) {
#ifdef TRKRECO_DEBUG
    std::cout << "... building rphi failure : circle fit error = ";
    std::cout << err << std::endl;
#endif
    return NULL;
    }
//  cout<<"radius "<<c.radius()<<" charge "<<c.charge()<<" center  "<<c.center()<<endl;
    //...Make a track...
    TTrack * t = new TTrack(c);
//   t->dump("hits sort flag", "    ");
    AList<TMLink> bad;
    err = _fitter.fit(* t);
    //cout<<"fit 1 done err: "<<err<<endl;
    // t->dump("hits sort flag", "    ");
    if (err < 0) goto discard;
    t->refine(bad, _maxSigma * 100.);
    //cout<<"refine 1 done"<<endl;
    //t->dump("hits sort flag", "    ");
    err = _fitter.fit(* t); 
    //cout<<"fit 2 done err: "<<err<<endl;
    //t->dump("hits sort flag", "    ");
    t->refine(bad, _maxSigma * 10.);
    //cout<<"refine 2 done"<<endl;
    //t->dump("hits sort flag", "    ");
    err = _fitter.fit(* t); 
    //cout<<"fit 3 done err: "<<err<<endl;
    //t->dump("hits sort flag", "    ");
    t->refine(bad, _maxSigma);
    //cout<<"refine 3 done"<<endl;
    //t->dump("hits sort flag", "    ");
    err = _fitter.fit(* t);
    //cout<<"fit 4 done err: "<<err<<endl;
    
    if (err < 0) goto discard;
#ifdef TRKRECO_DEBUG_DETAIL
    c.dump("detail", "    ccl> ");
    t->dump("detail", "    1st> ");
#endif
 
    //...Try to append non-core hits...
#ifdef TRKRECO_DEBUG
    std::cout << "... appending non-core hits : # = " << nonCores.length();
    std::cout << std::endl;
#endif
    t->appendByApproach(nonCores, _salvageLevel);
#ifdef TRKRECO_DEBUG
    t->dump("hits sort flag", "    ");
#endif
//   t->dump("hits sort flag", "    ");
    return t;
    
    //...Something happened...
discard:
#ifdef TRKRECO_DEBUG
    std::cout << "... building rphi failure : helix fit error = ";
    std::cout << err << std::endl;
#endif
    delete t;
    return NULL;
}
 
void
TBuilder::salvage(TTrack & t, AList<TMLink> & hits) const {
 
#ifdef TRKRECO_DEBUG
    std::cout << "... salvaging(TBuilder) : # of given hits=" << hits.length();
    std::cout << ", salvage level=" << _salvageLevel << std::endl;
    Dump(hits, "hits sort flag");
#endif
 
    unsigned nHits = hits.length();
    if (nHits == 0) return;
 
    //...Try to append this hit...
    t.appendByApproach(hits, _salvageLevel);
    _fitter.fit(t);
}
 
void
TBuilder::salvage(TTrack & t, AList<TSegment> & list) const {
 
#ifdef TRKRECO_DEBUG
    std::cout << "... salvaging by segments : # of segments = ";
    std::cout << list.length() << std::endl;
    for (unsigned i = 0; i < list.length(); i++)
    list[i]->dump("hits sort flag", "    ");
#endif
 
    //...Pick up links...
    AList<TMLink> links;
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++)
    links.append(list[i]->links());
 
    salvage(t, links);
}
 
TMLine *
TBuilder::initialLine1(const TTrack & t,
               const AList<TSegment> & list,
               const AList<TMLink> & lList) const {
 
    //...Check input...
    if (list.length() < 2) return NULL;
 
    //...Check super layer pattern...
    AList<TSegment> sl[5];
    AList<TMLink> tl[5];
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i ++) {
    unsigned j = list[i]->superLayerId() / 2;
    sl[j].append(list[i]);
    tl[j].append(lList[i]);
    }
#ifdef TRKRECO_DEBUG
    std::cout << "        ... initialLine1 : super layer ptn = ";
    for (unsigned i = 0; i < 5; i++) std::cout << sl[i].length();
    std::cout << std::endl;
#endif
 
    //...Count single segment layer...
    unsigned nSingle = 0;
    for (unsigned i = 0; i < 5; i++)
    if (sl[i].length() == 1) ++nSingle;
#ifdef TRKRECO_DEBUG
    std::cout << "        ... # of single segment layer = " << nSingle << std::endl;
#endif
    if (nSingle < 2) return NULL;
 
    //...Happy case...
    AList<TSegment> bestCombination;
    AList<TMLink> forLine;
    for (unsigned i = 0; i < 5; i++) {
    if (sl[i].length() != 1) continue;
    bestCombination.append(sl[i]);
    forLine.append(tl[i]);
    }
    TMLine & line = * new TMLine(forLine);
    int err = line.fit();
    if (err) {
    delete & line;
    return NULL;    
    }
    return & line;
}
 
TMLine *
TBuilder::initialLine2(const TTrack & t,
               const AList<TMLink> & lList) const {
#ifdef TRKRECO_DEBUG
    std::cout << "        ... initlialLine2 : # of links = " << lList.length();
    std::cout << std::endl;
#endif
 
    //...Static object...
    static TRobustLineFitter fitter("can you work?");
    TMLine * line = new TMLine(lList);
    int err = fitter.fit(* line);
    if (err) {
    delete line;
    return NULL;
    }
    return line;
}
 
TMLine *
TBuilder::initialLine(const TTrack & t, AList<TSegment> & list) const {
 
    //...Select good segments in acceptance...
    AList<TSegment> bad;
    AList<TMLink> lList;
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
    TMLink & l = * new TMLink();
    int err = t.szPosition(* list[i], l);
    if (err) {
        bad.append(list[i]);
        delete & l;
        continue;
    }
 
    const HepPoint3D & p = l.position();
    if (p.y() > l.wire()->forwardPosition().z() ||
        p.y() < l.wire()->backwardPosition().z()) {
        bad.append(list[i]);
        delete & l;
        continue;
    }
    lList.append(l);
    }
    list.remove(bad);
    n = list.length();
 
    //...Cal. expected # of super layers...
    unsigned nAMaxSL = 0;
    unsigned nAMinSL = 10;
    const AList<TSegment> & segments = t.segments();
    unsigned nA = segments.length();
    for (unsigned i = 0; i < nA; i++) {
    unsigned sl = segments[i]->superLayerId();
    if (sl > nAMaxSL) nAMaxSL = sl;
    if (sl < nAMinSL) nAMinSL = sl;
    }
    unsigned nExpected = (nAMaxSL - nAMinSL) / 2;
#ifdef TRKRECO_DEBUG
    std::cout << "    ... initialLine : axial super layer usage = " << nAMinSL;
    std::cout << " ~ " << nAMaxSL << std::endl;
    std::cout << "                    : expected stereo super layers = ";
    std::cout << nExpected << std::endl;
#endif
 
    TMLine * line = NULL;
    AList<TSegment> list0 = list;
    AList<TMLink> lList0 = lList;
    while (1) {
    bool last = (NSuperLayers(lList) <= nExpected);
 
    TMLine * line = initialLine2(t, lList);
    if (line) {
        AList<TSegment> tmp = selectStereoSegment(* line, list0, lList0);
#ifdef TRKRECO_DEBUG
        std::cout << "    ... initialLine2 : # of selected segments = ";
        std::cout << tmp.length() << std::endl;
#endif
        if ((tmp.length() >= nExpected) ||
        (tmp.length() >= 4)         ||
        last) {
        list = tmp;
        return line;
        }
        delete line;
    }
 
    line = initialLine1(t, list, lList);
    if (line) {
        AList<TSegment> tmp = selectStereoSegment(* line, list0, lList0);
        bool ok = (tmp.length() >= nExpected);
#ifdef TRKRECO_DEBUG
        std::cout << "    ... initialLine2 : # of selected segments = ";
        std::cout << tmp.length() << " : ok, last = ";
        std::cout << ok << ", " << last << std::endl;
#endif
        if (ok) {
        list = tmp;
        return line;
        }
        removeFarSegment(* line, list, lList);
        delete line;
    }
    else {
        return NULL;
    }
    }
}
 
TMLine *
TBuilder::initialLineOld(const TTrack & t, AList<TSegment> & list) const {
 
    //...Select good segments in acceptance...
    AList<TSegment> bad;
    AList<TMLink> lList;
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
    TMLink & l = * new TMLink();
    int err = t.szPosition(* list[i], l);
    if (err) {
        bad.append(list[i]);
        delete & l;
        continue;
    }
 
    const HepPoint3D & p = l.position();
    if (p.y() > l.wire()->forwardPosition().z() ||
        p.y() < l.wire()->backwardPosition().z()) {
        bad.append(list[i]);
        delete & l;
        continue;
    }
    lList.append(l);
    }
    list.remove(bad);
    n = list.length();
 
    //...Cal. expected # of super layers...
    unsigned nAMaxSL = 0;
    unsigned nAMinSL = 10;
    const AList<TSegment> & segments = t.segments();
    unsigned nA = segments.length();
    for (unsigned i = 0; i < nA; i++) {
    unsigned sl = segments[i]->superLayerId();
    if (sl > nAMaxSL) nAMaxSL = sl;
    if (sl < nAMinSL) nAMinSL = sl;
    }
    unsigned nExpected = (nAMaxSL - nAMinSL) / 2;
#ifdef TRKRECO_DEBUG
    std::cout << "    ... initialLine : axial super layer usage = " << nAMinSL;
    std::cout << " ~ " << nAMaxSL << std::endl;
    std::cout << "                    : expected stereo super layers = ";
    std::cout << nExpected << std::endl;
#endif
 
    TMLine * line = NULL;
    AList<TSegment> list0 = list;
    AList<TMLink> lList0 = lList;
    while (1) {
    TMLine * line = initialLine1(t, list, lList);
    if (line) {
        AList<TSegment> tmp = selectStereoSegment(* line, list0, lList0);
#ifdef TRKRECO_DEBUG
        std::cout << "    ... initialLine1 : # of selected segments = ";
        std::cout << tmp.length() << std::endl;
#endif
        if (tmp.length() >= nExpected) {
        lList0.remove(line->links());
        HepAListDeleteAll(lList0);
        list = tmp;
        return line;
        }
        delete line;
    }
 
    line = initialLine2(t, lList);
    if (line) {
        AList<TSegment> tmp = selectStereoSegment(* line, list0, lList0);
        bool ok = (tmp.length() >= nExpected);
        bool last = (NSuperLayers(lList) <= nExpected);
#ifdef TRKRECO_DEBUG
        std::cout << "    ... initialLine2 : # of selected segments = ";
        std::cout << tmp.length() << " : ok, last = ";
        std::cout << ok << ", " << last << std::endl;
#endif
        if (ok || last) {
        lList0.remove(line->links());
        HepAListDeleteAll(lList0);
        list = tmp;
        return line;
        }
        removeFarSegment(* line, list, lList);
        delete line;
    }
    else {
        HepAListDeleteAll(lList0);
        return NULL;
    }
    }
}
 
AList<TSegment>
TBuilder::selectStereoSegment(const TMLine & line,
                  const AList<TSegment> & list,
                  const AList<TMLink> & szList) const {
    AList<TSegment> outList;
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
    float distance = line.distance(* szList[i]);
    if (distance < _szSegmentDistance) outList.append(list[i]);
    }
    return outList;
}
 
void
TBuilder::removeFarSegment(const TMLine & line,
               AList<TSegment> & list,
               AList<TMLink> & szList) const {
    unsigned n = list.length();
    float maxDistance = 0.;
    unsigned maxId = 0;
    for (unsigned i = 0; i < n; i++) {
    float distance = line.distance(* szList[i]);
    if (distance > maxDistance) maxId = i;
    }
    list.remove(maxId);
    szList.remove(maxId);
}
 
TTrack *
TBuilder::buildStereo(TTrack & track,
              TMLine & line,
              const AList<TMLink> & list) const {
#ifdef TRKRECO_DEBUG
    std::cout << "... building stereo by links : # of links = ";
    std::cout << list.length() << std::endl;
#endif
 
    //...Classify TMLink's...
    AList<TMLink> cores;
    AList<TMLink> nonCores;
    SeparateCores(list, cores, nonCores);
    //...Check # of links...
    unsigned nCores = cores.length();
    if (nCores < _minNCores) {
#ifdef TRKRECO_DEBUG
    std::cout << "... stereo building failure : # of cores(=" << nCores;
    std::cout << ") is less then " << _minNCores << std::endl;
#endif
    return NULL;
    }
#ifdef TRKRECO_WINDOW
    sz.appendSz(track, cores, leda_brown);
#endif
 
    //...Cal. left and right position...
    AList<TMLink> forNewLine;
    for (unsigned i = 0; i < nCores; i++) {
    TMLink & t = * cores[i];
    for (unsigned i = 0; i < 2; i++) {
        TMLink & tt = * new TMLink(t);
        tt.leftRight(i);
        int err = track.szPosition(tt);
        if (err) {
        delete & tt;
        }
        else {
        if (line.distance(tt) < _szLinkDistance) {
            tt.link(& t);
            forNewLine.append(tt);
        }
        else {
            delete & tt;
        }
        }
    }
    }
 
    //...Create new line...
#ifdef TRKRECO_DEBUG
    std::cout << "    ... creating a new line" << std::endl;
#endif  
    unsigned nNewLine = forNewLine.length();
    TMLine newLine(forNewLine);
    static TRobustLineFitter fitter("abc");
    int err = fitter.fit(newLine);
    // int err = newLine.fit();
    if (err < 0) {
    HepAListDeleteAll(forNewLine);
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    ... 2nd linear fit failure. nLinks(";
    std::cout << forNewLine.length() << ")" << std::endl;
#endif  
    return NULL;
    }
 
#ifdef TRKRECO_WINDOW
    sz.append(newLine, leda_green);
    sz.wait();
#endif
 
#ifdef TRKRECO_DEBUG_DETAIL
    Dump(forNewLine, "sort hits stereo", "        ");
#endif
 
    //...Decide left/right...
    AList<TMLink> toRemove;
    TMLink * last = NULL;
    for (unsigned i = 0; i < nNewLine; i++) {
    if (last == NULL) last = forNewLine[i]->link();
    else {
        if (last == forNewLine[i]->link()) {
        if (newLine.distance(* forNewLine[i - 1]) >
            newLine.distance(* forNewLine[i]))
            toRemove.append(forNewLine[i - 1]);
        else
            toRemove.append(forNewLine[i]);
        last = NULL;
        }
        else {
        last = forNewLine[i]->link();
        }
    }
    }
    forNewLine.remove(toRemove);
    nNewLine = forNewLine.length();
 
#ifdef TRKRECO_DEBUG_DETAIL
    Dump(toRemove, "sort hits stereo", "      x ");
    Dump(forNewLine, "sort hits stereo", "        ");
#endif
 
    //...3D fit...
    for (unsigned i = 0; i < nNewLine; i++)
    track.append(* forNewLine[i]->link());
    Vector a(5);
    a = track.helix().a();
    a[3] = newLine.b();
    a[4] = track.charge() * newLine.a();
    track._helix->a(a);
 
    //...Refine...
    AList<TMLink> bad;
    err = _fitter.fit(track);
    track.refine(bad, _maxSigmaStereo * 100.);
    err = _fitter.fit(track);
    track.refine(bad, _maxSigmaStereo * 10.);
    err = _fitter.fit(track);
    track.refine(bad, _maxSigmaStereo);
    err = _fitter.fit(track);
 
#ifdef TRKRECO_WINDOW
    sz.text("stereo finished");
    sz.oneShot(track, leda_blue);
#endif
 
    //...Termination...
    HepAListDeleteAll(toRemove);
    HepAListDeleteAll(forNewLine);
    return & track;
}
 
TTrack *
TBuilder::buildStereo(const TTrack & t0, AList<TSegment> & segments) const {
 
    TTrack & t = * new TTrack(t0);
 
#ifdef TRKRECO_DEBUG
    std::cout << "... building stereo by segments : # of segments = ";
    std::cout << segments.length() << std::endl;
    for (unsigned i = 0; i < segments.length(); i++)
    segments[i]->dump("hits sort flag", "    ");
#endif
#ifdef TRKRECO_WINDOW
    sz.clear();
    sz.skip(false);
    sz.mode(2);
    sz.appendSz(t, segments, leda_black);
    AList<TSegment> tmps = segments;
    std::string s = "# of segments = " + itostring(int(segments.length()));
#endif
 
    //...Find initial line...
    TMLine * line = initialLineOld(t, segments);
    if (! line) {
#ifdef TRKRECO_DEBUG
    std::cout << "... building stereo failure : no initial line found" << std::endl;
#endif
#ifdef TRKRECO_WINDOW
    s = "no initial line found : " + s;
    sz.text(s);
    sz.wait();
#endif
    return NULL;
    }
#ifdef TRKRECO_WINDOW
    sz.append(* line, leda_red);
    sz.text(s);
    sz.wait();
#endif
 
    //...Pick up links...
    AList<TMLink> links;
    unsigned n = segments.length();
    for (unsigned i = 0; i < n; i++)
    links.append(segments[i]->links());
 
    //...Main funtion...
    TTrack * ts = buildStereo(t, * line, links);
 
    //...Check used segments...
    if (ts) {
    AList<TMLink> usedLinks = StereoHits(t.links());
    for (unsigned i = 0; i < n; i++) {
        TSegment & segment = * segments[i];
        AList<TMLink> used = Links(segment, t);
        if (used.length()) {
        t.segments().append(segment);
        segment.tracks().append(t);
        }
    }
    }
 
    HepAListDeleteAll((AList<TMLink> &) line->links());
    delete line;
    return ts;
}
 
AList<TMLine>
TBuilder::searchInitialLines(unsigned nSuperLayerMax) const {
#ifdef TRKRECO_DEBUG
    std::cout << "    ... searchInitialLines : # of segments in sl : ";
    for (unsigned i = 0; i < 5; i++)
    std::cout << _links[i].length() << ",";
    std::cout << std::endl
     << "                           : max # of super layers="
     << nSuperLayerMax << std::endl;
#endif
 
    AList<TMLine> lines;
    if (nSuperLayerMax > 5)
        lines.append(searchLines6());
    else if (nSuperLayerMax > 4)
    lines.append(searchLines5());
    else if (nSuperLayerMax > 3)
    lines.append(searchLines4());
    else if (nSuperLayerMax > 2)
    lines.append(searchLines3());
    else if (nSuperLayerMax > 1)
    lines.append(searchLines2());
//      else
//      lines.append(searchLines1());
 
    lines.sort(SortByB);
    return lines;
}
 
AList<TMLine>
TBuilder::searchLines6(void) const {
    unsigned n[6];
    for (unsigned i = 0; i < 6; i++)
        n[i] = _links[i].length();
                                                                                                                             
    AList<TMLine> lines;
    for (unsigned i0 = 0; i0 < n[0]; i0++) {
        for (unsigned i1 = 0; i1 < n[1]; i1++) {
            for (unsigned i2 = 0; i2 < n[2]; i2++) {
                for (unsigned i3 = 0; i3 < n[3]; i3++) {
                    for (unsigned i4 = 0; i4 < n[4]; i4++) {
            for (unsigned i5 = 0; i5 < n[5]; i5++) {
                                                                                                                             
                        AList<TMLink> forLine;
                        forLine.append(_links[0][i0]);
                        forLine.append(_links[1][i1]);
                        forLine.append(_links[2][i2]);
                        forLine.append(_links[3][i3]);
                        forLine.append(_links[4][i4]);
                        forLine.append(_links[5][i5]);
                                                                                                                             
                        TMLine & line = * new TMLine(forLine);
                        int err = line.fit();
                        if (err) {
                            delete & line;
                            continue;
                        }
                                                                                                                             
                        lines.append(line);
            }
                    }
                }
            }
        }
    }
                                                                                                                             
#ifdef TRKRECO_DEBUG
    std::cout << "        ... searchLines5 : # of lines found = "
         << lines.length() << std::endl;
#endif
                                                                                                                             
    return lines;
}
 
AList<TMLine>
TBuilder::searchLines5(void) const {
    AList<TMLine> lines;
 
    const AList<TMLink> * l[5];
    for (unsigned skip = 0; skip < 6; skip++) {
    if (skip == 0) {
        l[0] = & _links[1];
        l[1] = & _links[2];
        l[2] = & _links[3];
        l[3] = & _links[4];
        l[4] = & _links[5];
    }
    else if (skip == 1) {
        l[0] = & _links[0];
    }
    else if (skip == 2) {
        l[1] = & _links[1];
    }
    else if (skip == 3) {
        l[2] = & _links[2];
    }
    else if (skip == 4) {
        l[3] = & _links[3];
    }
    else if (skip == 5) {
        l[4] = & _links[4];
    }
 
    unsigned n[5];
    for (unsigned i = 0; i < 5; i++)
    n[i] = l[i]->length();
 
    for (unsigned i0 = 0; i0 < n[0]; i0++) {
    for (unsigned i1 = 0; i1 < n[1]; i1++) {
        for (unsigned i2 = 0; i2 < n[2]; i2++) {
        for (unsigned i3 = 0; i3 < n[3]; i3++) {
            for (unsigned i4 = 0; i4 < n[4]; i4++) {
 
            AList<TMLink> forLine;
            forLine.append((* l[0])[i0]);
            forLine.append((* l[1])[i1]);
            forLine.append((* l[2])[i2]);
            forLine.append((* l[3])[i3]);
            forLine.append((* l[4])[i4]);
 
            TMLine & line = * new TMLine(forLine);
            int err = line.fit();
            if (err) {
                delete & line;
                continue;
            }
 
            lines.append(line);
            }
        }
        }
    }
    }
    }
 
#ifdef TRKRECO_DEBUG
    std::cout << "        ... searchLines5 : # of lines found = "
     << lines.length() << std::endl;
#endif
 
    return lines;
}
 
AList<TMLine>
TBuilder::searchLines4(void) const {
    AList<TMLine> lines;
 
    const AList<TMLink> * l[4];
    for (unsigned skip = 0; skip < 15; skip++) {
    if (skip == 0) {
        l[0] = & _links[0];
        l[1] = & _links[1];
        l[2] = & _links[2];
        l[3] = & _links[3];
    }
    else if (skip == 1) {
        l[3] = & _links[4];
    }
    else if (skip == 2) {
        l[3] = & _links[5];
    }
    else if (skip == 3) {
        l[2] = & _links[3];
        l[3] = & _links[4];
    }
    else if (skip == 4) {
        l[3] = & _links[5];
    }
    else if (skip == 5) {
        l[2] = & _links[4];
    }
    else if (skip == 6) {
        l[1] = & _links[2];
        l[2] = & _links[3];
        l[3] = & _links[4];
    }
    else if (skip == 7) {
        l[3] = & _links[5];
    }
    else if (skip == 8) {
        l[2] = & _links[4];
    }
    else if (skip == 9) {
        l[1] = & _links[3];
    }
    else if (skip == 10) {
        l[0] = & _links[1];
        l[1] = & _links[2];
        l[2] = & _links[3];
        l[3] = & _links[4];
    }
    else if (skip == 11) {
        l[3] = & _links[5];
    }
    else if (skip == 12) {
        l[2] = & _links[4];
    }
    else if (skip == 13) {
        l[1] = & _links[3];
    }
    else if (skip == 14) {
        l[0] = & _links[2];
    }
 
    unsigned n[4];
    for (unsigned i = 0; i < 4; i++)
        n[i] = l[i]->length();
 
    for (unsigned i0 = 0; i0 < n[0]; i0++) {
        for (unsigned i1 = 0; i1 < n[1]; i1++) {
        for (unsigned i2 = 0; i2 < n[2]; i2++) {
            for (unsigned i3 = 0; i3 < n[3]; i3++) {
            AList<TMLink> forLine;
            forLine.append((* l[0])[i0]);
            forLine.append((* l[1])[i1]);
            forLine.append((* l[2])[i2]);
            forLine.append((* l[3])[i3]);
 
            TMLine & line = * new TMLine(forLine);
            int err = line.fit();
            if (err) {
                delete & line;
                continue;
            }
 
            lines.append(line);
            }
        }
        }
    }
    }
 
#ifdef TRKRECO_DEBUG
    std::cout << "        ... searchLines4 : # of lines found = "
     << lines.length() << std::endl;
#endif
 
    return lines;
}
 
AList<TMLine>
TBuilder::searchLines3(void) const {
    AList<TMLine> lines;
 
    const AList<TMLink> * l[3];
    for (unsigned skip = 0; skip < 20; skip++) {
    if (skip == 0) {
        l[0] = & _links[0];
        l[1] = & _links[1];
        l[2] = & _links[2];
    }
    else if (skip == 1) {
        l[2] = & _links[3];
    }
    else if (skip == 2) {
        l[2] = & _links[4];
    }
    else if (skip == 3) {
        l[2] = & _links[5];
    }
    else if (skip == 4) {
        l[1] = & _links[2];
    }
    else if (skip == 5) {
        l[2] = & _links[4];
    }
    else if (skip == 6) {
        l[2] = & _links[3];
    }
    else if (skip == 7) {
        l[1] = & _links[3];
        l[2] = & _links[4];
    }
    else if (skip == 8) {
        l[2] = & _links[5];
    }
    else if (skip == 9) {
        l[1] = & _links[4];
    }
    else if (skip == 10) {
        l[0] = & _links[1];
    }
    else if (skip == 11) {
        l[1] = & _links[3];
    }
    else if (skip == 12) {
        l[2] = & _links[4];
    }
    else if (skip == 13) {
        l[1] = & _links[2];
    }
    else if (skip == 14) {
        l[2] = & _links[3];
    }
    else if (skip == 15) {
        l[2] = & _links[5];
    }
    else if (skip == 16) {
        l[0] = & _links[2];
        l[1] = & _links[3];
    }
    else if (skip == 17) {
        l[2] = & _links[3];
    }
    else if (skip == 18) {
        l[1] = & _links[4];
        l[2] = & _links[5];
    }
    else if (skip == 19) {
        l[0] = & _links[3];
    }
 
    unsigned n[3];
    for (unsigned i = 0; i < 3; i++)
        n[i] = l[i]->length();
 
    for (unsigned i0 = 0; i0 < n[0]; i0++) {
        for (unsigned i1 = 0; i1 < n[1]; i1++) {
        for (unsigned i2 = 0; i2 < n[2]; i2++) {
            AList<TMLink> forLine;
            forLine.append((* l[0])[i0]);
            forLine.append((* l[1])[i1]);
            forLine.append((* l[2])[i2]);
 
            TMLine & line = * new TMLine(forLine);
            int err = line.fit();
            if (err) {
            delete & line;
            continue;
            }
 
            lines.append(line);
        }
        }
    }
    }
 
#ifdef TRKRECO_DEBUG
    std::cout << "        ... searchLines3 : # of lines found = "
     << lines.length() << std::endl;
#endif
 
    return lines;
}
 
AList<TMLine>
TBuilder::searchLines2(void) const {
    AList<TMLine> lines;
 
    const AList<TMLink> * l[2];
    for (unsigned skip = 0; skip < 15; skip++) {
    if (skip == 0) {
        l[0] = & _links[0];
        l[1] = & _links[1];
    }
    else if (skip == 1) {
        l[1] = & _links[2];
    }
    else if (skip == 2) {
        l[1] = & _links[3];
    }
    else if (skip == 3) {
        l[1] = & _links[4];
    }
    else if (skip == 4) {
        l[1] = & _links[5];
    }
    else if (skip == 5) {
        l[0] = & _links[1];
        l[1] = & _links[2];
    }
    else if (skip == 6) {
        l[1] = & _links[3];
    }
    else if (skip == 7) {
        l[1] = & _links[4];
    }
    else if (skip == 8) {
        l[1] = & _links[5];
    }
    else if (skip == 9) {
        l[0] = & _links[2];
        l[1] = & _links[3];
    }
    else if (skip == 10) {
        l[1] = & _links[4];
    }
    else if (skip == 11) {
        l[1] = & _links[5];
    }
    else if (skip == 12) {
        l[0] = & _links[3];
        l[1] = & _links[4];
    }
    else if (skip == 13) {
        l[1] = & _links[5];
    }
    else if (skip == 14) {
        l[0] = & _links[4];
    }
 
    unsigned n[2];
    for (unsigned i = 0; i < 2; i++)
        n[i] = l[i]->length();
 
    for (unsigned i0 = 0; i0 < n[0]; i0++) {
        for (unsigned i1 = 0; i1 < n[1]; i1++) {
        AList<TMLink> forLine;
        forLine.append((* l[0])[i0]);
        forLine.append((* l[1])[i1]);
 
        TMLine & line = * new TMLine(forLine);
        int err = line.fit();
        if (err) {
            delete & line;
            continue;
        }
 
        lines.append(line);
        }
    }
    }
 
#ifdef TRKRECO_DEBUG
    std::cout << "        ... searchLines2 : # of lines found = "
     << lines.length() << std::endl;
#endif
 
    return lines;
}
 
AList<TMLine>
TBuilder::searchLines1(void) const {
    AList<TMLine> lines;
 
    TMLine & line = * new TMLine(_forLine);
    int err = line.fit();
    if (err) {
    delete & line;
    }
    else {
    lines.append(line);
    }
 
#ifdef TRKRECO_DEBUG
    std::cout << "        ... searchLines1 : # of lines found = "
     << lines.length() << std::endl;
#endif
 
    return lines;
}
 
TMLine
TBuilder::searchLine(const TMLine & initialLine) const {
#ifdef TRKRECO_WINDOW
    std::string old = sz.text();
#endif
 
    static float _szLinkMaxDistance = 5;
 
    unsigned nLinks = _forLine.length();
    TMLine line0 = initialLine;
    line0.remove(initialLine.links());
#ifdef TRKRECO_DEBUG
    line0.fitted(true);
#endif    
    unsigned nGoodLinks0 = 0;
    AList<TMLink> targets;
    unsigned nLoops = 0;
    while (1) {
    ++nLoops;
 
    unsigned nGoodLinks = 0;
    for (unsigned i = 0; i < nLinks; i++) {
        float distance = line0.distance(* _forLine[i]);
        if (distance < _szLinkDistance)
        targets.append(_forLine[i]);
        if (distance < _szLinkMaxDistance) {
        ++nGoodLinks;
        }
    }
 
#ifdef TRKRECO_DEBUG
    std::cout << "    ... searchLine : # of close hits(last) = " << nGoodLinks0
         << ", # of close hits = " << nGoodLinks << std::endl;
#endif
 
    //...Check condition...
    if (nGoodLinks0 == nGoodLinks) break;
 
    //...Do fit...
    TMLine newLine(targets);
    static TRobustLineFitter fitter("robust fitter");
    int err = fitter.fit(newLine);
    if (err) {
#ifdef TRKRECO_WINDOW
        std::cout << "    ... searchLine : failed to fit" << std::endl;
        std::string s = "line search failed";
        sz.text(s);
        sz.wait();
#endif
        break;
    }
 
    //...To protect infinite loop...
    if (nLoops > 10) {
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "        reached to max # of loops(10) : break";
#endif
        break;
    }
 
#ifdef TRKRECO_WINDOW
    std::string ss = ", # of close hits(last) = " + itostring(int(nGoodLinks0))
        + ", # of close hits = " + itostring(int(nGoodLinks));
    sz.append(line0, leda_brown);
    sz.append(newLine, leda_green);
    sz.text(old + ", nLoops = " + itostring(int(nLoops)) + ss);
    sz.wait();
    sz.remove(line0);
    sz.remove(newLine);
#endif
 
    line0 = newLine;
    nGoodLinks0 = nGoodLinks;
    targets.removeAll();
    }
 
    return line0;
}
 
TTrack *
TBuilder::build(TTrack & t, const TMLine & l) const {
    AList<TMLink> links = l.links();
    AList<TMLink> toRemove;
    unsigned n = links.length();
    for (unsigned i = 1; i < n; i++) {
    if (links[i - 1]->link() == links[i]->link()) {
        toRemove.append(links[i]);
        continue;
    }
    if (i < 2)
        continue;
    if (links[i - 2]->link() == links[i]->link())
        toRemove.append(links[i]);
    }
    links.remove(toRemove);
 
    //...Pick up links...
    n = links.length();
    for (unsigned i = 0; i < n; i++)
    t.append(* links[i]->link());
 
    Vector a(5);
    a = t.helix().a();
    a[3] = l.b();
    //a[4] = t.charge() * l.a();//yzhang delete 2012-05-08 
    //    //yzhang add 2012-05-08 
    if(_fitter.getMagneticFieldPointer() == NULL){
      std::cout<<"   "<<__FILE__<<"   "<<__LINE__<<" Could not get MagneticFieldSvc  "<<std::endl;
    }
    const double Bz = -1000 * _fitter.getMagneticFieldPointer()->getReferField();
    if(Bz>0){
      a[4] = t.charge() * l.a();//yzhang add 2012-05-08 ???
    }else{
      a[4] = -1. * t.charge() * l.a();//yzhang add 2012-05-08  ???
    }
    //zhangy
    t._helix->a(a);
 
//#ifdef TRKRECO_WINDOW
//    TTrack * td = new TTrack(t);
//    sz.append(* td, leda_green);
//#endif
#ifdef TRKRECO_DEBUG_DETAIL
    static unsigned nTrk = 0;
    ++nTrk;
//      for (unsigned i = 0; i < t.nCores(); i++)
//      if (t.cores()[i]->wire()->name() == "43=220") {
//          std::cout << "43=220 removed" << std::endl;
//          t.remove(* t.cores()[i]);
//      }
    HepVector3D pp;
    if (links[0]->hit()->mc())
    pp = links[0]->hit()->mc()->hep()->p().vect();
    const HepVector3D p0 = pp;
//      if (nTrk == 1) {
//      Helix tmph(Point3D(0, 0, 0), p0, +1);
//      * t._helix = tmph;
//      }
    t.dump("detail", "before fit");
    std::cout << "Pdif mag=" << (t.p() - p0).mag() << std::endl;
//      t.links().removeAll();
//      for (unsigned i = 0; i < BsCouTab(RECMDC_WIRHIT
//    }
#endif
 
    //...Refine...
    AList<TMLink> bad;
    int err = _fitter.fit(t);
 
#ifdef TRKRECO_DEBUG_DETAIL
    t.dump("detail", "after fit");
    std::cout << "Pdif mag=" << (t.p() - p0).mag() << std::endl;
#endif
 
    t.refine(bad, _maxSigmaStereo * 100.);
    err = _fitter.fit(t);
    t.refine(bad, _maxSigmaStereo * 10.);
    err = _fitter.fit(t);
    t.refine(bad, _maxSigmaStereo);
#ifdef TRKRECO_DEBUG_DETAIL
//    if (nTrk == 2) {
//      Helix tmph(Point3D(0, 0, 0), p0, +1);
//      * t._helix = tmph;
//      std::cout << "initial mom" << p0 << " is set" << std::endl;
//    }
#endif
    err = _fitter.fit(t);
 
#ifdef TRKRECO_DEBUG_DETAIL
    t.dump("detail","         ");
    std::cout << "Pdif mag=" << (t.p() - p0).mag() << std::endl;
#endif
 
    //...Termination...
    return & t;
}
 
bool
TBuilder::initializeForStereo(const TTrack & t,
                  const AList<TSegment> & segments,
                  const AList<TSegment> & badSegments) const {
    _nSuperLayers = 0;
    for (unsigned i = 0; i < 5; i++)
    _nHits[i] = 0;
    //...sz position of segments...
    unsigned nSegments = segments.length();
//  cout<<" used Segments in stereo = "<<nSegments<<endl;
    for (unsigned i = 0; i < nSegments; i++) {
    TMLink & l = * new TMLink();
    int err = t.szPosition(* segments[i], l);
    //...Remove if sz can not be calculcated...
    if (err) {
        delete & l;
        continue;
    }
//  _links[l.wire()->superLayerId() / 2].append(l);
    if (!_links[l.wire()->axialStereoLayerId()/4].hasMember(l))
        _links[l.wire()->axialStereoLayerId()/4].append(l);      
//? _allLinks.append(segments[i]->links());
//  _allLinks.append(segments[i]->cores());
    for (unsigned j = 0; j < segments[i]->cores().length(); ++j) {
        TMLink & cL = * segments[i]->cores()[j];
        if(!_allLinks.hasMember(cL)) _allLinks.append(cL);
    }
    }
 
    //...Count a number of super layers...
    for (unsigned i = 0; i < 6; i++) {
    if (_links[i].length() > 0) {
        ++_nSuperLayers;
    }
//      else {
//          _allLinks.append(badSegments[i]->links());
//  #ifdef TRKRECO_DEBUG_DETAIL
//          std::cout << "    ... stereo super layer " << i * 2 + 1
//           << " has no link,"
//           << badSegments[i]->links().length() << " (bad) links added"
//           << std::endl;
//          badSegments[i]->dump("hits sort flag", "        ");
//  #endif
//      }
    }
 
    //...Append bad links also...
    if (badSegments.length()) {
    for (unsigned i = 0; i < 6; i++) {
        if (badSegments[i]->links().length()) {
        _allLinks.append(badSegments[i]->links());
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "    ... bad links added for stereo super layer "
             << i * 2 + 1 << std::endl;
        badSegments[i]->dump("hits sort flag", "        ");
#endif
        }
    }
    }
 
    //...sz position of links...
    unsigned nCores = _allLinks.length();
    if (nCores < _minNCores) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    ... initializeForStereo : # of cores(=" << nCores
         << ") is less then " << _minNCores << std::endl;
#endif
    return false;
    }
    for (unsigned j = 0; j < nCores; j++) {
    TMLink & link = * _allLinks[j];
    for (unsigned i = 0; i < 2; i++) {
        TMLink & tt = * new TMLink(link);
        unsigned stateR=(link.hit()->state())&WireHitPatternRight;
        unsigned stateL=(link.hit()->state())&WireHitPatternLeft;  
        bool boolR=(stateR==WireHitPatternRight);
        bool boolL=(stateL==WireHitPatternLeft); 
//      std::cout<<"zsl TBuilder::initializeForStereo::boolR:"<<boolR<<", boolL:"<<boolL<<endl; 
 
        tt.leftRight(i);
        int err = t.szPosition(tt);
//      std::cout<<"zsl TBuilder::initializeForStereo::err:"<<err<<endl;
//            std::cout<<"----------------------"<<endl;
        
        if (err) {
        delete & tt;
        continue;
        }
        tt.link(& link);
        _forLine.append(tt);
    }
    }
 
    //...Count # of axial hits...
    unsigned nA = t.cores().length();
    for (unsigned i = 0; i < nA; i++)
//  ++_nHits[t.cores()[i]->wire()->superLayerId() / 2];
        ++_nHits[t.cores()[i]->wire()->axialStereoLayerId()/4];
    
#ifdef TRKRECO_DEBUG
    std::cout << "    ... initializeForStereo : axial super layer usage = ";
    for (unsigned i = 0; i < 5; i++)
    std::cout << _nHits[i] << ",";
    std::cout << std::endl
     << "                            : # of stereo super layers="
     << _nSuperLayers << std::endl;
#endif
 
    return true;
}
 
unsigned
TBuilder::stereoQuality(const AList<TMLink> & links) const {
    unsigned n[6] = {0, 0, 0, 0, 0, 0};
 
    //...Check superlayers...
    unsigned nLinks = links.length();
    for (unsigned i = 0; i < nLinks; i++) {
    const TMLink & l = * links[i];
    if (l.wire()->axial()) continue;
//  ++n[l.wire()->superLayerId() / 2];
//        cout<<l.wire()->axialStereoLayerId()<<endl;   
    ++n[l.wire()->axialStereoLayerId()/4];
    }
 
    unsigned nTotal = 0;
    unsigned nMissing = 0;
    unsigned innermostStereo = 5;
    unsigned outermostStereo = 0;
    unsigned innermostAxial = 4;
    unsigned outermostAxial = 0;
    unsigned innerSL = 6;
    unsigned outerSL = 0;
    for (unsigned i = 0; i < 6; i++) {
    if (_links[i].length() > 1) {
        if(i > outermostStereo) outermostStereo = i;
        if(i < innermostStereo) innermostStereo = i;
    }
    if (n[i] > 1) ++nTotal;
    }
    for (unsigned i = 0; i < 5; i++) {
    if (_nHits[i] > 1) {
        if (i > outermostAxial) outermostAxial = i;
        if (i < innermostAxial) innermostAxial = i;
    }
    }
 
    if(innermostStereo > 1 && innermostAxial < 3) innerSL = 2;
    else innerSL = innermostStereo;
    if(outermostStereo > 1 && outermostAxial > 2) outerSL = 5;
    else outerSL = outermostStereo;
 
    for (unsigned i = 0; i < 6; i++) {
    if (_links[i].length() > 0) {
//zsl 051014        if ((_nHits[i] > 1) && (_nHits[i + 1] > 1))
        if (i <= outerSL && i >= innerSL)
        if (n[i] < 2)
            ++nMissing;
    }
    if (n[i] > 1) ++nTotal;
    }
    unsigned toBeReturn = 0;
    if (nMissing <= 1) toBeReturn = 2;
    else if (nMissing == 2) toBeReturn  = 1;
    if (nTotal < 2) toBeReturn = 0;
 
#ifdef TRKRECO_DEBUG
    std::cout << "    ... stereoQuality : axial ";
    for (unsigned i = 0; i < 6; i++)
    std::cout << _nHits[i] << ",";
    std::cout << std::endl
     << "                      : stereo ";
    for (unsigned i = 0; i < 5; i++)
    std::cout << n[i] << ",";
    std::cout << " : total " << nTotal << " super layers, "
     << " : quality=" << toBeReturn << std::endl;
#endif
 
    return toBeReturn;
}
 
 
TTrack *
TBuilder::buildStereoNew(const TTrack & t,
             AList<TSegment> & segments,
             AList<TSegment> & badSegments) const {
 
#ifdef TRKRECO_DEBUG
    std::cout << "... building stereo by segments(new) : # of segments = ";
    std::cout << segments.length() << std::endl;
    for (unsigned i = 0; i < segments.length(); i++)
    segments[i]->dump("hits sort flag", "    ");
#endif
#ifdef TRKRECO_WINDOW
    sz.clear();
#endif
    TTrack * bestCandidate = NULL;
    AList< AList<TMLink> > poorSeeds;
    bool ok = initializeForStereo(t, segments, badSegments);
    unsigned nSuperLayers = _nSuperLayers + 1;
//        std::cout<<"nSuperLayers in buildStereoNew:"<<_nSuperLayers<<endl;  
//    cout<<"ok of initializeForStereo: "<<ok<<endl;
 
    if (! ok) goto endOfBuilding;
    //...Main loop...
    while (--nSuperLayers) {
 
    //...Initial line search...
    AList<TMLine> initialLines = searchInitialLines(nSuperLayers);
#ifdef TRKRECO_WINDOW
    sz.clear();
    sz.skip(false);
    sz.mode(2);
    sz.appendSz(t, segments, leda_black);
    AList<TSegment> tmps = segments;
    std::string s = "# of segments = " + itostring(int(segments.length()));
    sz.appendSz(t, _allLinks, leda_black);
    sz.append(_forLine, leda_brown);
    s = "nSprLyr=" + itostring(int(nSuperLayers)) +
        ", # of initial lines = " + itostring(int(initialLines.length()));
    for (unsigned i = 0; i < initialLines.length(); i++)
        sz.append(* initialLines[i], leda_red);
    sz.text(s);
    sz.wait();
#endif
//  cout<<"length of initialLines = "<<initialLines.length()<<endl;
    if (initialLines.length() == 0) continue;
    //...Line loop...
    bool found = false;
    unsigned nInitialLines = initialLines.length();
//        cout<<"nInitialLines:"<<nInitialLines<<endl;
    
    for (unsigned i = 0; i < nInitialLines; i++) {
 
        //...Linear fit...
        const TMLine & line = * initialLines[i];
        TMLine newLine = searchLine(line);
 
        //...Skip if this is a seed of the poor result...
        if (poorSeeds.length()) {
        bool poorCase = false;
        for (unsigned j = 0; j < poorSeeds.length(); j++) {
            if (poorSeeds[j]->length() == newLine.links().length()) {
            AList<TMLink> tmp = * poorSeeds[j];
            tmp.remove(newLine.links());
            if (tmp.length() == 0) {
#ifdef TRKRECO_DEBUG_DETAIL
                std::cout << "    ... This is a poor seed :"
                      << " skipped"
                      << " : # of poor seeds = "
                      << poorSeeds.length() << std::endl;
#endif
                poorCase = true;
                break;
            }
            }
        }
        if (poorCase) continue;
        }
        
        //...Is this a good line...
        if (! stereoQuality(newLine.links()))
        continue;
 
        //...3D fit...
        TTrack * t3d = new TTrack(t);
        t3d = build(* t3d, newLine);
        if (t3d == NULL) continue;
 
        //...Check quality...
        unsigned quality = stereoQuality(t3d->links());
//      cout<<"point3, quality:"<<quality<<endl;
 
        //...Best case...
        if (quality == 2) {
#ifdef TRKRECO_WINDOW
        sz.text("stereo finished");
        sz.oneShot(* t3d, leda_blue);
#endif
        if (bestCandidate) delete bestCandidate;
        bestCandidate = t3d;
        found = true;
        break;
        }
        
        //...Poor case...
        AList<TMLink> * tmpL = new AList<TMLink>();
        tmpL->append(newLine.links());
        poorSeeds.append(tmpL);
        if (quality == 0) {
#ifdef TRKRECO_WINDOW
        sz.text("this candidate discarded");
        sz.oneShot(* t3d, leda_black);
#endif
        delete t3d;
        continue;
        }
        
        //...Not enough...
        if (bestCandidate) {
        if (bestCandidate->cores().length() <
            t3d->cores().length()) {
#ifdef TRKRECO_WINDOW
            sz.text("new candidate");
            sz.append(* bestCandidate, leda_brown);
            sz.oneShot(* t3d, leda_green);
            sz.remove(* bestCandidate);
#endif
            delete bestCandidate;
            bestCandidate = t3d;
        }
        else {
#ifdef TRKRECO_WINDOW
            sz.text("this candidate discarded");
            sz.oneShot(* t3d, leda_black);
#endif
            delete t3d;
        }
        }
        else {
        bestCandidate = t3d;
#ifdef TRKRECO_WINDOW
        sz.text("new candidate");
        sz.oneShot(* bestCandidate, leda_green);
        sz.remove(* t3d);
#endif
        }
    }
    
    //...Termination of a loop...
    HepAListDeleteAll(initialLines);
    if (found) break;
    }
 
endOfBuilding:
    _allLinks.removeAll();
    for (unsigned i = 0; i < 6; i++)
    HepAListDeleteAll(_links[i]);
    HepAListDeleteAll(_forLine);
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    ... # of poor seeds = " << poorSeeds.length() << std::endl;
#endif
#ifdef TRKRECO_WINDOW
    if (bestCandidate == NULL) {
    sz.text("3D failed");
    sz.wait();
    }
#endif
        
    //...Check used segments...
    if (bestCandidate) {
    AList<TMLink> usedLinks = StereoHits(bestCandidate->links());
    for (unsigned i = 0; i < segments.length(); i++) {
        TSegment & segment = * segments[i];
        AList<TMLink> used = Links(segment, * bestCandidate);
        if (used.length()) {
        bestCandidate->segments().append(segment);
        segment.tracks().append(bestCandidate);
        }
    }
    }
    
    if (poorSeeds.length())
    HepAListDeleteAll(poorSeeds);
 
#ifdef TRASA_DEBUG_DETAIL
    if (bestCandidate == NULL)
    std::cout << "... building stereo(new) failed" << std::endl;
    else
    std::cout << "... building stereo(new) ok" << std::endl;
#endif
    return bestCandidate;
}