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

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TConformalFinder0.cxx,v 1.14 2010/03/31 09:58:59 liucy Exp $
//-----------------------------------------------------------------------------
// Filename : TConformalFinder0.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class to find tracks with the conformal method.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
// TConformalLink removed. TSegment added.
//-----------------------------------------------------------------------------

#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TConformalFinder0.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/THistogram.h"
#include "TrkReco/TCircle.h"
#include "TrkReco/TTrack.h"
#include "TrkReco/TSegment0.h"
//#include "cdc/Range.h"
#include "TrkReco/Range.h"

std::string
TConformalFinder0::version(void) const {
    return "1.63";
}

TConformalFinder0::TConformalFinder0(float maxSigma,
                                     float fraction,
                                     float stereoZ3,
                                     float stereoZ4,
                                     float stereoChisq3,
                                     float stereoChisq4,
                                     float stereoMaxSigma,
                                     unsigned fittingCorrections,
                                     float salvageLevel,
                                     bool cosmic)
: TFinderBase(),
  _builder(0),
  _doStereo(true),
  _doSalvage(true),  //liuqg
  _fraction(fraction) {

    //...Parameters for a track...
    _trackSelector.nLinks(4);
    _trackSelector.nSuperLayers(2);
    _trackSelector.minPt(0.05);
    _trackSelector.maxImpact(100.);
    _trackSelector.maxSigma(maxSigma);
    _trackSelector.nLinksStereo(3);
    _trackSelector.maxDistance(30.);

    //...Make a builder...
    if (cosmic) _builder = new TBuilderCosmic("cosmic builder", salvageLevel);
    else        _builder = new TBuilder0("conformal builder",
                                         stereoZ3,
                                         stereoZ4,
                                         stereoChisq3,
                                         stereoChisq4,
                                         stereoMaxSigma,
                                         fittingCorrections,
                                         salvageLevel);

    //...Set up TBuilder...
    _builder->trackSelector(_trackSelector);
}

TConformalFinder0::~TConformalFinder0() {
    delete _builder;
}

void
TConformalFinder0::dump(const std::string & msg, const std::string & pre) const {
    std::cout << pre;
    TFinderBase::dump(msg);
    std::cout << pre;
    if (msg.find("state") != std::string::npos) {
        std::cout << "#axialConfPos=" << _axialConfLinks.length();
        std::cout << ",#stereoConfPos=" << _stereoConfLinks.length();
    }
}

void
TConformalFinder0::clear(void) {
    HepAListDeleteAll(_axialConfLinks);
    HepAListDeleteAll(_stereoConfLinks);
    _unusedAxialConfLinks.removeAll();
    _unusedStereoConfLinks.removeAll();
    _goodAxialConfLinks.removeAll();
    HepAListDeleteAll(_circles);
    _tracks.removeAll();
}

void
TConformalFinder0::conformalTransformation(const HepPoint3D & center,
                                          const AList<TMDCWireHit> & hits,
                                          AList<TMLink> & links) {

    unsigned nHits = hits.length();
    if (center == ORIGIN) {
        for (unsigned i = 0; i < nHits; i++) {
            TMDCWireHit * h = hits[i];
            const HepPoint3D & p = h->xyPosition();
            
//zsl       HepPoint3D cp(2. * p.x() / p.mag2(), 2. * p.y() / p.mag2());
            HepPoint3D cp(2. * p.x() / p.mag2(), 2. * p.y() / p.mag2(), 0.);
//          std::cout<<" xypo "<<p<<" pos "<<cp<<std::endl;
            links.append(new TMLink(0, h, cp));
        }
    }
    else {
        for (unsigned i = 0; i < nHits; i++) {
            TMDCWireHit * h = hits[i];
            HepPoint3D p(h->xyPosition() - center);
            HepPoint3D cp(2. * p.x() / p.mag2(), 2. * p.y() / p.mag2(), 0.);
            links.append(new TMLink(0, h, cp));
        }       
    }
}

void
TConformalFinder0::conformalTransformationDriftCircle(const HepPoint3D & center,
                                                      const AList<TMDCWireHit> & hits,
                                                      AList<TMLink> & links) { //added by Liuqg for Tsf
    unsigned nHits = hits.length();
    if (center == ORIGIN) {
        for (unsigned i = 0; i < nHits; i++) {
            TMDCWireHit * h = hits[i];
            const HepPoint3D & p = h->xyPosition();
                                                                                                                             
            const double r = 0.5*(h->drift(0) + h->drift(1));
            HepPoint3D cp(2. * p.x() / p.mag2(), 2. * p.y() / p.mag2(), 0.);
//          HepPoint3D cp2(100 * 2. * p.x() / (p.mag2() - r*r), 100 * 2. * p.y() / (p.mag2() - r*r), 0.);
            HepPoint3D cp2(2. * p.x() / (p.mag2() - r*r), 2. * p.y() / (p.mag2() - r*r), 0.);
//          float cDrift = 100 * 2. * r / (p.mag2() - r*r);
            double cDrift = 2. * r / (p.mag2() - r*r);
                                                                                                                             
            links.append(new TMLink(0, h, cp, cp2, cDrift));
        }
    }
    else {
        for (unsigned i = 0; i < nHits; i++) {
            TMDCWireHit * h = hits[i];
            HepPoint3D p(h->xyPosition() - center);
                                                                                                                             
            const double r = 0.5*(h->drift(0) + h->drift(1));
            HepPoint3D cp(2. * p.x() / p.mag2(), 2. * p.y() / p.mag2(), 0.);
            //unit of the following is km-1, origin is cm.
//          HepPoint3D cp2(100 * 2. * p.x() / (p.mag2() - r*r), 100 * 2. * p.y() / (p.mag2() - r*r), 0.);
            HepPoint3D cp2(2. * p.x() / (p.mag2() - r*r), 2. * p.y() / (p.mag2() - r*r), 0.);
//          float cDrift = 100 * 2. * r / (p.mag2() - r*r);
            double cDrift = 2. * r / (p.mag2() - r*r);
                                                                                                                             
            links.append(new TMLink(0, h, cp, cp2, cDrift));
        }
    }
}

void
TConformalFinder0::conformalTransformationRphi(const HepPoint3D & center,
                                              const AList<TMDCWireHit> & hits,
                                              AList<TMLink> & links) {

    unsigned nHits = hits.length();
    if (center == ORIGIN) {
        for (unsigned i = 0; i < nHits; i++) {
            TMDCWireHit * h = hits[i];
            const HepPoint3D & p = h->xyPosition();
            HepPoint3D cp(2. * p.x() / p.mag2(), 2. * p.y() / p.mag2(), 0.);
            double r = log(cp.mag()) + 4.;
            double phi = atan2(cp.y(), cp.x()) + M_PI;
            HepPoint3D cpt(phi, r, 0.);
            links.append(new TMLink(0, h, cpt));
        }
    }
    else {
        for (unsigned i = 0; i < nHits; i++) {
            TMDCWireHit * h = hits[i];
            HepPoint3D p(h->xyPosition() - center);
            HepPoint3D cp(2. * p.x() / p.mag2(), 2. * p.y() / p.mag2(), 0.);
            double r = log(cp.mag()) + 4.;
            double phi = atan2(cp.y(), cp.x()) + M_PI;
            HepPoint3D cpt(phi, r, 0.);
            links.append(new TMLink(0, h, cpt));
        }       
    }
}

AList<TSegment0>
TConformalFinder0::findClusters(const THistogram & hist) const {

    //...Obtain raw clusters...
    AList<TSegment0> list = hist.clusters0();
    unsigned n = list.length();
    if (n == 0) return list;

#ifdef TRKRECO_DEBUG_DETAIL
    // static TChecker chk0("clusters before splitting");
    // chk0.check(list);
    // chk0.dump("detail", "    ");
#endif

    //...Examine each cluster...
    AList<TSegment0> splitted;
    for (unsigned i = 0; i < n; i++) {
        TSegment0 * c = list[i];

        AList<TSegment0> newClusters = c->split();
        if (newClusters.length() == 0) {
            c->solveDualHits();
            continue;
        }

        list.append(newClusters);
        splitted.append(c);
#ifdef TRKRECO_DEBUG_DETAIL
        c->dump("hits", "    ");
        std::cout << "    ... splitted as" << std::endl;
        for (unsigned j = 0; j < newClusters.length(); j++) {
            std::cout << "    " << j << " : ";
            newClusters[j]->dump("hits");
        }
#endif
    }
    list.remove(splitted);
    HepAListDeleteAll(splitted);

#ifdef TRKRECO_DEBUG_DETAIL
    // static TChecker chk1("clusters after splitting");
    // chk1.check(list);
    // chk1.dump("detail", "    ");
#endif

    return list;
}

AList<TSegment0>
TConformalFinder0::findClusters2(const THistogram & hist) const {

    //...Obtain raw clusters...
    AList<TSegment0> list = hist.clusters0();
    unsigned n = list.length();
    if (n == 0) return list;

#ifdef TRKRECO_DEBUG_DETAIL
    // static TChecker chk0("clusters before splitting (2)");
    // chk0.check(list);
    // chk0.dump("detail", "    ");
#endif

    //...Examine each cluster...
    for (unsigned i = 0; i < n; i++) {
        TSegment0 * c = list[i];
        unsigned type = c->clusterType();

        if ((type == 1) || (type == 2)) {
            c->dump("hits mc", "    ");
            c->solveDualHits();
        }
    }

#ifdef TRKRECO_DEBUG_DETAIL
    // static TChecker chk1("clusters after splitting (2)");
    // chk1.check(list);
    // chk1.dump("detail", "    ");
#endif

    return list;
}

AList<TMLink>
TConformalFinder0::findCloseHits(const AList<TMLink> & links,
                                const TTrack & track) const {
    //
    // Coded by J.Suzuki
    //
    AList<TMLink> list;

    //...Check condition...
    if (track.links().length() == 0) {
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "TConformalFinder0::findCloseHits !!! ";
        std::cout << " no links found in a track : This should not be happened";
        std::cout << std::endl;
#endif

        return list;
    }

    //...Parameters...
    //    float dRcut[11] = {0, 3.5, 0., 5.5, 0., 6.5, 0., 7.5, 0., 9.5, 0.};
//    float dRcut[11] = {0., 4.3, 0., 6.5, 0., 7.5, 0., 8.0, 0., 9.5, 0.};
    float dRcut[11] = {4.3, 6.5, 0., 0., 0., 7.5, 8.0, 9.5, 11.0, 0., 0.};  //Liuqg, Bes

    //...Select Stereo hits associated to the current r-phi curve...
    double R0 = track.helix().curv();
    double xInnerWire = track.links()[0]->wire()->xyPosition().x();
    double yInnerWire = track.links()[0]->wire()->xyPosition().y();
    unsigned nall = links.length();
    for (unsigned j = 0; j < nall; j++) {
        TMLink & t = * links[j];
        const TMDCWire & w = * t.wire();
        HepVector3D X = w.xyPosition() - track.helix().center();
        double Rmag2 = X.mag2();
        double DR = fabs(sqrt(Rmag2) - fabs(R0));
        t.zStatus(-10);
        t.zPair(0);
        if (DR < dRcut[w.superLayerId()] &&
            (xInnerWire*w.xyPosition().x()+yInnerWire*w.xyPosition().y())>0.){
            list.append(t);
        }
    }

    //cout<<"stereo hits: "<<list.length()<<endl;
    return list;
}

TSegment0 *
TConformalFinder0::findBestLink(const TSegment0 & base,
                               const AList<TSegment0> & candidates) const {
    //...Parameters...
    double minAngle = 0.80;
    double maxDistance = 0.3;
//    double minAngle = 0.40;
//    double maxDistance = 0.5;

    //...Candidate loop...
    unsigned n = candidates.length();
    double minDistance = 999.;
    TSegment0 * best = NULL;
    for (unsigned j = 0; j < n; j++) {
        TSegment0 * current = candidates[j];
        if (current->nLinks() < 2) continue;

        float angle = base.direction().dot(current->direction());
#ifdef TRKRECO_DEBUG_DETAIL
        current->dump("vector hits mc", "    ");
        std::cout << "        angle=" << angle;
        if (angle < minAngle) std::cout << std::endl;
#endif
//      cout<<"             angle = "<<angle<<endl;
//      if(angle < 0.3) cout<<" base.direction(): "<<base.direction()
//      <<"  current->direction(): "<<current->direction()<<endl;
        if (angle < minAngle) continue;

        float distance = base.distance(* current);
        if (distance < minDistance) {
            minDistance = distance;
            best = current;
        }
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << ",dist=" << distance << std::endl;
#endif
    }

//    if (minDistance < maxDistance) return best;
    return best;
//    return NULL;  //Liuqg.....
}

TSegment0 *
TConformalFinder0::appendCluster(TTrack & t, AList<TSegment0> & list) const {

    //...Candidate loop...
    unsigned n = list.length();
    TSegment0 * best = NULL;
    unsigned nBest = 0;
    for (unsigned j = 0; j < n; j++) {
        TSegment0 * c = list[j];

        unsigned nOk = t.testByApproach(c->links(),
                                        _trackSelector.maxSigma());
        if (nOk > nBest) {
            nBest = nOk;
            best = c;
        }
    }

    //...Try to append...
    if (best) {
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "    ... appending a cluster" << std::endl;
        best->dump("hits mc", "        ");
#endif
        AList<TMLink> links(best->links());
        t.appendByApproach(links, _trackSelector.maxSigma());
        return best;
    }

    return NULL;
}

AList<TSegment0>
TConformalFinder0::findClusterLink(TSegment0 & base,
                                  const AList<TSegment0> * const list) const {

#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << name() << " ... finding cluster linkage" << std::endl;
    if (base.links().length() == 0)
        std::cout << name() << " !!! base doesn't have any TMLink." << std::endl;
    std::cout << "... base cluster" << std::endl;
    base.dump("cluster hits mc", "  ->");
#endif

    //...Preparation of return value...
    AList<TSegment0> seeds;
    seeds.append(base);

    //...Which super layer?...
//    unsigned outerMost = (base.links())[0]->wire()->superLayerId() / 2;
    unsigned outerMost = (base.links())[0]->wire()->axialStereoLayerId() / 4;

    //...Inner super layer loop...
    int next = outerMost;
    TSegment0 * last = & base;
    while (next) {
        --next;
        const AList<TSegment0> & candidates = list[next];
        if (candidates.length() == 0) continue;
        
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "... clusters in super layer " << next << std::endl;
#endif

        //...Find best match...
        TSegment0 * best = findBestLink(* last, candidates);
        if (best != NULL) {
            seeds.append(best);
            last = best;
#ifdef TRKRECO_DEBUG_DETAIL
            std::cout << "  ->Best is ";
            std::cout << best->position() << " ";
            best->dump("hits mc");
#endif
        }
    }

    return seeds;
}

TTrack *
TConformalFinder0::makeTrack(const AList<TSegment0> & list) const {
    AList<TMLink> links;
    for (unsigned i = 0; i < list.length(); i++) {
        const AList<TMLink> & tmp = list[i]->links();
        unsigned n = tmp.length();
        for (unsigned j = 0; j < n; j++) {
            if (tmp[j]->hit()->track()) continue;
            links.append(tmp[j]);
        }
    }

    TTrack * t = _builder->buildRphi(links);

    return t;
}

AList<TSegment0>
TConformalFinder0::findCloseClusters(const TTrack & track,
                                    const AList<TSegment0> & list,
                                    double maxDistance) const {

    //...Cal. direction of rotation of track...
    double radius = fabs(track.helix().radius());
    const HepPoint3D & center = track.helix().center();
    int rotation =
        (center.cross(track.links()[0]->xyPosition()).z() > 0.) ? 1 : -1;

#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << name() << " ... finding close clusters:maxDistance=";
    std::cout << maxDistance << std::endl;
    std::cout << "    radius,center,rotation=" << radius << ",";
    std::cout << center << "," << rotation << std::endl;
#endif

    //...Cluster loop...
    AList<TSegment0> close;
    unsigned n = list.length();
    for (unsigned i = 0; i < n; i++) {
        TSegment0 & c = * list[i];

        //...Cal. position of cluster in the real plane...
        HepPoint3D position(ORIGIN);
        unsigned m = c.links().length();
        for (unsigned j = 0; j < m; j++) {
            position += c.links()[j]->xyPosition();
        }
        position *= 1. / double(m);

#ifdef TRKRECO_DEBUG_DETAIL
        c.dump("cluster hits mc", "    ");
        std::cout << "        position=" << position;
        std::cout << ",diff=" << (position - center).mag() - radius << std::endl;
#endif
        //...Cal. distance to a track...
        HepVector3D diff = position - center;
        if ((diff.mag() - radius) < maxDistance) {

            //...Same side?...
            int direction = (center.cross(position).z() > 0.) ? 1 : -1;
            if (direction == rotation) close.append(c);
        }
    }

#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    found clusters" << std::endl;
    for (unsigned i = 0; i < close.length(); i++) {
        close[i]->dump("hits mc", "    ");
    }
#endif
    return close;
}

void
TConformalFinder0::appendClusters2(TTrack & track,
                                 AList<TSegment0> & list) const {

#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << name() << " ... appending clusters remained" << std::endl;
    std::cout << "    clusters to be tested : " << std::endl;
    for (unsigned i = 0; i < list.length(); i++) {
        list[i]->dump("cluster hits mc", "        ");
    }
#endif

    unsigned n = list.length();
    if (n == 0) return;

    AList<TSegment0> closer;
    closer.append(findCloseClusters(track, list, 1.));

#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    found clusters" << std::endl;
    for (unsigned i = 0; i < closer.length(); i++) {
        closer[i]->dump("cluster hits mc", "        ");
    }
#endif
    n = closer.length();
    if (closer.length() == 0) return;

    //...Append them...
    AList<TMLink> candidates;
    for (unsigned i = 0; i < n; i++)
        candidates.append(closer[i]->links());
    _builder->appendClusters(track, candidates);

    //...Remove TMLinks from clusters...
    for (unsigned i = 0; i < n; i++) {
        closer[i]->TTrackBase::remove(track.links());
        if (closer[i]->nLinks() == 0) list.remove(closer[i]);
    }
}

int
TConformalFinder0::doit(const AList<TMDCWireHit> & axialHits,
                        const AList<TMDCWireHit> & stereoHits,
                        AList<TTrack> & tracks,
                        AList<TTrack> &) {

    //...For debug...
    if (debugLevel()) {
        std::cout << name() << " ... processing" << std::endl;
        std::cout << "    axialHits=" << axialHits.length();
        std::cout << ",stereoHits=" << stereoHits.length();
        std::cout << ",tracks=" << tracks.length();
        std::cout << std::endl;

        if (debugLevel() > 1)
            std::cout << name() << " ... conformal transformation0" << std::endl;
    }

    //...Conformal transformation with IP constraint...
    conformalTransformationRphi(ORIGIN, axialHits, _axialConfLinks);
    conformalTransformationRphi(ORIGIN, stereoHits, _stereoConfLinks);
    _unusedAxialConfLinks.append(_axialConfLinks);
    _unusedStereoConfLinks.append(_stereoConfLinks);
    AList<TMLink> unusedConfLinks;
    if (_doSalvage) {
        unusedConfLinks.append(_axialConfLinks);
        unusedConfLinks.append(_stereoConfLinks);
    }

    //...For debug...
    if (debugLevel() > 1)
        std::cout << name() << " ... selecting good hits" << std::endl;

    //...Select good axial hits...
    AList<TMLink> goodHits;
    int nLinks = _axialConfLinks.length();
    for (unsigned i = 0; i < nLinks; i++) {
        TMLink * l = _axialConfLinks[i];
        const TMDCWireHit & h = * l->hit();
//liuqg if ((h.state() & WireHitIsolated) &&
//          (h.state() & WireHitContinuous))
            goodHits.append(l);
    }
    //...Main algorithm...
    standardFinding(goodHits, unusedConfLinks, _fraction);

    //...Main algorithm for second trial...
    specialFinding(goodHits, unusedConfLinks, _fraction);

    //...For debug...
    if (debugLevel()) {
        std::cout << name() << " ... processed : ";
        std::cout << "good hits=" << goodHits.length();
        std::cout << ",tracks=" << _tracks.length();
        std::cout << std::endl;
    }

    tracks.append(_tracks);
    return 0;
}

void
TConformalFinder0::standardFinding(AList<TMLink> & list,
                                  AList<TMLink> & unusedLinks,
                                  double fraction) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << name() << " ... standard finding with salvage : given hits :";
    Dump(list, "sort");
#endif

    //...Find segments...
    AList< AList<TSegment0> > segments = findSegments(list);
    AList<TSegment0> segmentList[5];
    AList<TSegment0> original[5];
    for (unsigned i = 0; i < 5; i++) {
        segmentList[i] = * segments[i];
        original[i] = * segments[i];
//      cout<< " segment: " << i << " : " << segmentList[i].length() << endl;
//      for (unsigned j = 0; j < segmentList[i].length(); j++)
//        cout<<"links in seg: "<<segmentList[i][j]->nLinks()<<endl;
    }

#ifdef TRKRECO_DEBUG_DETAIL
    for (unsigned i = 0; i < 5; i++) {
        std::cout << "... clusters in super layer " << i << std::endl;
        for (unsigned j = 0; j < segmentList[i].length(); j++) {
            segmentList[i][j]->dump("", "    ");
            segmentList[i][j]->dump("hits mc", "      ");
        }
    }
#endif

    //...Main loop...
    AList<TSegment0> retryList;
    AList<TTrack> salvageList;
    unsigned outerMost = 4;
    while (outerMost) {

        while (TSegment0 * base = segmentList[outerMost][0]) {

            //...Get linked clusters...
            AList<TSegment0> clusters = findClusterLink(* base, segmentList);

            //...Make a track...
            TTrack * t = makeTrack(clusters); // don't change 'clusters'
            if (t == NULL) {
                retryList.append(base);
                segmentList[outerMost].remove(base);
                continue;
            }

            //...Check track quality...
            // double f = float(t->nLinks()) / float(nTLinks(clusters));
            double f = float(t->nCores()) / float(NCoreLinks(clusters));
            if (f < fraction) {
#ifdef TRKRECO_DEBUG_DETAIL
                std::cout << "... fraction too low:" << f << std::endl;
                std::cout << "    used cores=" << t->nCores();
                std::cout << ", candidate cores=" << NCoreLinks(clusters);
                std::cout << " ... retry later" << std::endl;
#endif
                retryList.append(base);
                segmentList[outerMost].remove(base);
                delete t;
                continue;
            }

            //...Append other hits...
            appendClusters2(* t, retryList);

#ifdef TRKRECO_DEBUG_DETAIL
            std::cout << name() << " ... 2D result :" << std::endl;
            t->dump("detail", "    ");
#endif

            //...Make it 3D...
            TTrack * ts = t;
            if (_doStereo) {
                ts = _builder->buildStereo(* t,
                          findCloseHits(_unusedStereoConfLinks, * t));
            }

            //...Check track quality...
            if (! ts) {
#ifdef TRKRECO_DEBUG_DETAIL
                std::cout << "... failed to make a track 3D" << std::endl;
#endif
                retryList.append(base);
                segmentList[outerMost].remove(base);
                delete t;
                continue;
            }

            //...Salvaging...
//          _builder->salvage(* t, unusedLinks);

            //...OK...
            t->assign(WireHitConformalFinder);
            t->finder(TrackOldConformalFinder);
//                    TrackOldConformalFinder | TrackValid | Track3D);
            _tracks.append(t);

            //...Remove used links...
            const AList<TMLink> & usedLinks = t->links();
            list.remove(usedLinks);
            unusedLinks.remove(usedLinks);
            _unusedStereoConfLinks.remove(usedLinks);
            for (unsigned i = 0; i <= outerMost; i++)
                segmentList[i].remove(clusters);

            //...For debug...
            if (debugLevel() > 1) {
                std::cout << name() << " ... track # " << _tracks.length() - 1;
                std::cout << " found" << std::endl;
                t->dump("detail", "    ");
            }
        }

        //...Loop termination...
        --outerMost;
    }

    //...Termination...
    for (unsigned i = 0; i < 5; i++) {
        HepAListDeleteAll(original[i]);
        delete segments[i];
    }
}

void
TConformalFinder0::specialFinding(AList<TMLink> & list,
                                 AList<TMLink> & unusedLinks,
                                 double fraction) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << name() << " ... standard finding with salvage : given hits :";
    Dump(list, "sort");
#endif

    //...Find segments...
    AList< AList<TSegment0> > segments = findSegments(list);
    AList<TSegment0> segmentList[5];
    AList<TSegment0> original[5];
    for (unsigned i = 0; i < 5; i++) {
        segmentList[i] = * segments[i];
        original[i] = * segments[i];
    }

#ifdef TRKRECO_DEBUG_DETAIL
    for (unsigned i = 0; i < 5; i++) {
        std::cout << "... clusters in super layer " << i << std::endl;
        for (unsigned j = 0; j < segmentList[i].length(); j++) {
            segmentList[i][j]->dump("", "    ");
            segmentList[i][j]->dump("hits mc", "      ");
        }
    }
#endif

    //...Main loop...
    AList<TSegment0> retryList;
    unsigned outerMost = 4;
    while (outerMost) {

        while (TSegment0 * base = segmentList[outerMost][0]) {

            //...Get linked clusters...
            AList<TSegment0> clusters = findClusterLink(* base, segmentList);

        again:;

            //...Check # of clusters...
            if (clusters.length() < 2) {
                segmentList[outerMost].remove(base);
                continue;
            }

            //...Make a track...
            TTrack * t = makeTrack(clusters); // don't change 'clusters'
            if (t == NULL) {
                clusters.remove(clusters.last());
                goto again;
            }

            //...Check track quality...
            // double f = float(t->nLinks()) / float(nTLinks(clusters));
            double f = float(t->nCores()) / float(NCoreLinks(clusters));
            if (f < fraction) {
#ifdef TRKRECO_DEBUG_DETAIL
                std::cout << "... fraction too low:" << f << std::endl;
                std::cout << "    retry later" << std::endl;
#endif
                delete t;
                clusters.remove(clusters.last());
                goto again;
            }

            //...Append other hits...
            appendClusters2(* t, retryList);

            //...Make it 3D...
            
            TTrack * ts = t;
            if (_doStereo) {
                ts = _builder->buildStereo(* t,
                          findCloseHits(_unusedStereoConfLinks, * t));
            }

            //...Check track quality...
            if (! ts) {
                clusters.remove(clusters.last());
                delete t;
                goto again;
            }

            //...Salvaging...
//          _builder->salvage(* t, unusedLinks);

            //...OK...
            t->assign(WireHitConformalFinder);
            t->finder(TrackOldConformalFinder);
//                    TrackOldConformalFinder | TrackValid | Track3D);
            _tracks.append(t);

            //...Remove used links...
            const AList<TMLink> & usedLinks = t->links();
            list.remove(usedLinks);
            unusedLinks.remove(usedLinks);
            _unusedStereoConfLinks.remove(usedLinks);
            for (unsigned i = 0; i <= outerMost; i++)
                segmentList[i].remove(clusters);

            //...For debug...
            if (debugLevel() > 1) {
                std::cout << name() << " ... track # " << _tracks.length() - 1;
                std::cout << " found" << std::endl;
                t->dump("detail", "    ");
            }
        }

        //...Loop termination...
        --outerMost;
    }

    //...Termination...
    for (unsigned i = 0; i < 5; i++) {
        HepAListDeleteAll(original[i]);
        delete segments[i];
    }
}

AList< AList<TSegment0> >
TConformalFinder0::findSegments(const AList<TMLink> & in) const {
    AList< AList<TSegment0> > a;

#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << name() << " ... finding segments : given hits =" << std::endl;
    Dump(in, "sort");
#endif

    //...Create lists of links for each super layer...
    AList<TMLink> links[5];
    unsigned n = in.length();
    for (unsigned i = 0; i < n; i++) {
        TMLink & l = * in[i];
//      links[l.wire()->superLayerId() / 2].append(l);
        links[l.wire()->axialStereoLayerId()/4].append(l);
    }

    //...Create phi hists and clusters for each super layer...
    THistogram * hist[5];
    hist[0] = new THistogram(76);
    hist[1] = new THistogram(100);
    hist[2] = new THistogram(128);
    hist[3] = new THistogram(256);
    hist[4] = new THistogram(288);
    for (unsigned i = 0; i < 5; i++) {
        hist[i]->fillX(links[i]);
        AList<TSegment0> * b = new AList<TSegment0>();
        a.append(b);
        b->append(findClusters(* hist[i]));
        delete hist[i];
    }

    return a;
}

---