TConformalFinder.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TConformalFinder.cxx,v 1.47 2012/05/28 05:16:29 maoh Exp $
//-----------------------------------------------------------------------------
// Filename : TConformalFinder.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : [email protected]
//-----------------------------------------------------------------------------
// Description : A class to find tracks with the conformal method.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
// TConformalLink removed. TSegment added.
//-----------------------------------------------------------------------------
 
#include <algorithm>
#include <string.h>
#include "CLHEP/String/Strings.h"
//#include "basf/basfshm.h"
#include "TrkReco/TMDC.h"
#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TTrackManager.h"
#include "TrkReco/TConformalFinder.h"
#include "TrkReco/TConformalFinder0.h"
#include "TrkReco/THistogram.h"
#include "TrkReco/TMDCTsf.h"
#include "TrkReco/TCircle.h"
#include "TrkReco/TTrack.h"
#include "TrkReco/TMLine.h"
#include "TrkReco/TTrackHEP.h"
#include "TrkReco/TMDCWireHitMC.h"
 
#ifdef TRKRECO_DEBUG
unsigned TConformalFinder::_stage = ConformalOutside;
#endif
#ifdef TRKRECO_WINDOW
#include "TrkReco/TWindow.h"
#endif
 
//extern BasfSharedMem * BASF_Sharedmem;
 
static const float PI2 = 2. * M_PI;
static const float WIDTH[] = {PI2 / 40,
                  PI2 / 64,           //the first two are stereo, tighter would be better for noise? Liu
                  PI2 / 76,
                  PI2 / 100,
                  PI2 / 128,
                  PI2 / 160,
                  PI2 / 176,
                  PI2 / 208,
                  PI2 / 240,
                  PI2 / 256,
                  PI2 / 288};
 
std::string
TConformalFinder::version(void) const {
    return "3.03";
}
 
TConformalFinder::TConformalFinder(unsigned fastFinder,
                   unsigned slowFinder,
                   unsigned perfectSegmentFinding,
                   float maxSigma,
                   float maxSigmaStereo,
                   float salvageLevel,
                   unsigned minNLinksForSegment,
                   unsigned minNCoreLinks,
                   unsigned minNSegments,
                   unsigned salvageLoadWidth,
                   unsigned stereoMode,
                   unsigned stereoLoadWidth,
                   float szSegmentDistance,
                   float szLinkDistance,
                   unsigned fittingFlag)
: _fastFinder(fastFinder),
  _slowFinder(slowFinder),
  _doT0Reset(false),
  _perfectSegmentFinding(perfectSegmentFinding),
  _segmentSeparation(4),
  _minNLinksForSegment(minNLinksForSegment),
  _minNLinksForSegmentInRefine(2),   //3, 20060307
  _maxNLinksForSegment(8),
  _maxWidthForSegment(4),
  _minNCoreLinks(minNCoreLinks),
  _minNSegments(minNSegments),
  _minUsedFractionSlow2D(0.5),
  _salvageLoadWidth(salvageLoadWidth),
  _stereoMode(stereoMode),
  _stereoLoadWidth(stereoLoadWidth),
  _appendLoad(2),
  _maxSigma2(maxSigma),
  _T0ResetDone(false),
  _builder("conformal builder",
       maxSigma,
       maxSigmaStereo,
       salvageLevel,
       szSegmentDistance,
       szLinkDistance,
       fittingFlag),
  _s(0)
#ifdef TRKRECO_WINDOW
  , _rphiWindow("rphi window")
#endif
{
    _linkMaxDistance[0] = 0.02;
    _linkMaxDistance[1] = 0.025;
    _linkMaxDistance[2] = 0.025;
    _linkMinDirAngle[0] = 0.97;   //0.98
    _linkMinDirAngle[1] = 0.95;   //0.97
    _linkMinDirAngle[2] = 0.95;   //0.97
    //    BASF_Sharedmem->allocate("TrkConfSum", sizeof(struct summary));
    
    //yzhang add 2012-05-03 
    StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF);
    if(scmgn!=StatusCode::SUCCESS) { 
      std::cout<< __FILE__<<"Unable to open Magnetic field service"<<std::endl;
    }
}
 
TConformalFinder::~TConformalFinder() {
}
 
void
TConformalFinder::dump(const std::string & msg, const std::string & pre) const {
    if (msg.find("state") != std::string::npos) {
    std::cout << pre;
    TFinderBase::dump(msg);
    std::cout << pre;
    std::cout << "#axial=" << _hits[0].length();
    std::cout << ",#stereo=" << _hits[1].length();
    }
    if (msg.find("summary") != std::string::npos || msg.find("detail") != std::string::npos) {
/*  struct summary s;
    //      bzero((char *) & s, sizeof(struct summary));
    memset((char *) & s, 0, sizeof(struct summary));
    for (int i = 0; i < BASF_Sharedmem->nprocess(); i++) {
        int size;
        struct summary & r = * (struct summary *)
        BASF_Sharedmem->get_pointer(i, "TrkConfSum", & size);
        s._nEvents += r._nEvents;
        s._nTracksFast3D += r._nTracksFast3D;
        s._nTracksSlow3D += r._nTracksSlow3D;
        s._nTracksFast2D += r._nTracksFast2D;
        s._nTracksSlow2D += r._nTracksSlow2D;
        s._nTracksFast2DBadQuality += r._nTracksFast2DBadQuality;
        s._nTracksSlow2DBadQuality += r._nTracksSlow2DBadQuality;
    }
 
    std::cout << pre
          << "all events : " << s._nEvents << std::endl;
    std::cout << pre
          << "fast 3D tracks              : "
          << s._nTracksFast3D << std::endl;
    std::cout << pre
          << "fast 2D tracks(line failed) : "
          << s._nTracksFast2D << std::endl;
    std::cout << pre
          << "fast 2D tracks(3D failed)   : "
          << s._nTracksFast2DBadQuality << std::endl;
    std::cout << pre
          << "slow 3D tracks              : "
          << s._nTracksSlow3D << std::endl;
    std::cout << pre
          << "slow 2D tracks(line failed) : "
          << s._nTracksSlow2D << std::endl;
    std::cout << pre
          << "slow 2D tracks(3D failed)   : "
          << s._nTracksSlow2DBadQuality << std::endl;
*/
   }
}
 
void
TConformalFinder::clear(void) {
    for (unsigned i = 0; i < 3; i++) {
    if (i == 2)
        HepAListDeleteAll(_allHits[i]);
    else
        _allHits[i].removeAll();
    _hits[i].removeAll();
    _unused[i].removeAll();
    }
    for (unsigned i = 0; i < 2; i++) {
    for (unsigned j = 0; j < 6; j++) {
#ifdef TRKRECO_DEBUG_DETAIL
    cout<<"_allSegments length = "<<_allSegments[i][j].length()<<"  _allUnused length = "<<_allUnused[i][j].length()<<endl;
#endif
        HepAListDeleteAll(_allSegments[i][j]);
        _allUnused[i][j].removeAll();
    }
    }
    _2DTracks.removeAll();
    _3DTracks.removeAll();
}
 
void
TConformalFinder::selectGoodHits(void) {
    for (unsigned i = 0; i < 2; i++) {
    unsigned n = _allHits[i].length();
    for (unsigned j = 0; j < n; j++) {
        unsigned state = _allHits[i][j]->hit()->state();
//      if (_allHits[i][j]->hit()->reccdc()->tdc>800) continue;
//      if ((state & WireHitIsolated) && (state & WireHitContinuous))  //tmp...
//      if ((state & WireHitIsolated))
               _hits[i].append(_allHits[i][j]);
//      else
        _unused[i].append(_allHits[i][j]);
    }
    }
    _hits[2].append(_hits[0]);
    _hits[2].append(_hits[1]);
    _unused[2].append(_unused[0]);
    _unused[2].append(_unused[1]);
 
#ifdef TRKRECO_WINDOW
    _rphiWindow.clear();
    _rphiWindow.skip(false);
    _rphiWindow.skipAllWindow(false);
    _rphiWindow.append(_allHits[2]);
    _rphiWindow.append(_hits[2], leda_pink);
//    _rphiWindow.skip(true);
//    _rphiWindow.skipAllWindow(true);
#endif
}
 
void
TConformalFinder::findSegments(void) {
 
    //...Create lists of links for each super layer...
    AList<TMLink> links[11];
    unsigned n = _hits[2].length();
    for (unsigned i = 0; i < n; i++) {
    TMLink & l = * _hits[2][i];
    links[l.wire()->superLayerId()].append(l);
    }
 
    //...Create phi hists...
    THistogram * hist[11];
/*    hist[0] = new THistogram(64);
    hist[1] = new THistogram(80);
    hist[2] = new THistogram(96);
    hist[3] = new THistogram(128);
    hist[4] = new THistogram(144);
    hist[5] = new THistogram(160);
    hist[6] = new THistogram(176);
    hist[7] = new THistogram(208);
    hist[8] = new THistogram(240);
    hist[9] = new THistogram(256);
    hist[10] = new THistogram(288);
*/
    hist[0] = new THistogram(56);
    hist[1] = new THistogram(80);
    hist[2] = new THistogram(88);
    hist[3] = new THistogram(112);
    hist[4] = new THistogram(140);
    hist[5] = new THistogram(160);
    hist[6] = new THistogram(176);
    hist[7] = new THistogram(208);
    hist[8] = new THistogram(240);
    hist[9] = new THistogram(256);
    hist[10] = new THistogram(288);
    
    for (unsigned i = 0; i < 11; i++) {
//        unsigned superLayerId = i / 2;
//  unsigned id = i % 2;
    unsigned superLayerId;
    unsigned id;
        switch(i){
        case 0: superLayerId = 0; id = 1; break;
        case 1: superLayerId = 1; id = 1; break;
        case 2: superLayerId = 0; id = 0; break;
        case 3: superLayerId = 1; id = 0; break;
        case 4: superLayerId = 2; id = 0; break;
        case 5: superLayerId = 2; id = 1; break;
        case 6: superLayerId = 3; id = 1; break;
        case 7: superLayerId = 4; id = 1; break;
        case 8: superLayerId = 5; id = 1; break;
        case 9: superLayerId = 3; id = 0; break;
        case 10: superLayerId = 4; id = 0; break;
        default: break;         
        }
    // hist[i]->fillX(links[i]);
    hist[i]->fillPhi(links[i]);
 
    //...Segment finding...
    AList<TSegment> tmp = hist[i]->segments();
 
    //...Remove bad segments...
    unsigned n = tmp.length();
    if (id == 0) {
        AList<TSegment> bad;
        for (unsigned j = 0; j < n; j++) {
        TSegment * s = tmp[j];
        if (s->links().length() < _minNLinksForSegment)
            bad.append(s);
        else if (s->links().length() > _maxNLinksForSegment)
            bad.append(s);
//      else if (Width(s->links()) > _maxWidthForSegment)
//          bad.append(s);
        }
        tmp.remove(bad);
        for (unsigned j = 0; j < bad.length(); j++) {
        _unused[id].append(bad[j]->links());
        _unused[2].append(bad[j]->links());
        }
        HepAListDeleteAll(bad);
        n = tmp.length();
    }
 
    //...Classify segments...
//      if (n > 1) {
//          for (unsigned k = 0; k < (n - 1); k++) {
//          TSegment & s0 = * tmp[k];
//          bool ok = true;
//          for (unsigned l = k + 1; l < n; l++) {
//              TSegment & s1 = * tmp[l];
//              int distance =
//              abs(InnerMost(s0.links())->wire()->localIdDifference(
//                  * InnerMost(s1.links())->wire()));
//              if (distance < _segmentSeparation) {
//              s0.state(s0.state() | TSegmentCrowd);
//              s1.state(s1.state() | TSegmentCrowd);
//              }
//          } 
//          }
//      }
 
    //...Store them...
    _allSegments[id][superLayerId].append(tmp);
    _allUnused[id][superLayerId] = _allSegments[id][superLayerId];
    delete hist[i];
    }    
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... segment finding finished" << std::endl;
#endif
}
 
void
TConformalFinder::linkSegments(unsigned level) {
    //...Clear old links...
    unsigned superLayer = 5;
    while (--superLayer) {
    AList<TSegment> & segments = _allUnused[0][superLayer];
    unsigned n = segments.length();
    //cout<<"SL: "<<superLayer<<"; segments.length:"<<n<<endl;
    for (unsigned i = 0; i < n; i++) {
        TSegment & base = * segments[i];
#ifdef TRKRECO_DEBUG
        if (base.tracks().length()) {
        std::cout << "TConformalFinder::linkSegments !!! segment has ";
        std::cout << "an assignment to track(s)" << std::endl;
        }
#endif
        base.innerLinks().removeAll();
        base.outerLinks().removeAll();
    }
    }
 
    //...Super layer loop...
    superLayer = 5;
    while (--superLayer) {
    AList<TSegment> & segments = _allUnused[0][superLayer];
    AList<TSegment> & targets = _allUnused[0][superLayer - 1];
    AList<TSegment> & targets2 = _allUnused[0][superLayer - 2];
    unsigned n = segments.length();
    for (unsigned i = 0; i < n; i++) {
        TSegment & base = * segments[i];
#ifdef TRKRECO_DEBUG
        if (base.tracks().length()) {
        std::cout << "TConformalFinder::linkSegments !!! segment has ";
        std::cout << "an assignment to track(s)" << std::endl;
        }
#endif
 
        const HepPoint3D & p = base.position();
        HepVector3D v = base.direction();
//      std::cout<<"superlayer "<<superLayer<<" link length "<<base.outerLinks().length()<<std::endl;
        if (base.outerLinks().length() == 1)
        v = p - OuterMostUniqueLink(base)->position();
//          if (base.outerLinks().length() == 1)
//          v = p - base.outerLinks()[0]->position();
//          if (base.outerLinks().length() == 1)
//          v = direction(base);
 
        AList<TSegment> & candidates = base.innerLinks();
        TSegment * best = link(base, p, v, targets, candidates, level);
        if ((best == NULL) && (level > 0) && (superLayer > 1)) {
        best = link(base, p, v, targets2, candidates, level);
        }
        if (best) candidates.insert(best);
//      base.dump("hits");
//      if(best)best->dump("hits");
 
        unsigned m = candidates.length();
        for (unsigned j = 0; j < m; j++)
        candidates[j]->outerLinks().append(base);
    }
    }
 
    //check innerLinks. Liuqg
/*    superLayer = 5;
    while(--superLayer){
      //cout<<"SuperLayer: "<<superLayer<<endl;
      AList<TSegment> & segments = _allUnused[0][superLayer];
      unsigned n = segments.length();
      for (unsigned k = 0; k < n; k++) {
    TSegment & base = * segments[k];
    int inners = base.innerLinks().length();
    if (inners<2) continue;
    AList<TSegment> badList;
    for(int i = 0; i < inners; ++i){
      for(int j = i+1; j < inners; ++j){
        AList<TMLink> links = base.innerLinks()[j]->links();
        AList<TMLink> multiLinks;
        for (int k = 0; k < links.length(); ++k){
        TMLink & l = * links[k];
        if (base.innerLinks()[i]->links().hasMember(l)) multiLinks.append(l);
        }
        if (multiLinks.length() < 2) continue;
        int nHits[4];
        int multiLayers = 0;
        for (int k = 0; k < multiLinks.length(); ++k)
        ++nHits[multiLinks[k]->hit()->wire()->localLayerId()];
        for (int k = 0; k < 4; ++k)
        if (nHits[k] > 0) ++multiLayers;
        if(multiLayers >= 2 && (links.length() > base.innerLinks()[i]->links().length())) 
          badList.append(base.innerLinks()[i]);
        else if (multiLayers >= 2) badList.append( base.innerLinks()[j]);
        else continue;
      }
    }
    base.innerLinks().remove(badList);
      }
    }*/
 
#ifdef TRKRECO_DEBUG
    superLayer = 5;
    while(--superLayer){
      cout<<"SuperLayer: "<<superLayer<<endl;
      AList<TSegment> & segments = _allUnused[0][superLayer];
      unsigned n = segments.length();
      for (unsigned i = 0; i < n; i++) {
    TSegment & base = * segments[i];
    cout<<"Layer: "<<base.links()[0]->hit()->wire()->layerId()
      <<"   Local: "<<base.links()[0]->hit()->wire()->localId()
      <<"   innerSeg: "<<base.innerLinks().length()<<endl;
      }
    }
#endif
 
#ifdef TRKRECO_WINDOW
    // _rphiWindow.skipAllWindow(false);
    displayStatus("Conf::linkSegments ... results");
    _rphiWindow.wait();
#endif    
}
 
void
TConformalFinder::resolveSegments(AList<TTrack> & trackCandidates) const {
#ifdef TRKRECO_DEBUG
    std::cout << "... resolving assignments" << std::endl;
#endif
 
    //...Pick up segments which has multiple assignments...
    AList<TSegment> multi;
    unsigned n = trackCandidates.length();
    for (unsigned i = 0; i < n; i++) {
    TTrack & t = * trackCandidates[i];
    AList<TSegment> & segments = t.segments();
    unsigned nS = segments.length();
    for (unsigned j = 0; j < nS; j++) {
        if (segments[j]->tracks().length() > 1)
        multi.append(segments[j]);
    }
    }
    multi.purge();
 
    //...Resolve assignments...
    n = multi.length();
    for (unsigned i = 0; i < n; i++) {
    TSegment & s = * multi[i];
    const AList<TTrack> & tracks = s.tracks();
    unsigned nT = tracks.length();
 
    //...Check TMLink overlap...
    AList<TMLink> multiLinks;
    const AList<TMLink> & links = s.links();
    unsigned nL = links.length();
    for (unsigned j = 0; j < nL; j++) {
        TMLink & l = * links[i];
        bool overlap = false;
        for (unsigned k = 0; k < nT; k++) {
        TTrack & t = * tracks[k];
        if (t.links().hasMember(l))
            overlap = true;
        }
        multiLinks.append(l);
    }
    multiLinks.purge();
 
#ifdef TRKRECO_DEBUG
    std::cout << "    segment : ";
    s.dump("hits sort flag");
    std::cout << "        # of assigned tracks = " << nT << std::endl;
    std::cout << "        # of overlap TMLinks = " << multiLinks.length();
    std::cout << std::endl;
#endif
    //...Select the closest...
    nL = multiLinks.length();
    for (unsigned j = 0; j < nL; j++) {
        TMLink & l = * multiLinks[j];
 
        float bestDiff = 999999999.;
        TTrack * best = NULL;
        for (unsigned k = 0; k < nT; k++) {
        TTrack & t = * tracks[k];
        t.approach(l);
        float distance = (l.positionOnWire() - l.positionOnTrack())
            .mag();
        float diff = fabs(distance - l.hit()->drift());
        if (diff < bestDiff) {
            bestDiff = diff;
            best = & t;
        }
        }
 
        for (unsigned k = 0; k < nT; k++) {
        TTrack & t = * tracks[k];
        if (& t == best) continue;
        t.remove(l);
        }       
 
#ifdef TRKRECO_DEBUG
        {
        std::cout << "            " << l.wire()->name() << " -> ";
        std::cout << best->name() << std::endl;
        }
#endif
    }
    }
}
 
AList<TSegment>
TConformalFinder::removeBadSegments(TTrack & t) const {
#ifdef TRKRECO_DEBUG
    std::cout << "... removing bad segments : # used seg = ";
    std::cout << t.segments().length() << std::endl;
#endif
#ifdef TRKRECO_DEBUG_DETAIL
    for (unsigned i = 0; i < t.segments().length(); i++)
    t.segments()[i]->dump("hits sort flag", "    ");
#endif
 
    const AList<TSegment> & segments = t.segments();
    AList<TSegment> bads;
    unsigned used = 0;
    TSegment * innerMost;
    unsigned n = segments.length();
    for (unsigned i = 0; i < n; i++) {
    TSegment & s = * segments[i];
    AList<TMLink> links = Links(s, t);
    unsigned nLinks = links.length();
 
    unsigned nCores = Cores(links).length();
    unsigned nCoresSegment = s.nCores();
//      if (nCores < _minNCoreLinks) {
//  if (nCores < (nCoresSegment / 2)) {
    if ((nCores < _minNCoreLinks) && (nCores < ((nCoresSegment + 1) / 2))){
        bads.append(s);
        continue;
    }
 
    unsigned id = segments[i]->superLayerId();
    used |= (1 << id);
    if (! id) innerMost = & s;
    }
 
//      //...Check used super layers...
//      n = segments.length();
//      for (unsigned i = 0; i < n; i++)
//      if ((used & 0x155) == 0x101)
//          bads.append(innerMost);
 
    if (! bads.length()) return bads;
 
#ifdef TRKRECO_DEBUG
    std::cout << "    bad segments : # = " << bads.length() << std::endl;
#endif
 
    n = bads.length();
    for (unsigned i = 0; i < n; i++) {
    TSegment & s = * bads[i];
 
#ifdef TRKRECO_DEBUG
    AList<TMLink> links = Links(s, t);
    unsigned nLinks = links.length();
    unsigned nCores = Cores(links).length();
    std::cout << "    # used links = " << nLinks;
    std::cout << ", # used cores = " << nCores << std::endl;
    s.dump("hits sort flag", "        ");
#endif
 
    s.tracks().remove(t);
    t.segments().remove(s);
    t.remove(s.links());
    }
 
    //...Refit...
    t.fit();
 
#ifdef TRKRECO_DEBUG
    std::cout << "... refitting" << std::endl;
    t.dump("detail", "2nd> ");
#endif
 
    return bads;
}
 
TSegment *
TConformalFinder::link(const TSegment & base,
               const HepPoint3D & p,
               const HepVector3D & v,
               const AList<TSegment> & candidates,
               AList<TSegment> & alternatives,
               unsigned level) const {
 
#ifdef TRKRECO_DEBUG
    std::cout << "    link : base = " << std::endl;
    base.dump("vector hits mc sort", "->  ");
    std::cout << "    p=" << p << ", v=" << v.unit() << std::endl;
#endif
 
    //...Parameters...
    //static const float maxDistance = 0.02;
    //static const float minDirAngle = 0.97;
 
    //...Candidate loop...
    unsigned n = candidates.length();
    float maxAngle = -999.;
    float minDistance = _linkMaxDistance[level];
    float maxDirAngle = _linkMinDirAngle[level];
 
    TSegment * best = NULL;
    for (unsigned j = 0; j < n; j++) {
    TSegment & current = * candidates[j];
 
#ifdef TRKRECO_DEBUG_DETAIL
    current.dump("vector hits mc sort", "    ");
    std::cout << "        dist,dirAngle,angle=" << current.distance(p, v);
    std::cout << "," << ((current.position() - p).unit()).dot(v.unit());
    std::cout << "," << v.dot(current.direction()) << std::endl;
#endif
    float distance = current.distance(p, v);
    //cout<<"current, Layer: "<<current.links()[0]->hit()->wire()->layerId()
    //  <<" local: "<<current.links()[0]->hit()->wire()->localId()<<endl;
    //cout<<" distance of segs = "<<distance<<" limit = "<<_linkMaxDistance[level]<<endl;
    if (distance > _linkMaxDistance[level]) continue;
 
    HepVector3D dir = (current.position() - p);
    if (dir.x() > M_PI)        dir.setX(dir.x() - PI2);
    else if (dir.x() < - M_PI) dir.setX(PI2 + dir.x());
 
    float dirAngle = dir.unit().dot(v.unit());
    //cout<<" dirAngle of segs = "<<dirAngle<<" limit = "<<_linkMinDirAngle[level]<<endl;
    if (dirAngle < _linkMinDirAngle[level]) continue;   //tmp for tsf
 
//      if (distance < minDistance) {
//          minDistance = distance;
//          best = & current;
//      }
    if (dirAngle > maxDirAngle) {
        maxDirAngle = dirAngle;
        best = & current;
    }
    alternatives.append(current);
    }
 
#ifdef TRKRECO_DEBUG
    std::cout << "    # of best + alternatives = " << alternatives.length() << std::endl;
    if (best) {
    std::cout << "        Best is ";
    best->dump("hits");
    }
#endif
 
    alternatives.remove(best);
    if (best) return best;
    return NULL;
}
 
void
TConformalFinder::deleteTrack(TTrack & t) const {
    const AList<TSegment> & segments = t.segments();
    unsigned n = segments.length();
    for (unsigned i = 0; i < n; i++) {
    TSegment & s = * segments[i];
    s.tracks().remove(t);
    }
 
    delete & t;
}
 
void
TConformalFinder::removeUsedSegments(const AList<TTrack> & tracks) {
    unsigned n = tracks.length();
    for (unsigned i = 0; i < n; i++) {
    TTrack & t = * tracks[i];
    AList<TSegment> & segments = t.segments();
    AList<TSegment> toBeRemoved;
    unsigned nS = segments.length();
    for (unsigned j = 0; j < nS; j++) {
        TSegment & s = * segments[j];
        unsigned sId = s.superLayerId();
//      unsigned as = sId % 2;
//      unsigned id = sId / 2;
        unsigned as;
        unsigned id;
        switch (sId) {
        case 0: as = 1; id = 0;break;
        case 1: as = 1; id = 1;break;
        case 2: as = 0; id = 0;break;
        case 3: as = 0; id = 1;break;
        case 4: as = 0; id = 2;break;
        case 5: as = 1; id = 2;break;
        case 6: as = 1; id = 3;break;
        case 7: as = 1; id = 4;break;
        case 8: as = 1; id = 5;break;
        case 9: as = 0; id = 3;break;
        case 10: as = 0; id = 4;break;
        default: break;
        }
 
        //...Check used links...
        AList<TMLink> links = Links(s, t);
        if (links.length() == 0) {
        s.tracks().remove(t);
        toBeRemoved.append(s);
#ifdef TRKRECO_DEBUG
        std::cout << "!!! why this happends???" << std::endl;
        std::cout << "    no used link" << std::endl;
#endif
        continue;
        }
 
        //...Remove from lists...
        _allUnused[as][id].remove(s);
 
        //...Remove incoming links...
        const AList<TSegment> & outers = s.outerLinks();
        unsigned nO = outers.length();
//          std::cout << " >>> removing outer links" << std::endl;
//          s.dump("hits", "    To ");
        for (unsigned i = 0; i < nO; i++) {
        outers[i]->innerLinks().remove(s);
//          outers[i]->dump("hits", "    From ");
        }
 
        //...Remaining hits...
        AList<TMLink> unused = s.links();
        unused.remove(links);
        s.remove(unused);
        unsigned nUnused = unused.length();
 
        //...Create new segment if too many remaining links...
/*      if (nUnused < _minNLinksForSegment) {
        for (unsigned k = 0; k < nUnused; k++) {
            _unused[as].append(unused[k]);
            _unused[2].append(unused[k]);
        }
        }
        else {
        TSegment * ss = new TSegment(unused);
        AList<TSegment> newSegments = ss->split();   //tmp for tsf
        if (newSegments.length() == 0) {
            ss->solveDualHits();
            newSegments.append(ss);
            for (unsigned k = 0; k < nO; k++)
            outers[k]->innerLinks().append(ss);
        }
        else {
            delete ss;
        }
        _allUnused[as][id].append(newSegments);
        _allSegments[as][id].append(newSegments);
        }
*/
    }
    segments.remove(toBeRemoved);
    }
}
 
void
TConformalFinder::updateTLinks(AList<TTrack> & tracks) {
    unsigned n = tracks.length();
    for (unsigned i = 0; i < n; i++) {
    const AList<TMLink> & links = tracks[i]->links();
    unsigned nL = links.length();
    for (unsigned j = 0; j < nL; j++) {
        TMLink & l = * links[j];
        tracks[i]->approach(l);
    }
    }
}
 
int
TConformalFinder::doit(const AList<TMDCWireHit> & axial,
               const AList<TMDCWireHit> & stereo,
               AList<TTrack> & tracks,
               AList<TTrack> & tracks2D) {
#ifdef TRKRECO_DEBUG
    _stage = ConformalInitialization;
#endif
 
    static bool first = true;
    if (first) {
    first = false;
    int size;
//  _s = (struct summary *)
//      BASF_Sharedmem->get_pointer(BASF_Sharedmem->get_id(),
//                  "TrkConfSum",
//                  & size);
    }
 
    //...For debug...
    if (debugLevel()) {
    std::cout << name() << " ... processing"
          << " axial=" << axial.length()
          << ",stereo=" << stereo.length()
          << ",tracks=" << tracks.length()
          << std::endl;
    }
 
    //...Create TMLinks with conformal position...
//    TConformalFinder0::conformalTransformation(ORIGIN, axial, _allHits[0]);
//    TConformalFinder0::conformalTransformation(ORIGIN, stereo, _allHits[1]);
    // for Tsf
    TConformalFinder0::conformalTransformationDriftCircle(ORIGIN, axial, _allHits[0]);
    TConformalFinder0::conformalTransformationDriftCircle(ORIGIN, stereo, _allHits[1]);
    _allHits[2].append(_allHits[0]);
    _allHits[2].append(_allHits[1]);
 
    //...Select good hits...
    selectGoodHits();
#ifdef TRKRECO_DEBUG
    cout<<"axial Hits:"<<_allHits[0].length()<<"   good axial hits:"<<_hits[0].length()
        <<endl<<"stereo Hits:"<<_allHits[1].length()<<"   good setero hits:"<<_hits[1].length()  
        <<endl;
#endif
        
    //...Segment finding...
    if (_perfectSegmentFinding)
    findSegmentsPerfect();
    else
//  findSegments();
    findSegmentsTsf();
 
#ifdef TRKRECO_DEBUG
    cout<<"axial Seg: "<<_allSegments[0][0].length()<<", "<<_allSegments[0][1].length()
      <<", "<<_allSegments[0][2].length()<<", "<<_allSegments[0][3].length() 
      <<", "<<_allSegments[0][4].length()<<", "<<_allSegments[0][5].length()  
      <<"   stereo Seg: "<<_allSegments[1][0].length()<<", "<<_allSegments[1][1].length() 
      <<", "<<_allSegments[1][2].length()<<", "<<_allSegments[1][3].length()  
      <<", "<<_allSegments[1][4].length()<<", "<<_allSegments[1][5].length()<<endl;
 
    for(int i=0; i<2; ++i) {
      if (i == 0) cout<<"Axial......"<<endl;
      else cout<<"Stereo......"<<endl;
      for(int j=0; j<6; ++j) {
    for(int k=0; k<_allSegments[i][j].length(); ++k){
      cout<<"SL(a/s): "<<j<<" seeds in Seg"<<k<<": "<<_allSegments[i][j][k]->links().length()<<endl;
      for(int kk=0; kk<_allSegments[i][j][k]->links().length(); ++kk)
        cout<<" layerId: "<<_allSegments[i][j][k]->links()[kk]->hit()->wire()->layerId()
          <<"   localId: "<<_allSegments[i][j][k]->links()[kk]->hit()->wire()->localId()<<endl;
    }
      }
    }
#endif
 
    //...Fast finding...
    unsigned level = 0;
    _T0ResetDone = false;
    if (_fastFinder) {
    linkSegments(level);
    fastFinding2D(level);
    updateTLinks(_2DTracks);
 
    //...T0 reset here...
    if (_doT0Reset) {
        std::cout << "TConformalFinder ... T0 reset is done" << std::endl;
        _T0ResetDone = true;
        return 0;
    }
 
#ifdef TRKRECO_WINDOW
    _rphiWindow.skip(false);
#endif
    fastFinding3D(level);
    updateTLinks(_2DTracks);
    updateTLinks(_3DTracks);
 
#ifdef TRKRECO_WINDOW
    _rphiWindow.skip(false);
    displayTracks(_2DTracks, _allUnused, "all 2D after fast level 0");
    displayTracks(_3DTracks, _allUnused, "all 3D after fast level 0");
#endif // TRKRECO_WINDOW
 
 
    //...Fast finding again...
    level = 1;
    linkSegments(level);
    fastFinding2D(level);
    updateTLinks(_2DTracks);
 
#ifdef TRKRECO_WINDOW
    _rphiWindow.skip(false);
#endif
    fastFinding3D(level);
    updateTLinks(_2DTracks);
    updateTLinks(_3DTracks);
 
#ifdef TRKRECO_WINDOW
    _rphiWindow.skip(false);
    displayTracks(_2DTracks, _allUnused, "all 2D after fast level 1");
    displayTracks(_3DTracks, _allUnused, "all 3D after fast level 1");
#endif //TRKRECO_WINDOW
       }
 
    //...Slow finding...
    if (_slowFinder) {
    level = 2;
    linkSegments(level);
    slowFinding2D(level);
    updateTLinks(_2DTracks);
#ifdef TRKRECO_WINDOW
    _rphiWindow.skip(false);
//  _rphiWindow.skipAllWindow(false);
#endif
    fastFinding3D(level);
    updateTLinks(_2DTracks);
    updateTLinks(_3DTracks);
 
#ifdef TRKRECO_WINDOW
    _rphiWindow.skip(false);
//  _rphiWindow.skipAllWindow(false);
    displayTracks(_2DTracks, _allUnused, "all 2D after slow level 2");
    displayTracks(_3DTracks, _allUnused, "all 3D after slow level 2");
#endif //TRKRECO_WINDOW
    }
 
#ifdef TRKRECO_DEBUG
    int zsltrk = _2DTracks.length();
    std::cout << name() << " ... # 3D tracks = " << _3DTracks.length()
             << ", # 2D tracks = " << _2DTracks.length() << std::endl;  
    for (unsigned i = 0; i < zsltrk; i++){
//        _2DTracks[i]->dump(); 
          cout<<"pt:"<<_2DTracks[i]->pt()<<"   impact:"<<_2DTracks[i]->impact()<<endl;
          //          nTrack->fill(0., 1.0);
    }
#endif
 
    //...Mask hits with bad chisq...
    TTrackManager::maskBadHits(_3DTracks, _maxSigma2);
 
    tracks2D = _2DTracks;
    tracks = _3DTracks;
    for(int i=0;i<_3DTracks.length();i++){
      _3DTracks[i]->finder(TrackOldConformalFinder);
      _3DTracks[i]->setFinderType(2);
    }
#ifdef TRKRECO_DEBUG
//    if (_3DTracks.length() == 0) {
//      cout<<"3D failed: 2D length = "<<_2DTracks.length()<<endl;
      for(int nn = 0; nn < tracks2D.length(); ++nn){
    cout<<"2D: "<<nn<<"  radius: "<<tracks2D[nn]->radius()<<endl;
    for (int mm = 0; mm < tracks2D[nn]->links().length(); ++mm)
      cout<<"layer: "<<tracks2D[nn]->links()[mm]->wire()->layerId()
        <<" local: "<<tracks2D[nn]->links()[mm]->wire()->localId()<<endl;
      }
      
      cout<<"unused axial Seg: "<<_allUnused[0][0].length()<<", "<<_allUnused[0][1].length()
    <<", "<<_allUnused[0][2].length()<<", "<<_allUnused[0][3].length()
    <<", "<<_allUnused[0][4].length()<<", "<<_allUnused[0][5].length()
    <<" unused stereo Seg: "<<_allUnused[1][0].length()<<", "<<_allUnused[1][1].length()
    <<", "<<_allUnused[1][2].length()<<", "<<_allUnused[1][3].length()
    <<", "<<_allUnused[1][4].length()<<", "<<_allUnused[1][5].length()
    <<endl;
 
      for(int i=0; i<2; ++i) {
    if (i == 0) cout<<"unused axial hits: "<<endl;
    else cout<<"unused stereo hits: "<<endl;
    for(int k = 0; k < _allHits[i].length(); ++k) {
      if (_allHits[i][k]->hit()->state() & WireHitUsed) continue;
      else cout<<"      layerId: "<<_allHits[i][k]->hit()->wire()->layerId()
        <<"     localId: "<<_allHits[i][k]->hit()->wire()->localId()<<endl;
    }
    if (i == 0) cout<<"unused axial link in segs: "<<endl;
    else cout<<"unused stereo link in segs: "<<endl;
    for(int j=0; j<6; ++j) {
      for(int k=0; k<_allUnused[i][j].length(); ++k){
        cout<<"sl: "<<i<<" "<<j<<"  length of seg "<<k<<": "<<_allUnused[i][j][k]->links().length()<<endl;
        for(int kk=0; kk<_allUnused[i][j][k]->links().length(); ++kk)
          cout<<"   layerId: "<<_allUnused[i][j][k]->links()[kk]->hit()->wire()->layerId()
        <<"     localId: "<<_allUnused[i][j][k]->links()[kk]->hit()->wire()->localId()<<endl;
      }
    }
      }
//    }
#endif
 
    //...Termination...
/*    tracks = _3DTracks;
    tracks2D = _2DTracks;
    ++_s->_nEvents;
    unsigned n3 = _3DTracks.length();
    for (unsigned i = 0; i < n3; i++)
    if (_3DTracks[i]->finder() & TrackFastFinder)
        ++_s->_nTracksFast3D;
    else if (_3DTracks[i]->finder() & TrackSlowFinder)
        ++_s->_nTracksSlow3D;
    unsigned n2 = _2DTracks.length();
    for (unsigned i = 0; i < n2; i++)
    if (_2DTracks[i]->finder() & TrackFastFinder)
        ++_s->_nTracksFast2D;
    else if (_2DTracks[i]->finder() & TrackSlowFinder)
        ++_s->_nTracksSlow2D;
 
    if (debugLevel()) {
    std::cout << name() << " ... # 3D tracks = " << _3DTracks.length()
          << ", # 2D tracks = " << _2DTracks.length() << std::endl;
    }
 
#ifdef TRKRECO_DEBUG
    _stage = ConformalOutside;
#endif
*/
    return 0;
}
 
void
TConformalFinder::fastFinding3D(unsigned level) {
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... fast finding (3D) : level = " << level << std::endl;
#endif
#ifdef TRKRECO_DEBUG
    if (_stage == ConformalSlow2D) {
    _stage = ConformalSlow3D;
    }
    else {
    _stage = ConformalFast3DLevel0;
    if (level == 1)
        _stage = ConformalFast3DLevel1;
    }
#endif
 
    AList<TTrack> tracks3D;
    AList<TTrack> touched;
    AList<TSegment> bads;
    unsigned n = _2DTracks.length();
    for (unsigned i = 0; i < n; i++) {
      
    const TTrack & t = * _2DTracks[i];
 
#ifdef TRKRECO_DEBUG
    cout<<"links in 2Dtracks: "<<t.links().length()<<endl;
    for (int ii = 0; ii < t.links().length(); ++ii) {
        cout<<"layerId: "<<t.links()[ii]->wire()->layerId()
          <<"  localId: "<<t.links()[ii]->wire()->localId()<<endl;
    }
    cout<<"center: "<<t.center()<<" radius = "<<t.radius()<<endl;
#endif
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "==> fast 3D building : " << t.name() << std::endl;
    t.dump("hits sort flag", "    2D track hits=");
#endif
    AList<TSegment> segments = stereoSegments(t);
        //jialk
        unsigned NS = segments.length();
        if ( NS >= 30 ) continue;
 
    //yuany
    //cout<<"stereo segments:"<<segments.length()<<endl;
    //unsigned n = segments.length();
    //for (unsigned i = 0; i < n; i++) {
    //segments[i]->dump("hits");
    //}
 
    AList<TSegment> badSegments;
    if (_stereoMode == 2)
        badSegments = stereoSegmentsFromBadHits(t);
    //yuany
    //cout<<"bad stereo segments:"<<endl;
    //n = badSegments.length();
    //for (unsigned i = 0; i < n; i++) {
    //badSegments[i]->dump("hits");
    //}
 
#ifdef TRKRECO_WINDOW
    displayStatus("Conf::fast3D ... seed");
    _rphiWindow.append(segments, leda_blue);
    _rphiWindow.append(badSegments, leda_red);
    _rphiWindow.oneShot(t, leda_green);
#endif
 
    //cout<<"stereo mode: "<<_stereoMode<<endl;
    
    //...Save a 2D track...
    TTrack * s = NULL;
    if (_stereoMode)
        s = _builder.buildStereoNew(t, segments, badSegments);
    else
        s = _builder.buildStereo(t, segments);
    HepAListDeleteAll(badSegments);
    //cout<<"after 3D build..."<<endl;
    if (! s) {
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "... 3D failure" << std::endl;
#endif
#ifdef TRKRECO_WINDOW
        displayStatus("Conf::fastd3D ... 3D failed");
        _rphiWindow.append(segments, leda_blue);
        _rphiWindow.oneShot(t, leda_red);
#endif      
        continue;
    }
    //...Quality check...
    if (! TTrackManager::goodTrack(* s)) {
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "... 3D failure (bad quality)" << std::endl;
#endif
#ifdef TRKRECO_WINDOW
        displayStatus("Conf::fastd3D ... 3D failed (bad quality)");
        _rphiWindow.append(segments, leda_blue);
        _rphiWindow.oneShot(* s, leda_red);
#endif      
// Remove following sentences because _s is NULL now 
// and it's wrong to use the pointer for some noise events... 05/10/15 zangsl       
//
//      if (s->finder() & TrackFastFinder)
//      ++_s->_nTracksFast2DBadQuality;
//      else if (s->finder() & TrackSlowFinder)
//      ++_s->_nTracksSlow2DBadQuality;
   
        deleteTrack(* s);
        continue;
    }
    
    //...New name...
    s->name(t.name() + "-3D");
 
#ifdef TRKRECO_WINDOW
    displayStatus("Conf::fastd3D ... 3D ok");
    _rphiWindow.append(segments, leda_blue);
    _rphiWindow.oneShot(* s, leda_green);
#endif      
 
    //...Salvage by segments...
    salvage(* s, 3, bads);
    //cout<<"3D finder done!"<<endl;
    tracks3D.append(s);
    touched.append(_2DTracks[i]);
 
    s->assign(WireHitConformalFinder);
    s->quality(0);
 
    //...Segment...
    static AList<TTrack> tmp;
    tmp.removeAll();
    tmp.append(s);
    removeUsedSegments(tmp);
        // liucy :add constraint of number of tracks
        if(tracks3D.length()>20)break;
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... 3D finished : " << s->name() << std::endl;
    s->dump("detail", "    ");
#endif
#ifdef TRKRECO_WINDOW
    displayStatus("Conf::fastd3D ... finished");
    _rphiWindow.oneShot(* s, leda_green);
#endif      
    }
 
    // resolveSegments(tracks3D);
    _3DTracks.append(tracks3D);
    _2DTracks.remove(tracks3D);
    _2DTracks.remove(touched);
    for (unsigned i = 0; i < touched.length(); i++) {
    deleteTrack(* touched[i]);
    // saved[i]->fit();
    }
}
 
#ifdef TRKRECO_WINDOW
void
TConformalFinder::displayStatus(const std::string & m) const {
    _rphiWindow.clear();
    _rphiWindow.text(m);
    _rphiWindow.append(_allHits[2], leda_pink);
    _rphiWindow.append(_hits[2], leda_cyan);
    _rphiWindow.append(_unused[2]);
    displayAppendSegments(_allSegments, leda_grey1);
    displayAppendSegments(_allUnused, leda_orange);
}
 
void
TConformalFinder::displayAppendSegments(const AList<TSegment> a[2][6],
                    leda_color c) const {
    for (unsigned i = 0; i < 2; i++) {
    for (unsigned j = 0; j < 6; j++) {
        const AList<TSegment> & segments = a[i][j];
        unsigned n = segments.length();
        for (unsigned k = 0; k < n; k++) {
        _rphiWindow.append(* segments[k], c);
        }
    }
    }
}
 
void
TConformalFinder::displayTracks(const AList<TTrack> & l,
                const AList<TSegment> segments[2][6],
                const std::string & c) const {
    displayStatus("Conf::display ... " + c);
    for (unsigned i = 0; i < l.length(); i++)
    _rphiWindow.append(* l[i], leda_green);
    _rphiWindow.wait();
}
#endif
 
bool
TConformalFinder::quality2D(TTrack & t, unsigned level) const {
#ifdef TRKRECO_WINDOW
    std::string s = "Conf::quality2D ... level " + itostring(level);
#endif
    //jialk caution origin is 15.
    //...Check # radius...Liuqg
    if(fabs(t.radius()) < 15.) return false;  //corresponding to pt = 300*1T*0.50/2(m).
 
    //...Check # of segments(superlayers)...
    // unsigned nSuperLayers = NSuperLayers(t.links(), 2);
    unsigned nSuperLayers = NSuperLayers(t.links());
    if (nSuperLayers < _minNSegments) {
#ifdef TRKRECO_DEBUG
    std::cout << "... quality check : level=" << level << " : bad" << std::endl;
    std::cout << "    reason : # segments(superlayer) =" << nSuperLayers;
    std::cout << " < " << _minNSegments << std::endl;
#endif
#ifdef TRKRECO_WINDOW
    s += " rejected because of few segments(superlayers)";
    displayStatus(s);
    _rphiWindow.oneShot(t, leda_red);
#endif
    return false;
    }
    if (level == 0) {
#ifdef TRKRECO_WINDOW
    s += " ok : # of used segments(superlayers) = ";
    s += itostring(nSuperLayers);
    s += " > " + itostring(_minNSegments);
    displayStatus(s);
    _rphiWindow.oneShot(t, leda_green);
#endif
 
    return true;
    }
 
    //...Check superlayer usage...
    unsigned n3 = NSuperLayers(t.links(), 3);
 
    //...Check super layer usage...
    if (n3 != nSuperLayers) {
    unsigned sl = SuperLayer(t.links(), 2);
    // unsigned sl = SuperLayer(t.links(), 3);
    if (sl == 0) {
#ifdef TRKRECO_DEBUG
        std::cout << "... quality check : level = " << level << " : bad";
        std::cout << std::endl;
        std::cout << "    reason : super layer pattern(n2) = 0 " << std::endl;
#endif
#ifdef TRKRECO_WINDOW
        s += " rejected because of bad super-layer pattern(n2=0)";
        displayStatus(s);
        _rphiWindow.oneShot(t, leda_red);
#endif
        return false;
    }
    //unsigned fl = 0;   //belle, liuqg
    //while ((sl & (1 << (fl * 2))) == 0) ++fl;
    //bool empty = false;
 
    /*unsigned fl[5] = {2, 3, 4, 9, 10};     //tmp....liuqg
    unsigned nfl = 0;
    for (unsigned i = 0; i < 5; i++) {
        if ((sl & (1 << fl[i])) == 1) {
        nfl = i;
        break;
        }
    }
        bool empty = false;
    for (unsigned i = nfl + 1; i < 5; i++) {
        bool thisLayer = (sl & (1 << fl[i]));
        if (thisLayer && empty) {
 
//#ifdef TRKRECO_DEBUG
        std::cout << "... quality check : level = " << level;
        std::cout << " : bad" << std::endl;
        std::cout << "    reason : super layer pattern = ";
//      for (unsigned j = 0; j < 6; j++) std::cout << (sl >> (j * 2)) % 2;
        std::cout << std::endl;
//#endif
#ifdef TRKRECO_WINDOW
        s += " rejected because of bad super-layer pattern";
        displayStatus(s);
        _rphiWindow.oneShot(t, leda_red);
#endif
        return false;
        }
        if (! thisLayer) empty = true;
    }*/
    }
 
#ifdef TRKRECO_DEBUG
    std::cout << "... quality check : level = " << level;
    std::cout << " : ok, super layer pattern = " << std::endl;
    std::string ptn;
    unsigned sl = SuperLayer(t.links(), 2);
    for (unsigned j = 0; j < 6; j++) {
    unsigned k = (sl >> (j * 2)) % 2;
    std::cout << k;
    if (k) ptn += "1";
    else   ptn += "0";
    }
    std::cout << std::endl;
#endif
#ifdef TRKRECO_WINDOW
    s += " ok : super layer ptn = " + ptn;
    displayStatus(s);
    _rphiWindow.oneShot(t, leda_green);
#endif
 
    return true;
}
 
void
TConformalFinder::fastFinding2D(unsigned level) {
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... fast finding (2D) : level = " << level << std::endl;
#endif
#ifdef TRKRECO_DEBUG
    _stage = ConformalFast2DLevel0;
    if (level == 1)
    _stage = ConformalFast2DLevel1;
#endif
 
    AList<TSegment> * segments = & _allUnused[0][0];
 
    unsigned idBase = _2DTracks.length() + _3DTracks.length();
    unsigned nTrial = 0;
    unsigned seedLayer = 5;
    //yuany
/*    cout<<"before loop "<<endl;
    for (unsigned i = 0; i < 5; i++) {
      for(unsigned j=0;j<_allUnused[0][i].length();j++){
    (_allUnused[0][i])[j]->dump("hits");
    (_allUnused[0][i])[j]->dump("link");
      }
    }*/
    
//   while ((seedLayer--) > 2) {  2006/04/27
    while ((seedLayer--) > 1) {
      
      //...Seed loop...
      AList<TTrack> trackCandidates;
      unsigned n = segments[seedLayer].length();
 
      for (unsigned i = 0; i < n; i++) {
    TSegment & seed = * segments[seedLayer][i];
    std::string name = "f" + itostring(nTrial + idBase);
    //yuany
    //seed.dump("hits");
    //seed.dump("link");
#ifdef TRKRECO_DEBUG
    std::cout << ">>> fast 2D : super layer " << seedLayer << ", ";
    std::cout << i << " : ";
    seed.dump("link");
#endif
    
    //...Check uniqueness...
    AList<TSegment> seeds;
    //cout<<"NUiqueLinks:"<<NUniqueLinks(seed)<<"; NMajorLinks:"<<NMajorLinks(seed)<<endl;
    //cout<<"innerLinks size:"<<seed.innerLinks().length()<<endl;
    if (NUniqueLinks(seed) > 0) {
      seeds = UniqueLinks(seed);
    }
    else if (NMajorLinks(seed) > 0) {
      seeds = MajorLinks(seed);
    }
    else {
      continue;
    }
    seeds.append(seed);
 
#ifdef TRKRECO_DEBUG
    cout<<"seeds in fastFinding2D: "<<seeds.length()<<endl;
    for (int kk = 0; kk < seeds.length(); ++kk)
      cout<<"LayerId: "<<seeds[kk]->links()[0]->wire()->layerId()
        <<" LocalId: "<<seeds[kk]->links()[0]->wire()->localId()<<endl;
#endif
 
        //...Refine...
    refine:
 
        //solve LR of seeds;
        for (int k = 0; k < seeds.length(); ++k)   //tsf
          seeds[k]->solveLR();
 
#ifdef TRKRECO_WINDOW
        displayStatus("Conf::fast2D ... seed segments");
        _rphiWindow.oneShot(seeds, leda_green);
#endif
 
        //...Try to build a track...
#ifdef TRKRECO_DEBUG
        std::cout << "==> fast building 2D : seed layer = " << seedLayer;
        std::cout << ", " << name << std::endl;
#endif
 
        TTrack * t = _builder.buildRphi(seeds);
        ++nTrial;
        //cout<<"buildRphi done."<<endl;
 
        //...Track check...
        if (! t) continue;
        t->name(name);
        bool ok = quality2D(* t, 0);
        //cout<<"track quality " <<ok<<endl;
        if (! ok) {
        deleteTrack(* t);
        continue;
        }
        //cout<<"length of seg in 2d build = "<<t->segments().length()<<"  cores in 2d = "<<t->nCores()<<endl;
 
        //...Bad segment rejection...
        AList<TSegment> bads = removeBadSegments(* t);
        //cout<<"test point no0;  bad segment length: "<<bads.length()<<endl;
 
        // This was taken out because the goto loop continues forever for some noise events. 2005/10/14
        if (bads.length()) {
          ok = quality2D(* t, 1);
          if (! ok) {
        seeds = refineSegments(* t);
        if (seeds.length() >= _minNSegments) {
          deleteTrack(* t);
          goto refine;
        }
        deleteTrack(* t);
        continue;
          }
        }
        //...Salvage by segments...
        salvage(* t, 1, bads);
        //cout <<"test point no1"<<endl;
        //refineLinks(* t, 3);   //change the MinN from 3 to 2, zangshilei 20060307
            refineLinks(* t, 2);
        //cout<<"test point no1.5"<<endl;
 
        //change the level from 2 to 0, 20060428; because the program stops(MCv5.1.0/run10237.dat)
        //release by Liuqg
        ok = quality2D(* t, 2); 
        if (! ok) {
        deleteTrack(* t);
        continue;
        }
        //cout<<"test point no2"<<endl;
        
        //...Append segments...
        if (NSuperLayers(t->links()) < 5) {
        t = expand(* t);
        }
        //cout<<"test point no3"<<endl;
 
        //t->dump("detail", "  ");
#ifdef TRKRECO_DEBUG
        std::cout << "... 2D finished : " << t->name() << std::endl;
        t->dump("hits sort flag pull", "    ");
#endif
#ifdef TRKRECO_WINDOW
        displayStatus("Conf::fast2D ... finished");
        _rphiWindow.oneShot(* t, leda_green);
#endif
        t->finder(TrackFastFinder);
        t->quality(TrackQuality2D);
        t->fitting(TrackFitGlobal);
        trackCandidates.append(t);
    }
 
    //...Resolve multi-track candidates...
    resolveSegments(trackCandidates);
 
    //...Remove used segments...
    //yuany
/*  cout<<"before remove "<<endl;
      for (unsigned i = 0; i < 5; i++) {
    for(unsigned j=0;j<_allUnused[0][i].length();j++){
      (_allUnused[0][i])[j]->dump("hits");
      (_allUnused[0][i])[j]->dump("link");
    }
      }
*/
    removeUsedSegments(trackCandidates);
    //yuany
//  cout<<"after remove "<<endl;
/*      for (unsigned i = 0; i < 5; i++) {
    for(unsigned j=0;j<_allUnused[0][i].length();j++){
      (_allUnused[0][i])[j]->dump("hits");
      (_allUnused[0][i])[j]->dump("link");
    }
      }
*/
    _2DTracks.append(trackCandidates);
    }
 
    resolveHits(_2DTracks);
}
 
const TMDCWire *
TConformalFinder::conformal2Wire(const HepPoint3D & p) {
    static const TMDC & cdc = * TMDC::getTMDC();
 
    std::cout << "p = " << p << std::endl;
    float r = sqrt(4. / p.mag2());
    float phi = p.phi();
    return cdc.wire(r, phi);
}
 
AList<TSegment>
TConformalFinder::pickUpSegmentsInConformal(float phi[12],
                        float loadWidth,
                        unsigned axialStereoSwitch) const {
    AList<TSegment> outList;
    HepPoint3D center(1., 1., 0.);
    center.setPhi(phi[0]);
 
    //...Search for segments...
    for (unsigned sl = 0; sl < 2; sl++) {
    if (sl == 0) {
        if (! (axialStereoSwitch & 1)) continue;
    }
    else {
        if (! (axialStereoSwitch & 2)) continue;
    }
    for (unsigned i = 0; i < 6; i++) {
        const AList<TSegment> & list = _allUnused[sl][i];
        unsigned n = list.length();
        for (unsigned j = 0; j < n; j++) {
        TSegment & s = * list[j];
 
        const HepPoint3D & p = s.position();
        if (center.dot(p) < 0.) continue;
 
        bool found = false;
        const AList<TMLink> & inners = s.inners();
        unsigned m = inners.length();
        for (unsigned k = 0; k < m; k++) {
            float dPhi = phi[i * 2 + sl] -
            inners[k]->position().phi();
            dPhi = fabs(fmod(dPhi, PI2));
            if (dPhi > (PI2 - dPhi)) dPhi = PI2 - dPhi;
            if (dPhi < WIDTH[i * 2 + sl] * loadWidth) {
            outList.append(s);
            found = true;
            break;
            }
        }
 
        if (found) continue;
        const AList<TMLink> & outers = s.outers();
        m = outers.length();
        for (unsigned k = 0; k < m; k++) {
            float dPhi = phi[i * 2 + sl + 1] -
            outers[k]->position().phi();
            dPhi = fabs(fmod(dPhi, PI2));
            if (dPhi > (PI2 - dPhi)) dPhi = PI2 - dPhi;
            if (dPhi < WIDTH[i * 2 + sl] * loadWidth) {
            outList.append(s);
            break;
            }
        }
        }
    }
    }
 
    return outList;
}
 
AList<TMLink>
TConformalFinder::pickUpLinksInConformal(float phi[12], 
                     float loadWidth,
                     unsigned axialStereoSwitch) const {
    HepPoint3D center(1., 1., 0.);
    center.setPhi(phi[0]);
 
    AList<TMLink> outList;
    unsigned nBad = _unused[2].length();
    for (unsigned i = 0; i < nBad; i++) {
    unsigned sl = _unused[2][i]->wire()->superLayerId();
    unsigned as = sl % 2;
    if (as == 0) {
        if (! (axialStereoSwitch & 1)) continue;
    }
    else {
        if (! (axialStereoSwitch & 2)) continue;
    }   
 
    const HepPoint3D & p = _unused[2][i]->position();
    if (center.dot(p) < 0.) continue;
 
    float dPhi = phi[sl] - _unused[2][i]->position().phi();
    dPhi = fabs(fmod(dPhi, PI2));
    if (dPhi > (PI2 - dPhi)) dPhi = PI2 - dPhi;
    if (dPhi < WIDTH[sl] * loadWidth) {
        outList.append(_unused[2][i]);
        continue;
    }
    dPhi = phi[sl + 1] - _unused[2][i]->position().phi();
    dPhi = fabs(fmod(dPhi, PI2));
    if (dPhi > (PI2 - dPhi)) dPhi = PI2 - dPhi;
    if (dPhi < WIDTH[sl] * loadWidth)
        outList.append(_unused[2][i]);
    }
    return outList;
}
 
int
TConformalFinder::crossPointsInConformal(const AList<TSegment> & inList,
                     HepPoint3D points[12]) const {
    //...Parameters...
    static const float confRadius2[] = {4. / ( 8.3 *  8.3),
                    4. / (16.9 * 16.9),
                    4. / (21.7 * 21.7),
                    4. / (31.3 * 31.3),
                    4. / (36.1 * 36.1),
                    4. / (44.1 * 44.1),
                    4. / (50.5 * 50.5),
                    4. / (58.5 * 58.5),
                    4. / (64.9 * 64.9),
                    4. / (72.9 * 72.9),
                    4. / (79.3 * 79.3),
                    4. / (87.4 * 87.4)};
 
    //...Get conformal points from segments as seeds...
    AList<TMLink> forLine;
    HepPoint3D center;
    unsigned n = inList.length();
    for (unsigned i = 0; i < n; i++) {
    const HepPoint3D & p = inList[i]->position();
    forLine.append(new TMLink(NULL, NULL, p));
    center += p;
    }
    center *= 1. / float(n);
 
    //...Make a line in the conformal plane...
    TMLine line(forLine);
    int err = line.fit();
    if (err) {
#ifdef TRKRECO_DEBUG
    std::cout << "crossPoints ... failed due to line fit failure" << std::endl;
#endif
    HepAListDeleteAll(forLine);
    return -1;
    }
 
    //...Calculate points...
    for (unsigned i = 0; i < 12; i++) {
    HepPoint3D p[2];
    float c = - line.a() * line.b();
    float d = 1. + line.a() * line.a();
    float e = sqrt(confRadius2[i] * d - line.b() * line.b());
    p[0].setX((c + e) / d);
    p[0].setY(line.a() * p[0].x() + line.b());
    p[1].setX((c - e) / d);
    p[1].setY(line.a() * p[1].x() + line.b());
 
    float mag0 = (center - p[0]).mag();
    float mag1 = (center - p[1]).mag();
 
    if (mag0 < mag1) points[i] = p[0];
    else             points[i] = p[1];
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "0,1  = " << p[0] << ", " << p[1] << std::endl;
#endif
    }
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... center is : " << center << std::endl;
    std::cout << "... cross points are : " << std::endl;
    for (unsigned i = 0; i < 12; i++) {
    std::cout << "    " << i << " : " << points[i] << std::endl;
    }
#endif
#ifdef TRKRECO_WINDOW
    displayStatus("Conf::crossPoints ... line to salvage for conf plane");
    _rphiWindow.append(inList, leda_green);
    _rphiWindow.oneShot(line, leda_blue);
#endif
 
    HepAListDeleteAll(forLine);
    return 0;
}
 
AList<TSegment>
TConformalFinder::stereoSegments(const TTrack & t) const {
#ifdef TRKRECO_DEBUG
    std::cout << "... finding stereo segments" << std::endl;
#endif
 
    const AList<TSegment> & bases = (const AList<TSegment> &) t.segments();
    AList<TSegment> seeds;
    unsigned n = bases.length();
    if (n == 0) return seeds;
 
    //...Calculate cross points in the conformal plane...
    TPoint2D points[12];
    int err = crossPoints(t, points);
    //yuany
//   for(unsigned i=0;i<12;i++)
//     cout<<points[i]<<endl;
    if (err) {
#ifdef TRKRECO_DEBUG
    std::cout << "... stereo segment finding failed" << std::endl;
#endif
    return seeds;
    }
 
    //...Pick up segments...
    // return pickUpSegments(points, float(_stereoLoadWidth), 2);
    AList<TSegment> list = pickUpSegments(points, float(_stereoLoadWidth), 2);
    //yuany
    //cout<<"list length "<<list.length()<<endl;
 
    //...Save direction of a segment of axial layers...
    TPoint2D dir[5];
    for (unsigned i = 0; i < n; i++) {
    const TSegment & s = * bases[i];
    AList<TMLink> ll = s.links();
    unsigned sl = ll[0]->wire()->axialStereoLayerId()/4;
    dir[sl] = TPoint2D(s.direction());
//  cout<<"in sl = "<<sl<<"   dir = "<<dir[sl]<<endl;
    }
 
    //...Cal. direction if empty...
    for (unsigned i = 0; i < 5; i++) {
    if (dir[i].mag() < .5) {
        unsigned j = i;
        while ((j < 5) && (dir[j].mag() < .5))
        ++j;
        if (j > 4) j = 4;
        if (dir[j].mag() < .5) {
        j = i;
        while ((j > 0) && (dir[j].mag() < .5))
            --j;
        }
        dir[i] = dir[j];
    }
    }
    //...Remove bad segments...
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    ... direction :" << std::endl;
    for (unsigned i = 0; i < 6; i++)
    std::cout << "        " << i << " : " << dir[i] << std::endl;
    std::cout << "    ... direction cos :" << std::endl;
#endif
    AList<TSegment> badList;
    AList<TSegment> toBeChecked[6];
    unsigned nL = list.length();
    for (unsigned i = 0; i < 6; ++i) 
    toBeChecked[i].removeAll();
 
    for (unsigned j = 0; j < nL; ++j) {
    TSegment & s = * list[j];
    AList<TMLink> ll = s.links();
    unsigned sl = ll[0]->wire()->axialStereoLayerId()/4;
    toBeChecked[sl].append(s);
    }
 
    for (unsigned i = 0; i < 6; ++i) {
    unsigned nToBeChecked = toBeChecked[i].length();
    if (nToBeChecked < 4) continue;
    unsigned nBadList = badList.length();
    for (unsigned j = 0; j < nToBeChecked; ++j) {
        if((badList.length() - nBadList) >= nToBeChecked - 4) break;
        TSegment & s = * toBeChecked[i][j];
        TPoint2D sDir = s.direction();
        unsigned sl = i;
        unsigned aSl;
 
        if (sl < 2) aSl = 0;
        else if (sl < 4) aSl = 2;
        else aSl = 3;
        if (dir[aSl].dot(sDir) < 0.7) badList.append(s);
        if (sl == 2 || sl == 3 || sl == 4) {
        if (dir[aSl - 1].dot(sDir) < 0.6) badList.append(s);
        }
        else {
        if (dir[aSl + 1].dot(sDir) < 0.6) badList.append(s);
        }
    }
    toBeChecked[i].remove(badList);
    if (toBeChecked[i].length() > 4) {
        AList<TSegment> tmp;  //with >= 4 links
        unsigned nSegs = toBeChecked[i].length();
        for (unsigned j = 0; j < nSegs; ++j) {
        TSegment & s = * toBeChecked[i][j];
        if (s.links().length() > 3) tmp.append(s);
        }
        toBeChecked[i].remove(tmp);
        for (unsigned j = 0; j < tmp.length(); ++j) {
        AList<TMLink> tmpLinks = tmp[j]->links();
        for (unsigned k = 0; k < toBeChecked[i].length(); ++k) {
            TSegment & ss = * toBeChecked[i][k];
            AList<TMLink> threeLinks = ss.links();
            for (unsigned kk = 0; kk < threeLinks.length(); ++kk) {
            TMLink & ll = * threeLinks[kk];
            if (tmpLinks.hasMember(ll)) {
                badList.append(ss);
                break;
            }
            }
        }
        }//j
    }
//jialk
       toBeChecked[i].remove(badList);
       if (toBeChecked[i].length() > 4) {//still with > 4 segs after above filter
             AList<TSegment> temp;
             unsigned nSegs = toBeChecked[i].length();
             for(unsigned k=0; k < nSegs; ++k){
                TSegment &s = *toBeChecked[i][k];
                temp.append(s);
             }
             //use bubble to get a sequence
             for(unsigned y=0; y<nSegs; y++){
                  for(unsigned x=1; x<nSegs-y; x++){
                      TSegment &s1 = *temp[x];
                      TSegment &s2 = *temp[x-1];
                      if(s1.links().length() < s2.links().length()){
                           TSegment &  s0 = *temp[x-1];
                           *temp[x-1] = *temp[x];
                           *temp[x] = s0;
                      }
                  }  
             }
           for(unsigned j=4; j<nSegs; ++j){
                TSegment &s4 = *temp[j];
                badList.append(s4);
           }
         }
 
           
    
    toBeChecked[i].remove(badList);
 
//    cout<<"**********the num of stereo segs after filter***********  "<<toBeChecked[i].length()<<endl;
    for(unsigned j = 0; j < toBeChecked[i].length(); ++j){
        TSegment &s5 = *toBeChecked[i][j];
//        cout<<"num of hits"<<s5.links().length()<<endl;
    }
 
   }
 
/*    for (unsigned i = 0; i < nL; i++) {
    TSegment & s = * list[i];
    AList<TMLink> ll = s.links();
    unsigned sl = ll[0]->wire()->axialStereoLayerId()/4;
    TPoint2D sDir = s.direction();
    unsigned aSl;
 
    if (sl < 2) aSl = 0;
    else if (sl < 4) aSl = 2;
    else aSl = 3;
    if (dir[aSl].dot(sDir) < 0.7) badList.append(s);
    if (sl == 2 || sl == 3 || sl == 4) {
        if (dir[aSl - 1].dot(sDir) < 0.6) badList.append(s);
    }
    else {
        if (dir[aSl + 1].dot(sDir) < 0.6) badList.append(s);
    }
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "        " << dir[sl].dot(sDir) << ", ";
    std::cout << dir[sl + 1].dot(sDir) << " : ";
    s.dump("hits sort");
#endif
    }
*/
    //...Width cut...
//      for (unsigned i = 0; i < nL; i++) {
//      TSegment & s = * list[i];
//      unsigned width = s.width();
//      if (width > 2) badList.append(s);
//      }
 
//    cout<<"stereoSegment: "<<list.length()<<"  badList: "<<badList.length()<<endl;
    list.remove(badList); 
//    cout<<"after remove bads: "<<list.length()<<endl;
 
#ifdef TRKRECO_DEBUG
    std::cout << "    ... bad segments :" << std::endl;
    for (unsigned i = 0; i < badList.length(); i++)
    badList[i]->dump("hits sort", "        ");
#endif
 
    return list;
}
 
HepVector3D
TConformalFinder::direction(const TSegment & segment) const {
    AList<TMLink> forLine;
    const TSegment * s = & segment;
    while (1) {
    if (s->outerLinks().length() != 1) break;
    const TSegment * o = s->outerLinks()[0];
    const HepPoint3D & p = o->position();
    forLine.append(new TMLink(NULL, NULL, p));  
    s = o;
    }
 
    if (forLine.length() == 0)
    return segment.direction();
    else if (forLine.length() < 2) {
    HepAListDeleteAll(forLine);
    return segment.direction();
    }
 
    //...Make a line in the conformal plane...
    TMLine line(forLine);
    int err = line.fit();
    if (err) {
#ifdef TRKRECO_DEBUG
    std::cout << "direction ... failed due to line fit failure" << std::endl;
#endif
    HepAListDeleteAll(forLine);
    return segment.direction();
    }
 
    HepAListDeleteAll(forLine);
    Vector3 tmp(-1., -(line.a() + line.b()));
    tmp.unit();
    return HepVector3D(tmp);
}
 
 
void
TConformalFinder::salvage(TTrack & track,
              unsigned asSwitch,
              const AList<TSegment> & bads) const {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... salvaging in real plane" << std::endl;
#endif
 
    //...Calculate cross points...
    TPoint2D points[12];
    int err = crossPoints(track, points);
    if (err) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... salvaging failed in calculation of intersection"
          << std::endl;
#endif
    return;
    }
 
    //...Pick up segments...
    AList<TSegment> toBeChecked = pickUpSegments(points,
                         float(_salvageLoadWidth),
                         asSwitch);
    toBeChecked.remove(bads);
    toBeChecked.remove(track.segments());
    toBeChecked = trackSide(track, toBeChecked);
 
    //...Pick up links...
    AList<TMLink> links;
    unsigned nB = toBeChecked.length();
    for (unsigned i = 0; i < nB; i++) {
    const AList<TMLink> & tmp = toBeChecked[i]->links();
    unsigned nL = tmp.length();
    for (unsigned j = 0; j < nL; j++) {
        if (tmp[j]->hit()->track()) continue;
        links.append(tmp[j]);
    }
    }
 
    //...Search in bad hits...
    AList<TMLink> all = pickUpLinks(points, float(_salvageLoadWidth), asSwitch);
#ifdef TRKRECO_DEBUG
    std::cout<<" pickUpLinks length = "<<all.length()<<std::endl;
#endif
 
    all = trackSide(track, all);
    all.remove(track.links());
#ifdef TRKRECO_DEBUG
    std::cout<<" after reomove pickUpLinks length = "<<all.length()<<std::endl;
#endif
    links.append(all);
 
#ifdef TRKRECO_WINDOW
    AList<TMLink> tmp = links;
#endif
 
    //...Ask builder...
    _builder.salvage(track, links);
 
    //...Check used segments...
    const AList<TMLink> & usedLinks = track.links();
    for (unsigned i = 0; i < nB; i++) {
    TSegment & segment = * toBeChecked[i];
    AList<TMLink> used = Links(segment, track);
    if (used.length() == 0) continue;
 
    track.segments().append(segment);
    segment.tracks().append(track);
    }
 
#ifdef TRKRECO_WINDOW
    displayStatus("Conf::salvage ... results");
    TTrackBase y(tmp);
    _rphiWindow.append(y, leda_red);
    _rphiWindow.append(toBeChecked, leda_red);
    _rphiWindow.oneShot(track, leda_green);
#endif
 
    //...Termination...
    return;
}
 
int
TConformalFinder::crossPoints(const TTrack & track,
                  TPoint2D points[12]) const {
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    cross points are :" << std::endl;
#endif
 
    /*   
    //...Parameters...
    static const float R[] = { 8.3,
                  16.9,
                  21.7,
                  31.3,
                  36.1,
                  44.1,
                  50.5,
                  58.5,
                  64.9,
                  72.9,
                  79.3,
                  87.4};
    static const float R2[] = {  8.3 *  8.3,
                16.9 * 16.9,
                21.7 * 21.7,
                31.3 * 31.3,
                36.1 * 36.1,
                44.1 * 44.1,
                50.5 * 50.5,
                58.5 * 58.5,
                64.9 * 64.9,
                72.9 * 72.9,
                79.3 * 79.3,
                87.4 * 87.4};
    */
    static const float R[] = { 7.3,
                  12.21,
                  18.9,
                  25.375,
                  31.86,
                  39.16,
                  45.685,
                  52.215,
                  58.665,
                  65.91,
                  72.37,
                  80.};
    static const float R2[] = { R[0]*R[0],
                                R[1]*R[1],
                                R[2]*R[2],
                                R[3]*R[3],
                                R[4]*R[4],
                                R[5]*R[5],
                                R[6]*R[6],
                                R[7]*R[7],
                                R[8]*R[8],
                                R[9]*R[9],
                                R[10]*R[10],
                                R[11]*R[11]};
 
    
    static const TPoint2D o(0., 0.);
    TPoint2D c = track.center();
    TPoint2D co = - c;
    double r = fabs(track.radius());
    double r2 = r * r;
    double d2 = c.mag2();
    double d = sqrt(d2);
    double sl = - c.x() / c.y();
    //...Calculate points...
    unsigned nOk = 0;
    for (unsigned i = 0; i < 12; i++) {
    double minR = r < R[i] ? r : R[i];
    double maxR = r < R[i] ? R[i] : r;
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    r,R ... " << minR << ", " << maxR << ", " << d << std::endl;
#endif
 
    if ((r + R[i] < d) || (minR + d < maxR)) { //no intersection of two circles
        points[i] = o;
        continue;
    }
    ++nOk;
    double a = R2[i] + d2 - r2;
    double s = sqrt(4. * R2[i] * d2 - a * a);
    double q = 0.5 * a / c.y();
    points[i].x(0.5 * (c.x() * a + c.y() * s) / d2);
    points[i].y(q + sl * points[i].x());
    const double Bz = -1000*(m_pmgnIMF->getReferField());
    if (co.cross(points[i] - c) * track.charge() * Bz > 0.) {
        points[i].x(0.5 * (c.x() * a - c.y() * s) / d2);
        points[i].y(q + sl * points[i].x());
    }
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    " << i << " : " << points[i] << std::endl;
    std::cout << "    chg=" << track.charge()<<std::endl;
    std::cout << ", c=" << c << ", co=" <<std::endl;
    std::cout << ", " << co.cross(points[i] - c) * track.charge() << std::endl;
    std::cout << "        " << 0.5 * (c.x() * a + c.y() * s) / d2;
    std::cout << ", " << q + sl * (0.5 * (c.x() * a + c.y() * s) / d2);
    std::cout << "        " << 0.5 * (c.x() * a - c.y() * s) / d2;
    std::cout << ", " << q + sl * (0.5 * (c.x() * a - c.y() * s) / d2);
    std::cout << std::endl;
#endif
 
    }
    if (nOk) return 0;
    else return -1;
}
 
AList<TSegment>
TConformalFinder::pickUpSegments(const TPoint2D x[12],
                 float loadWidth,
                 unsigned axialStereoSwitch) const {
    static const TPoint2D O(0., 0.);
    AList<TSegment> outList;
 
    //yuany
    //...Search for segments...
    for (unsigned sl = 0; sl < 2; sl++) {
    unsigned nMax = 5;
    if (sl == 0) {
        if (! (axialStereoSwitch & 1)) continue;
    }
    else {
        nMax = 6;
        if (! (axialStereoSwitch & 2)) continue;
    }
    for (unsigned i = 0; i < nMax; i++) {
 
      unsigned layerId;
      switch(sl){
      case 0:
        switch(i){
        case 0:layerId = 2;break;
        case 1:layerId = 3;break;
        case 2:layerId = 4;break;
        case 3:layerId = 9;break;
        case 4:layerId = 10;break;
        }
      case 1:
        switch(i){
        case 0:layerId = 0;break;
        case 1:layerId = 1;break;
        case 2:layerId = 5;break;
        case 3:layerId = 6;break;
        case 4:layerId = 7;break;
        case 5:layerId = 8;break;
        }
      }
     //       cout<<"sl "<<sl<<" i "<<i<<" layerId "<<layerId<<endl;
 
        if (x[layerId] == O) continue;
        float a = WIDTH[layerId] * loadWidth;
        float phi0 = x[layerId].phi();
        float phi1 = x[layerId+1].phi();//Liuqg
        float phiC = (phi0 + phi1) / 2.;
        float phiD = fabs(phi0 - phiC);
        if (phiD > M_PI_2) {
        phiC += M_PI;
        phiD = M_PI - phiD;
        }
        phiD += a;
 
        const AList<TSegment> & list = _allUnused[sl][i];
        unsigned n = list.length();
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "    pickUpSegments ... phi0,1,C,D=" << phi0 << ",";
        std::cout << phi1 << "," << phiC << "," << phiD <<" nhits "<<n<< std::endl;
#endif
        for (unsigned j = 0; j < n; j++) {
        TSegment & s = * list[j];
 
#ifdef TRKRECO_DEBUG_DETAIL
        s.dump("hits sort", "        ");
#endif
 
        bool found = false;
        const AList<TMLink> & inners = s.inners();
        unsigned m = inners.length();
        for (unsigned k = 0; k < m; k++) {
            float phi = inners[k]->position().phi();
            if (phi < 0.) phi += PI2;
            float dPhi = fabs(phi - phiC);
            if (dPhi > M_PI) dPhi = PI2 - dPhi;
#ifdef TRKRECO_DEBUG_DETAIL
            std::cout << "            " << inners[k]->wire()->name();
            std::cout << ", phi,dPhi=" << phi << "," << dPhi << std::endl;
#endif
            if (dPhi < phiD) {
            outList.append(s);
            found = true;
            break;
            }
        }
 
        if (found) continue;
        const AList<TMLink> & outers = s.outers();
        m = outers.length();
        for (unsigned k = 0; k < m; k++) {
            float phi = outers[k]->position().phi();
            if (phi < 0.) phi += PI2;
            float dPhi = fabs(phi - phiC);
            if (dPhi > M_PI) dPhi = PI2 - dPhi;
#ifdef TRKRECO_DEBUG_DETAIL
            std::cout << "            " << outers[k]->wire()->name();
            std::cout << ", phi,dPhi=" << phi << "," << dPhi << std::endl;
#endif
            if (dPhi < phiD) {
            outList.append(s);
            found = true;
            break;
            }
        }
        }
    }
    }
    return outList;
}
 
AList<TMLink>
TConformalFinder::pickUpLinks(const TPoint2D x[12],
                  float loadWidth,
                  unsigned axialStereoSwitch) const {
 
    static const TPoint2D O(0., 0.);
    AList<TMLink> outList;
 
    unsigned nBad = _unused[2].length();
    for (unsigned i = 0; i < nBad; i++) {
    unsigned sl = _unused[2][i]->wire()->superLayerId();
    unsigned as = 1;
    if (_unused[2][i]->wire()->axial()) as = 0;
    if (as == 0) {
        if (! (axialStereoSwitch & 1)) continue;
    }
    else {
        if (! (axialStereoSwitch & 2)) continue;
    }
 
    float a = WIDTH[sl] * loadWidth;
    float phi0 = x[sl].phi();
    float phi1 = x[sl + 1].phi();
    float phi = _unused[2][i]->position().phi();
 
    if (phi < 0.) phi += PI2;
    if (phi1 < phi0) {
        phi1 = phi0;
        phi0 = x[sl + 1].phi();
    }
    float dPhi = phi1 - phi0;
    if (dPhi < M_PI) {
        phi0 -= a;
        phi1 += a;
        if (phi > phi0 && phi < phi1) outList.append(_unused[2][i]);
    }
    else {
        phi0 += a;
        phi1 -= a;
        if (phi < phi0 || phi > phi1) outList.append(_unused[2][i]);
    }
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "TConformalFinder:: pickUpLinks " << std::endl;
    std::cout << "a "<<a <<" phi0 "<<phi0<<" phi1 "<<phi1<<" phi "<< phi<<" dPhi "<<dPhi << std::endl;
#endif
    }
    return outList;
}
 
void
TConformalFinder::slowFinding2D(unsigned level) {
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... slow finding (2D) : level = " << level << std::endl;
#endif
#ifdef TRKRECO_DEBUG
    _stage = ConformalSlow2D;
#endif
 
    AList<TSegment> * segments = & _allUnused[0][0];
 
    unsigned id = _2DTracks.length() + _3DTracks.length();
    unsigned seedLayer = 5;
    while (seedLayer--) {
 
    //...Seed loop...
    AList<TSegment> tmp = segments[seedLayer];
    unsigned n = tmp.length();
    for (unsigned i = 0; i < n; i++) {
        AList<TTrack> trackCandidates;
        TSegment & current = * tmp[i];
        if (current.links().length() < 3) continue;
        if (current.innerLinks().length() != 1) continue;
        if (current.innerLinks()[0]->links().length() < 3) continue;
        AList<TSegment> seeds;
        seeds.append(current);
        seeds.append(current.innerLinks()[0]);
 
        std::string name = "s" + itostring(id);
 
#ifdef TRKRECO_DEBUG
        std::cout << "==> slow building : " << name << std::endl;
#endif
 
        //...Try to build a track...
        for (int k = 0; k < seeds.length(); ++k)   //tsf
        seeds[k]->solveLR();
 
        TTrack * t = expand(seeds);
        if (t) {
        AList<TSegment> bads;
        t->name(name);
        salvage(* t, 1, bads);
        if (! trackQuality(* t)) {
            deleteTrack(* t);
            continue;
        }
        t->finder(TrackSlowFinder);
        t->quality(TrackQuality2D);
        t->fitting(TrackFitGlobal);
        trackCandidates.append(t);
 
#ifdef TRKRECO_DEBUG
        std::cout << "... 2D finished : " << t->name() << std::endl;
        t->dump("hits sort flag pull", "    ");
#endif
#ifdef TRKRECO_WINDOW
        displayStatus("Conf::expand ... finished");
        _rphiWindow.oneShot(* t, leda_green);
#endif      
        removeUsedSegments(trackCandidates);
        _2DTracks.append(trackCandidates);
        id = _2DTracks.length() + _3DTracks.length();
        }
    }
 
    //...Resolve multi-track candidates...
    // resolveSegments(trackCandidates);
 
    //...Remove used segments...
    }
}
 
TTrack *
TConformalFinder::expand(AList<TSegment> & seeds) const {
#ifdef TRKRECO_WINDOW
    displayStatus("Conf::expand ... seeds");
    _rphiWindow.oneShot(seeds, leda_green);
#endif      
#ifdef TRKRECO_DEBUG
    std::cout << "... expand : # of seeds = " << seeds.length() << std::endl;
#endif
    TTrack * t = NULL;
    if (seeds.length() > 2) {
    TTrack * t = _builder.buildRphi(seeds);
    if (t) {
        if (! trackQuality(* t)) {
        deleteTrack(* t);
        t = NULL;
        }
    }
    }
    if (t == NULL) {
    AList<TMLink> links = Links(seeds);
    TCircle c(links);
    c.fit();
    t = new TTrack(c);
    t->fit();
    t->segments().append(seeds);
    }
    return expand(* t);
}
 
TTrack *
TConformalFinder::expand(TTrack & ti) const {
#ifdef TRKRECO_WINDOW
    displayStatus("Conf::expand ... seed track");
    _rphiWindow.oneShot(ti, leda_green);
#endif      
#ifdef TRKRECO_DEBUG
    std::cout << "... expand : by a track" << std::endl;
#endif
 
    TTrack * t = & ti;
    AList<TSegment> seeds = t->segments();
    unsigned nSegments = seeds.length();
    unsigned nCores = t->nCores();
    unsigned nTrial = 0;
    while (1) {
 
    //...Target superlayer...
    unsigned sl = SuperLayer(t->cores());
    unsigned inner;
    unsigned outer;
    targetSuperLayer(sl, inner, outer);
//cout<<"sl: "<<sl<<"  inner: "<<inner<<"  outer: "<<outer<<endl;
    unsigned target = inner;
    if (target == 11)
        target = outer;
    if (target == 11)
        break;
 
    //...Get target segments...
    AList<TSegment> segments = targetSegments(* t, target);
    if (segments.length() == 0) {
        unsigned fl[5] = {2, 3, 4, 9, 10};
        unsigned inner2 = 0;
        unsigned outer2 = 5;
        for (int i = 0; i < 5; i++) {
        if (fl[i] == inner) inner2 = i;
        if (fl[i] == outer) outer2 = i;
        }
        if (inner > 2&& inner2>0) target = fl[inner2 - 1];
        if (target == 9) {
        target = fl[outer2 + 1];
        }
        if (target == 13) break;
        segments = targetSegments(* t, target);
    }
#ifdef TRKRECO_DEBUG
    std::cout << "... expand : segments to be checked = ";
    std::cout << segments.length() << std::endl;
#endif
#ifdef TRKRECO_WINDOW
    displayStatus("Conf::expand ... candidate segments");
    _rphiWindow.append(segments, leda_blue);
    _rphiWindow.oneShot(* t, leda_green);
#endif      
 
    //...Create candidate tracks...
    unsigned n = segments.length();
    AList<TTrack> tracks;
    for (unsigned i = 0; i < n; i++) {
        segments[i]->solveLR();
        AList<TSegment> seed0 = seeds;
        seed0.append(segments[i]);
        TTrack * t0 = _builder.buildRphi(seed0);
        if (! t0) continue;
        if ((t0->segments().length() <= nSegments) ||
        (t0->nCores() <= nCores)               ||
        (SuperLayer(t0->cores()) <= sl)) {
        deleteTrack(* t0);
        continue;
#ifdef TRKRECO_DEBUG
        std::cout << "... expand : candidate deleted" << std::endl;
#endif
        }
        tracks.append(t0);
    }
 
#ifdef TRKRECO_DEBUG
    std::cout << "... expand : target superlayer = " << target;
    std::cout << ", # of track candidates = " << tracks.length()
          << std::endl;
#endif
 
    //...Expand by links...
    if (tracks.length() == 0) {
        AList<TMLink> links = targetLinks(* t, target);
        _builder.salvage(* t, links);
        unsigned newSl = SuperLayer(t->cores());
        if (newSl & (1 << (target))) {
#ifdef TRKRECO_DEBUG
        std::cout << "... expand : links to be appended = ";
        std::cout << links.length() << std::endl;
#endif
#ifdef TRKRECO_WINDOW
        displayStatus("Conf::expand ... candidate links");
        _rphiWindow.append(links, leda_blue);
        _rphiWindow.oneShot(* t, leda_green);
#endif
        TTrack * t0 = new TTrack(* t);
        if (! t0) break;
        if ((t0->nCores() <= nCores) ||
            (SuperLayer(t0->cores()) <= sl)) {
            deleteTrack(* t0);
            break;
#ifdef TRKRECO_DEBUG
            std::cout << "... expand : candidate deleted" << std::endl;
#endif
        }
        tracks.append(t0);
        }
        else
        break;
    }
 
#ifdef TRKRECO_DEBUG
    std::cout << "... expand : # of track candidates = " << tracks.length()
          << std::endl;
#endif
 
    //...Select best...
    unsigned nTracks = tracks.length();
    if (nTracks == 0) break;
    TTrack * tNew = tracks[0];
    unsigned nBestCores = tNew->nCores();
    for (unsigned i = 1; i < nTracks; i++) {
        if (tracks[i]->nCores() > nBestCores) {
        deleteTrack(* tNew);
        tNew = tracks[i];
        nBestCores = tNew->nCores();
        }
        else {
        deleteTrack(* tracks[i]);
        }
    }
    nSegments = tNew->segments().length();
    nCores = nBestCores;
    seeds = tNew->segments();
 
    //...Quality check...
    if (! trackQuality(* tNew)) {
        deleteTrack(* tNew);
        break;
    }
    deleteTrack(* t);
    t = tNew;
    }
 
#ifdef TRKRECO_DEBUG
    if (t == & ti) std::cout << "... expand : failed" << std::endl;
    else std::cout << "...  expand : track expanded" << std::endl;
#endif
#ifdef TRKRECO_WINDOW
    displayStatus("Conf::expand ... finished");
    _rphiWindow.oneShot(* t, leda_green);
#endif
    
    return t;
}
 
void
TConformalFinder::targetSuperLayer(unsigned sl,
                   unsigned & inner,
                   unsigned & outer) const {
    inner = 11;
    outer = 11;
    bool innerFound = false;
    bool outerFound = false;
    for (unsigned i = 0; i < 5; i++) {             //changed to BESIII,Liuqg
    unsigned fl[5] ={2, 3, 4, 9, 10};
    if (! innerFound) {
        if (sl & (1 << (fl[i]))) innerFound = true;
        else inner = fl[i];
    }
    if (! outerFound) {
        if (sl & (1 << (fl[4 - i]))) outerFound = true;
        else outer = fl[4 - i];
    }
    }
}
 
AList<TSegment>
TConformalFinder::targetSegments(const TTrack & t, unsigned sl) const {
    AList<TSegment> targets;
 
    //...Calculate cross points...
    TPoint2D points[12];
    int err = crossPoints(t, points);
    if (err) {
#ifdef TRKRECO_DEBUG
    std::cout << "    ... no target found" << std::endl;
#endif
    return targets;
    }
 
    //...Pick up segments...
    AList<TSegment> toBeChecked = pickUpSegments(points, 3, 1);
    unsigned n = toBeChecked.length();
    for (unsigned i = 0; i < n; i++) {
    if (toBeChecked[i]->superLayerId() == sl)
        targets.append(toBeChecked[i]);
    }
 
    return targets;
}
 
AList<TMLink>
TConformalFinder::targetLinks(const TTrack & t, unsigned sl) const {
    AList<TMLink> targets;
 
    //...Calculate cross points...
    TPoint2D points[12];
    int err = crossPoints(t, points);
    if (err) {
#ifdef TRKRECO_DEBUG
    std::cout << "    ... no target found" << std::endl;
#endif
    return targets;
    }
 
    //...Pick up segments...
    AList<TMLink> toBeChecked = pickUpLinks(points, 3, 1);
    unsigned n = toBeChecked.length();
    for (unsigned i = 0; i < n; i++) {
    if (toBeChecked[i]->wire()->superLayerId() == sl)
        targets.append(toBeChecked[i]);
    }
 
    return targets;
}
 
AList<TSegment>
TConformalFinder::refineSegments(const TTrack & t) const {
    const AList<TSegment> & original = t.segments();
    AList<TSegment> outList;
    unsigned n = original.length();
    for (unsigned i = 0; i < n; i++) {
    AList<TMLink> tmp = Links(* original[i], t);
    if (tmp.length() >= _minNLinksForSegmentInRefine)
        outList.append(original[i]);
    }
    unsigned sl = SuperLayer(outList);
    unsigned nCMax = 0;
    unsigned nStart = 0;
    unsigned nC = 0;
    unsigned nS = 0;
    unsigned fl[5] = {2, 3, 4, 9, 10};
    for (unsigned i = 0; i < 5; i++) {
        if (sl & (1 << fl[i])) {
        ++nC;
        if (nC == 1) nS = i;
        if (nC > nCMax) {
        nCMax = nC;
        nStart = nS;
        }
    }
    else {
        nC = 0;
    }
    }
//    nStart *= 2;
//    nCMax *= 2;
 
    outList.removeAll();
    if (nCMax >= _minNSegments) {
    for (unsigned i = 0; i < n; i++) {
        unsigned id = original[i]->superLayerId();
        if ((id >= fl[nStart]) && (id < fl[nStart + nCMax]))
        outList.append(original[i]);
    }
    }
 
#ifdef TRKRECO_DEBUG
    std::cout << "... refine segments : orignal sl = ";
    bitDisplay(sl, 11, 0);
    std::cout << ", output sl = ";
    bitDisplay(SuperLayer(outList), 11, 0);
    std::cout << std::endl;
#endif
 
    return outList;
}
 
bool
TConformalFinder::trackQuality(const TTrack & t) const {
    unsigned n1 = NSuperLayers(t.links());
    unsigned n3 = NSuperLayers(t.links(), 2);
 
    //...check # radius...Liuqg
    if(fabs(t.radius())<15.) return false;
 
#ifdef TRKRECO_WINDOW
    std::cout << "... trackQuality : n1,n3,nMissing=" << n1 << "," << n3;
    std::cout << "," << NMissingAxialSuperLayers(t.links()) << std::endl;
#endif
#ifdef TRKRECO_WINDOW
    displayStatus("trackQuality");
    if ((n3 < 2) || (NMissingAxialSuperLayers(t.links()) > 1))
    _rphiWindow.oneShot(t, leda_red);
    else
    _rphiWindow.oneShot(t, leda_green);
#endif
    if (n3 < 2) return false;
//    if (NMissingAxialSuperLayers(t.links()) > 1) return false;
 
    return true;
}
 
void
TConformalFinder::refineLinks(TTrack & t, unsigned minN) const {
    const AList<TMLink> & links = t.links();
    AList<TMLink> sl[11];
    unsigned n = links.length();
    for (unsigned i = 0; i < n; i++)
    sl[links[i]->wire()->superLayerId()].append(links[i]);
    AList<TMLink> toBeRemoved;
    for (unsigned i = 0; i < 11; i++)
    if (sl[i].length() < minN)
        toBeRemoved.append(sl[i]);
    t.remove(toBeRemoved);
#ifdef TRKRECO_DEBUG
    std::cout << "... refining links : removed links = " << toBeRemoved.length();
    std::cout << std::endl;
#endif
}
 
AList<TMLink>
TConformalFinder::trackSide(const TTrack & t, const AList<TMLink> & a) const {
    static const TPoint2D o(0., 0.);
    TPoint2D c = t.center();
    TPoint2D co = - c;
 
    AList<TMLink> tmp;
    unsigned n = a.length();
    for (unsigned i = 0; i < n; i++) {
    TPoint2D x = a[i]->wire()->xyPosition();
    const double Bz = -1000*(m_pmgnIMF->getReferField());
    if (co.cross(x - c) * t.charge() * Bz < 0.)
        tmp.append(a[i]);
    }
    return tmp;
}
 
AList<TSegment>
TConformalFinder::trackSide(const TTrack & t,
                const AList<TSegment> & a) const {
    static const TPoint2D o(0., 0.);
    TPoint2D c = t.center();
    TPoint2D co = - c;
 
    AList<TSegment> tmp;
    unsigned n = a.length();
    for (unsigned i = 0; i < n; i++) {
    const AList<TMLink> & b = a[i]->links();
    bool bad = false;
    unsigned nB = b.length();
    for (unsigned j = 0; j < nB; j++) {
        TPoint2D x = b[j]->wire()->xyPosition();
        //yzhang 2012-05-03 
        const double Bz = -1000*(m_pmgnIMF->getReferField());
        if (co.cross(x - c) * t.charge()* Bz > 0.) {
        bad = true;
        break;
        }
    }
    if (bad) continue;
    tmp.append(a[i]);
    }
    return tmp;
}
 
void
TConformalFinder::resolveHits(AList<TTrack> & tracks) const {
    unsigned nTracks = tracks.length();
    if (nTracks < 2) return;
 
    for (unsigned i = 0; i < nTracks - 1; i++) {
    TTrack & t0 = * tracks[i];
    const AList<TMLink> & links0 = t0.links();
    unsigned n0 = links0.length();
    AList<TMLink> remove0;
    for (unsigned j = i + 1; j < nTracks; j++) {
        TTrack & t1 = * tracks[j];
        const AList<TMLink> & links1 = t1.links();
        unsigned n1 = links1.length();
        AList<TMLink> remove1;
 
        //...Check overlap...
        for (unsigned k = 0; k < n0; k++) {
        TMLink & l = * links0[k];
 
        //...Decide which is better...
        if (links1.hasMember(l)) {
            TMLink l0(NULL, l.hit());
            TMLink l1(NULL, l.hit());
            t0.approach(l0);
            t1.approach(l1);
            if (l0.pull() > l1.pull())
            remove0.append(l);
            else
            remove1.append(l);
        }
        }
 
        if (remove1.length()) {
        t1.remove(remove1);
        t1.fit();
        }
    }
 
    if (remove0.length()) {
        t0.remove(remove0);
        t0.fit();
    }
    }
}
 
AList<TSegment>
TConformalFinder::stereoSegmentsFromBadHits(const TTrack & t) const {
 
    AList<TSegment> output;
    for (unsigned i = 0; i < 6; i++)
    output.append(new TSegment());
 
    //...Calculate cross points...
    TPoint2D points[12];
    int err = crossPoints(t, points);
    if (err) {
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... conf::stereoSegBadHits : "
          << "error in calculation of intersection" << std::endl;
#endif
    return output;
    }
 
    static const TPoint2D O(0., 0.);
    unsigned nBad = _unused[2].length();
    for (unsigned i = 0; i < nBad; i++) {
    unsigned sl = _unused[2][i]->wire()->superLayerId();
    unsigned asl = _unused[2][i]->wire()->axialStereoLayerId()/4;
    unsigned as;
    if (_unused[2][i]->wire()->axial())
        as = 0;
    else as = 1;
//  unsigned as = sl % 2;
    if (as == 0) continue;
 
    float a = WIDTH[sl] * _salvageLoadWidth;
    float phi0 = points[sl].phi();
    float phi1 = points[sl + 1].phi();
    float phi = _unused[2][i]->position().phi();
    if (phi < 0.) phi += PI2;
    if (phi1 < phi0) {
        phi1 = phi0;
        phi0 = points[sl + 1].phi();
    }
    float dPhi = phi1 - phi0;
    if (dPhi < M_PI) {
        phi0 -= a;
        phi1 += a;
        if (phi > phi0 && phi < phi1)
//      output[sl / 2]->append(* _unused[2][i]);
        output[asl / 4]->append(* _unused[2][i]);
    }
    else {
        phi0 += a;
        phi1 -= a;
        if (phi < phi0 || phi > phi1)
        output[asl / 4]->append(* _unused[2][i]);
    }
    }
 
    return output;
}
 
void
TConformalFinder::findSegmentsPerfect(void) {
 
    //...Create lists of links for each super layer...
    AList<TMLink> links[11];
    unsigned nHits = _hits[2].length();
    for (unsigned i = 0; i < nHits; i++) {
    TMLink & l = * _hits[2][i];
    links[l.wire()->superLayerId()].append(l);
    }
 
    //...MC tracks...
    const unsigned nHep = TTrackHEP::list().length();
 
    //...Loop over each super layer...
    for (unsigned i = 0; i < 11; i++) {
    unsigned superLayerId = i / 2;
    unsigned id = i % 2;
 
    //...Counters...
    AList<TMLink> ** seeds =
        (AList<TMLink> **) malloc(nHep * sizeof(AList<TMLink> *));
    for (unsigned j = 0; j < nHep; j++)
        seeds[j] = NULL;
 
    //...Hit loop...
    unsigned n = links[i].length();
    for (unsigned j = 0; j < n; j++) {
        TMLink & l = * links[i][j];
        const TMDCWireHitMC * mc = l.hit()->mc();
        if (! l.hit()->mc()) {
        std::cout << "TConformalFinder::findSegmentsPerfect !!! "
              << "no MC info. found ... aborted" << std::endl;
        return;
        }
        const unsigned id = mc->hep()->id();
 
        if (! seeds[id])
        seeds[id] = new AList<TMLink>();
        seeds[id]->append(l);
    }
 
    //...Create segments...
    AList<TSegment> segments;
    for (unsigned j = 0; j < nHep; j++) {
        if (! seeds[j]) continue;
 
        const unsigned length = seeds[j]->length();
        if (length < _minNLinksForSegment) continue;
        if (length > _maxNLinksForSegment) continue;
        if (Width(* seeds[j]) > _maxWidthForSegment) continue;
 
        TSegment & s = * new TSegment();
        s.append(* seeds[j]);
        segments.append(s);
    }
    _allSegments[id][superLayerId].append(segments);
    _allUnused[id][superLayerId] = _allSegments[id][superLayerId];
 
    //...Clear...
    for (unsigned j = 0; j < nHep; j++) {
        if (seeds[j])
        delete seeds[j];
    }
    free(seeds);
    }
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "... segment finding(perfect) finished" << std::endl;
#endif
}
 
void
TConformalFinder::findSegmentsTsf(void) {
    //...Create lists of links for each super layer...
    AList<TMLink> links[11];
    unsigned n = _hits[2].length();
    for (unsigned i = 0; i < n; i++) {
    TMLink & l = * _hits[2][i];
    links[l.wire()->superLayerId()].append(l);
    }
 
    //find Segments by tsf
    TMDCTsf * tsf[11];
 
    //initial
    for (unsigned i = 0; i < 11; ++i)
    tsf[i] = new TMDCTsf(i);
 
    for (unsigned i = 0; i < 11; i++) {
    unsigned superLayerId;
    unsigned id;
        switch(i){
        case 0: superLayerId = 0; id = 1; break;
        case 1: superLayerId = 1; id = 1; break;
        case 2: superLayerId = 0; id = 0; break;
        case 3: superLayerId = 1; id = 0; break;
        case 4: superLayerId = 2; id = 0; break;
        case 5: superLayerId = 2; id = 1; break;
        case 6: superLayerId = 3; id = 1; break;
        case 7: superLayerId = 4; id = 1; break;
        case 8: superLayerId = 5; id = 1; break;
        case 9: superLayerId = 3; id = 0; break;
        case 10: superLayerId = 4; id = 0; break;
        default: break;         
        }
    //find segments
    AList<TSegment> tmp = tsf[i]->segments(links[i]);
 
    //...Store them...
    _allSegments[id][superLayerId].append(tmp);
    _allUnused[id][superLayerId] = _allSegments[id][superLayerId];
    delete tsf[i];
 
    //...remove used in seg...
//  for (int j = 0; j < _allUnused[id][superLayerId].length(); ++j){
//      AList<TMLink> used = _allUnused[id][superLayerId][j]->links();
//      _unused[id].remove(used);
//      _unused[2].remove(used);
//  }
    }
}