TConformalFinder.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TConformalFinder.cxx,v 1.49 2022/01/28 22:04:42 maqm Exp $
//-----------------------------------------------------------------------------
// Filename : TConformalFinder.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 <algorithm>
#include <string.h>
// #include "CLHEP/String/Strings.h"
// #include "basf/basfshm.h"
#include "TrkReco/TCircle.h"
#include "TrkReco/TConformalFinder.h"
#include "TrkReco/TConformalFinder0.h"
#include "TrkReco/THistogram.h"
#include "TrkReco/TMDC.h"
#include "TrkReco/TMDCTsf.h"
#include "TrkReco/TMDCUtil.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TMLine.h"
#include "TrkReco/TTrack.h"
#include "TrkReco/TTrackHEP.h"
#include "TrkReco/TTrackManager.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 )
    , m_pmgnIMF( nullptr )
#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));
}
 
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[j];
      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];
      ////std::cout <<"multiLinks "<<nL<<" "<<j<<" "<<&l<< std::endl;
      // if( nullptr == &l ) {
      //    std::cout <<"yzhang debug 720 skip TMLink == 0 multiLinks "<< std::endl;
      //   continue;
      // }
 
      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 {
  std::cout << std::fixed << std::setprecision( 6 ); // yzhang debug for 720
 
#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 );
  std::cout << std::defaultfloat; // yzhang debug for 720
  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 = 1; // add initialize 25-05-15
      unsigned id = 0; // add initialize 25-05-15
      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() > 1 )
  {
    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 " + std::to_string( 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 += std::to_string( nSuperLayers );
    s += " > " + std::to_string( _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" + std::to_string( 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 * ( getPmgnIMF()->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" + std::to_string( 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 * ( getPmgnIMF()->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 * ( getPmgnIMF()->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);
    //  }
  }
}
 
IBesMagFieldSvc* TConformalFinder::getPmgnIMF() const {
  if ( !m_pmgnIMF )
  {
    StatusCode scmgn = Gaudi::svcLocator()->service( "MagneticFieldSvc", m_pmgnIMF );
    if ( scmgn != StatusCode::SUCCESS )
    { std::cout << "Unable to open Magnetic field service" << std::endl; }
  }
  return m_pmgnIMF;
}