TTrackManager.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TTrackManager.cxx,v 1.70 2022/01/28 22:04:42 maqm Exp $
//-----------------------------------------------------------------------------
// Filename : TTrackManager.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A manager of TTrack information to make outputs as Reccdc_trk.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
/* for isnan */
#if defined( __sparc )
#  if defined( __EXTENSIONS__ )
#    include <cmath>
#  else
#    define __EXTENSIONS__
#    include <cmath>
#    undef __EXTENSIONS__
#  endif
#elif defined( __GNUC__ )
#  if defined( _XOPEN_SOURCE )
#    include <cmath>
#  else
#    define _XOPEN_SOURCE
#    include <cmath>
#    undef _XOPEN_SOURCE
#  endif
#endif
 
#define TTRACKMANAGER_INLINE_DEFINE_HERE
#include <cstdlib>
#include <cstring>
#include <values.h>
 
// #include "belle.h"
//  #include "CLHEP/String/Strings.h"
// #include "basf/basfshm.h"
#include "MdcTables/MdcTables.h"
#include "MdcTables/TrkTables.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TTrack.h"
#include "TrkReco/TTrackHEP.h"
#include "TrkReco/TTrackMC.h"
#include "TrkReco/TTrackManager.h"
// #include "TrkReco/TrkReco.h"
 
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IDataManagerSvc.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/StatusCode.h"
 
#include "EventModel/EventModel.h"
#include "GaudiKernel/ContainedObject.h"
#include "GaudiKernel/ObjectVector.h"
#include "GaudiKernel/SmartRef.h"
#include "GaudiKernel/SmartRefVector.h"
 
#include "EvTimeEvent/RecEsTime.h"
 
#include "Identifier/Identifier.h"
#include "Identifier/MdcID.h"
 
// #include "TrkReco/TSvdAssociator.h" // jtanaka 000925
 
// extern BasfSharedMem * BASF_Sharedmem;
 
TTrackManager::TTrackManager()
    : _maxMomentum( 10. )
    , _sigmaCurlerMergeTest( sqrt( 100. ) )
    , _nCurlerMergeTest( 4 )
    , _debugLevel( 0 )
    , _fitter( "TTrackManager Fitter" )
    , _cFitter( "TTrackManager 2D Fitter" )
    , _s( 0 )
    , m_pmgnIMF( nullptr )
    , m_rawDataProviderSvc( nullptr ) {}
 
TTrackManager::~TTrackManager() {}
 
std::string TTrackManager::version( void ) const { return std::string( "2.27" ); }
 
void TTrackManager::dump( const std::string& msg, const std::string& pre ) const {
  bool def = ( msg == "" ) ? true : false;
  /*
      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, "TrkMgrSum", & size);
              s._nEvents += r._nEvents;
              s._nTracks += r._nTracks;
              s._nTracksAll += r._nTracksAll;
              s._nTracks2D += r._nTracks2D;
              s._nTracksFinal += r._nTracksFinal;
              s._nSuperMoms += r._nSuperMoms;
              s._nToBeMerged += r._nToBeMerged;
              s._nToBeMergedMoreThanTwo += r._nToBeMergedMoreThanTwo;
          }
 
          std::cout << pre;
          std::cout << "all events : " << s._nEvents << std::endl;
          std::cout << pre;
          std::cout << "all tracks : " << s._nTracksAll << std::endl;
          std::cout << pre;
          std::cout << "    good tracks  : " << s._nTracks << std::endl;
          std::cout << pre;
          std::cout << "    2D tracks    : " << s._nTracks2D << std::endl;
          std::cout << pre;
          std::cout << "    final tracks : " << s._nTracksFinal << std::endl;
          std::cout << pre;
          std::cout << "    super mom.   : " << s._nSuperMoms << std::endl;
          std::cout << pre;
          std::cout << "    to be mreged : " << s._nToBeMerged << std::endl;
          std::cout << pre;
          std::cout << "    to be mreged2: " << s._nToBeMergedMoreThanTwo
                    << std::endl;
      }
  */
  if ( def || msg.find( "eventSummary" ) != std::string::npos ||
       msg.find( "detail" ) != std::string::npos )
  {
    std::cout << pre << "tracks reconstructed : " << _tracksAll.length() << std::endl;
    std::cout << pre << "good tracks          : " << _tracks.length() << std::endl;
    std::cout << pre << "2D tracks            : " << _tracks2D.length() << std::endl;
    std::cout << pre << "Track list:" << std::endl;
 
    std::string tab = pre;
    std::string spc = " ";
    for ( unsigned i = 0; i < _tracksAll.length(); i++ )
    {
      std::cout << tab << TrackDump( *_tracksAll[i] ) << std::endl;
      if ( msg.find( "hits" ) != std::string::npos ||
           msg.find( "detail" ) != std::string::npos )
        Dump( _tracksAll[i]->links(), "hits sort flag" );
    }
  }
}
 
void TTrackManager::maskCurlHits( const AList<TMDCWireHit>& axial,
                                  const AList<TMDCWireHit>& stereo,
                                  const AList<TTrack>&      tracks ) const {
  //...Coded by jtanaka...
 
  int i = 0;
  while ( const TTrack* t = tracks[i++] )
  {
    int j = 0;
    while ( TMDCWireHit* a = axial[j++] )
    {
      double x = t->helix().center().x() - a->xyPosition().x();
      double y = t->helix().center().y() - a->xyPosition().y();
      double r = sqrt( x * x + y * y );
      double R = fabs( t->helix().radius() );
      double q = t->helix().center().x() * a->xyPosition().y() -
                 t->helix().center().y() * a->xyPosition().x();
      double qq = q * t->charge();
      if ( R - 2. < r && r < R + 2. && qq > 0. ) { a->state( a->state() | WireHitUsed ); }
    }
    j = 0;
    while ( TMDCWireHit* s = stereo[j++] )
    {
      double x = t->helix().center().x() - s->xyPosition().x();
      double y = t->helix().center().y() - s->xyPosition().y();
      double r = sqrt( x * x + y * y );
      double R = fabs( t->helix().radius() );
      double q = t->helix().center().x() * s->xyPosition().y() -
                 t->helix().center().y() * s->xyPosition().x();
      double qq = q * t->charge();
      if ( R - 2.5 < r && r < R + 2.5 && qq > 0. ) { s->state( s->state() | WireHitUsed ); }
    }
  }
}
 
void TTrackManager::salvage( const AList<TMDCWireHit>& hits ) const {
 
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << name() << " ... salvaging" << std::endl;
  std::cout << "    # of given hits=" << hits.length() << std::endl;
#endif
 
  //...Check arguments...
  unsigned nTracks = _tracks.length();
  if ( nTracks == 0 ) return;
  unsigned nHits = hits.length();
  if ( nHits == 0 ) return;
 
  //...Hit loop...
  for ( unsigned i = 0; i < nHits; i++ )
  {
    TMDCWireHit& h = *hits[i];
 
    //...Already used?...
    if ( h.state() & WireHitUsed ) continue;
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    checking " << h.wire()->name() << std::endl;
#endif
 
    //...Select the closest track to a hit...
    TTrack* best = closest( _tracks, h );
#ifdef TRKRECO_DEBUG_DETAIL
    if ( !best )
    {
      std::cout << "        no track candidate returned";
      std::cout << "by TTrackManager::closest" << std::endl;
    }
#endif
    if ( !best ) continue;
 
    //...Try to append this hit...
    AList<TMLink> link;
    link.append( new TMLink( 0, &h ) );
    best->appendByApproach( link, 30. );
    // best->assign(WireHitConformalFinder);
    best->finder( TrackTrackManager );
  }
}
 
TTrack* TTrackManager::closest( const AList<TTrack>& tracks, const TMDCWireHit& hit ) const {
 
  TMLink t;
  t.hit( &hit );
  unsigned n           = tracks.length();
  double   minDistance = MAXDOUBLE;
  TTrack*  minTrk      = NULL;
 
  //...Loop over all tracks...
  for ( unsigned i = 0; i < n; i++ )
  {
    TTrack& trk = *tracks[i];
    int     err = trk.approach( t );
    if ( err < 0 ) continue;
    if ( minDistance > t.distance() )
    {
      minDistance = t.distance();
      minTrk      = &trk;
    }
  }
 
  return minTrk;
}
 
StatusCode
// TTrackManager::makeTds(bool doClear, int tkStat) { //yzhang change interface
TTrackManager::makeTds( RecMdcTrackCol* trkcol, RecMdcHitCol* hitcol, int tkStat, int brunge,
                        int cal ) { // yzhang change interface 2010-05-14
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
 
  MsgStream log( msgSvc, "TrkReco" );
 
  IDataProviderSvc* eventSvc = NULL;
  Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
 
#ifdef TRKRECO_DEBUG
  if ( eventSvc ) { log << MSG::INFO << "makeTds:event Svc has been found" << endmsg; }
  else
  {
    log << MSG::FATAL << "makeTds:Could not find eventSvc" << endmsg;
    return StatusCode::FAILURE;
  }
#endif
  // check whether  Recon  already registered
  //    DataObject *aReconEvent;
  //    eventSvc->findObject("/Event/Recon",aReconEvent);
  //    if(!aReconEvent) {
  //       ReconEvent* recevt = new ReconEvent;
  //       StatusCode sc = eventSvc->registerObject("/Event/Recon",recevt );
  //       if(sc!=StatusCode::SUCCESS) {
  //        log << MSG::FATAL << "Could not register ReconEvent" <<endmsg;
  //        return( StatusCode::FAILURE);
  //      }
  //   }
  //
  //   /// Unregister Tracks
  //   IDataManagerSvc *dataManSvc;
  //   if(doClear){//yzhang add, do NOT clear Tds when associate rec
  //     dataManSvc= dynamic_cast<IDataManagerSvc*>(eventSvc);
  //     DataObject *aTrackCol;
  //     eventSvc->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
  //     if(aTrackCol != NULL) {
  //       dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
  //       eventSvc->unregisterObject("/Event/Recon/RecMdcTrackCol");
  //     }
  //     DataObject *aRecHitCol;
  //     eventSvc->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
  //     if(aRecHitCol != NULL) {
  //       dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
  //       eventSvc->unregisterObject("/Event/Recon/RecMdcHitCol");
  //     }
  //   }
  //   /// Writing
  //   RecMdcTrackCol* trkcol = new RecMdcTrackCol;
  //   RecMdcHitCol* hitcol = new RecMdcHitCol;
 
  unsigned n      = _tracks.length();
  int      nTdsTk = trkcol->size();
  for ( int i = 0; i < n; i++ )
  {
    RecMdcTrack* trk        = new RecMdcTrack;
    TTrack*      t          = _tracks[i];
    int          trackindex = i + nTdsTk; // for combination
 
    HitRefVec     hitrefvec;
    AList<TMLink> hits  = t->links();
    float         chisq = 0;
 
    HepPoint3D pos    = t->helix().pivot();
    int        charge = -1;
    HepVector  m_a( 5, 0 );
    m_a            = t->helix().a();
    int statfinder = t->getFinderType();
    if ( abs( statfinder ) > 1000 && brunge ) statfinder = 103;
    if ( abs( statfinder ) > 1000 && !brunge ) statfinder = 3;
    // be cautious
    if ( !brunge ) m_a[2] = -m_a[2];
    if ( abs( m_a[1] ) > 999 ) m_a[1] = 0;
    while ( m_a[1] < 0 ) { m_a[1] += 2 * pi; }
    while ( m_a[1] > 2 * pi ) { m_a[1] -= 2 * pi; }
    /// added by Jike Wang
    const double x0 = t->helix().pivot().x();
    const double y0 = t->helix().pivot().y();
 
    const double dr    = m_a[0];
    const double phi0  = m_a[1];
    const double kappa = m_a[2];
    const double dz    = m_a[3];
    const double tanl  = m_a[4];
 
    const double Bz    = -1000 * ( getPmgnIMF()->getReferField() );
    const double alpha = 333.564095 / Bz;
 
    const double cox = x0 + dr * cos( phi0 ) - alpha * cos( phi0 ) / kappa;
    const double coy = y0 + dr * sin( phi0 ) - alpha * sin( phi0 ) / kappa;
 
    unsigned nHits    = t->links().length();
    unsigned nStereos = 0;
    unsigned firstLyr = 44;
    unsigned lastLyr  = 0;
    for ( unsigned j = 0; j < nHits; j++ )
    {
 
      TMLink* l = hits[j];
 
      HepPoint3D ontrack = l->positionOnTrack();
      HepPoint3D onwire  = l->positionOnWire();
      HepPoint3D dir( ontrack.y() - coy, cox - ontrack.x(), 0 );
      double     pos_phi = onwire.phi();
      double     dir_phi = dir.phi();
      while ( pos_phi > pi ) { pos_phi -= pi; }
      while ( pos_phi < 0 ) { pos_phi += pi; }
      while ( dir_phi > pi ) { dir_phi -= pi; }
      while ( dir_phi < 0 ) { dir_phi += pi; }
      double entrangle = dir_phi - pos_phi;
      while ( entrangle > pi / 2 ) entrangle -= pi;
      while ( entrangle < ( -1 ) * pi / 2 ) entrangle += pi;
 
      // jialk setFltLen to tds
      int    imass   = 3;
      float  tl      = t->helix().a()[4];
      float  f       = sqrt( 1. + tl * tl );
      float  s       = fabs( t->helix().curv() ) * fabs( l->dPhi() ) * f;
      float  p       = f / fabs( t->helix().a()[2] );
      float  Mass[5] = { 0.000511, 0.105, 0.140, 0.493677, 0.98327 };
      double tof = s / 29.98 * sqrt( 1.0 + ( Mass[imass - 1] / p ) * ( Mass[imass - 1] / p ) );
 
      RecMdcHit* hit = new RecMdcHit;
      hit->setId( l->hit()->reccdc()->id );
      //       hit->setTrkId(i);
      hit->setTrkId( trackindex ); // for combination
      hit->setDriftDistLeft( l->drift( 0 ) );
      hit->setDriftDistRight( l->drift( 1 ) );
      hit->setErrDriftDistLeft( l->dDrift( 0 ) );
      hit->setErrDriftDistRight( l->dDrift( 1 ) );
      hit->setChisqAdd( l->pull() );
      hit->setFlagLR( l->leftRight() );
      //  hit->setStat(l->hit()->state());
      hit->setStat( 1 );
      //  std::cout<<"state :"<< l->hit()->state() << std::endl;
      hit->setAdc( l->hit()->reccdc()->adc );
      //           hit->setTdc((l->hit()->reccdc()->tdc + t0bes)/0.09375); //jialk
      hit->setTdc( l->hit()->reccdc()->timechannel );
      hit->setFltLen( tof * 30 ); // jialk
      // std::cout<<"tdc :"<< l->hit()->reccdc()->tdc <<"  "<<"t0bes : "<< t0bes <<std::endl;
      if ( cal ) { hit->setDoca( l->distancenew() ); }
      else { hit->setDoca( l->distance() ); }
      hit->setZhit( l->positionOnTrack().z() );
 
      /// add by jialk: returns driftTime prop time correction & entra angle
      hit->setEntra( entrangle );
      hit->setDriftT( l->getDriftTime() );
 
      Identifier hitIdf = MdcID::wire_id( l->wire()->layerId(), l->wire()->localId() );
      hit->setMdcId( hitIdf );
      hitcol->push_back( hit );
 
      SmartRef<RecMdcHit> refhit( hit );
      hitrefvec.push_back( refhit );
      chisq += l->pull();
      if ( l->wire()->stereo() ) ++nStereos;
      if ( firstLyr > l->wire()->layerId() ) firstLyr = l->wire()->layerId();
      if ( lastLyr < l->wire()->layerId() ) lastLyr = l->wire()->layerId();
    }
 
    HepSymMatrix m_ea = t->helix().Ea();
    double       errorMat[15];
    int          k = 0;
    for ( int ie = 0; ie < 5; ie++ )
    {
      for ( int je = ie; je < 5; je++ )
      {
        errorMat[k] = m_ea[ie][je];
        k++;
      }
    }
 
    trk->setTrackId( trackindex );
    trk->setHelix( m_a );
    trk->setPxy( t->pt() );
    trk->setPx( t->pt() * ( -sin( t->helix().phi0() ) ) );
    trk->setPy( t->pt() * cos( t->helix().phi0() ) );
    trk->setPz( t->pz() );
    trk->setP( t->ptot() );
 
    // jialk
    double theta = acos( ( t->pz() ) / t->ptot() );
    trk->setTheta( theta );
 
    double poca_dr   = t->helix().dr();
    double poca_phi0 = t->helix().phi0();
    trk->setPhi( poca_phi0 );
    HepPoint3D poca( pos.x() + poca_dr * cos( poca_phi0 ),
                     pos.y() + poca_dr * sin( poca_phi0 ), pos.z() + t->helix().dz() );
    trk->setPoca( poca );
    trk->setX( poca.x() );
    trk->setY( poca.y() );
    trk->setZ( poca.z() );
    trk->setR( sqrt( poca.x() * poca.x() + poca.y() * poca.y() ) );
    trk->setPivot( pos );
    trk->setVX0( pos.x() );
    trk->setVY0( pos.y() );
    trk->setVZ0( pos.z() );
 
    trk->setError( m_ea );
    //    trk->setError(errorMat);   //...............2
 
    //    double poca_dr = t->helix().dr();
    //    double poca_phi0 = t->helix().phi0();
    //    HepPoint3D poca(poca_dr*cos(poca_phi0),poca_dr*sin(poca_phi0),t->helix().dz());
    //    trk->setVX0(pos.x()+poca_dr*cos(poca_phi0));
    //    trk->setVY0(pos.y()+poca_dr*sin(poca_phi0));
    //    trk->setVZ0(pos.z()+t->helix().dz());
    //    cout<<"pivot:("<<pos.x()<<","<<pos.y()<<","<<pos.z()<<")"<<endl;
    //    cout<<"poca:("<<trk->getVX0()<<","<<trk->getVY0()<<","<<trk->getVZ0()<<")"<<endl;
 
    trk->setChi2( chisq );
    trk->setNdof( nHits - 5 );
    if ( t->quality() & TrackQuality2D ) trk->setNdof( nHits - 3 );
 
    TMLink* last = OuterMost( hits );
    t->approach( *last );
    trk->setFiTerm( last->dPhi() );
 
    trk->setNhits( nHits );
    trk->setNster( nStereos );
    trk->setStat( statfinder ); // yzhang add stat: ConformalFinder=2, CurlFiner=3
    // xuqn 2013-02-26
    // trk->setCharge(int(t->charge())*(-1));
    if ( !brunge ) trk->setCharge( int( t->charge() ) * ( -1 ) );
    else trk->setCharge( t->charge() );
    trk->setVecHits( hitrefvec );
    trk->setFirstLayer( firstLyr );
    trk->setLastLayer( lastLyr );
    // cout<<"first: "<<firstLyr<<"  last: "<<lastLyr<<endl;
    trkcol->push_back( trk );
  }
 
  // StatusCode trksc;
  // trksc = eventSvc->registerObject("/Event/Recon/RecMdcTrackCol", trkcol);
  // if( trksc.isFailure() ) {
  //   log << MSG::FATAL << "Could not register MdcTrack" << endmsg;
  //   return trksc;
  // }
  //
  // StatusCode hitsc = eventSvc->registerObject("/Event/Recon/RecMdcHitCol", hitcol);
  // if ( hitsc.isFailure() ) {
  //   log << MSG::FATAL << "Could not register MdcRecHit" << endmsg;
  //   return hitsc;
  // }
  // log << MSG::INFO << "MdcTrackCol registered successfully!" <<endmsg;
 
  /// check the result:MdcTrackCol
#ifdef TRKRECO_DEBUG
  SmartDataPtr<RecMdcTrackCol> newtrkCol( eventSvc, "/Event/Recon/RecMdcTrackCol" );
  if ( !newtrkCol )
  {
    log << MSG::FATAL << "Could not find RecMdcTrackCol" << endmsg;
    return ( StatusCode::FAILURE );
  }
  log << MSG::INFO << "Begin to check RecMdcTrackCol" << endmsg;
  std::cout << std::fixed << std::setprecision( 6 ); // yzhang debug for 720
  RecMdcTrackCol::iterator iter_trk = newtrkCol->begin();
  for ( ; iter_trk != newtrkCol->end(); iter_trk++ )
  {
    const RecMdcTrack* tk = *iter_trk;
    std::cout << "//==== " << name() << " print RecMdcTrack No." << tk->trackId() << " :"
              << std::endl;
    cout << " dr " << tk->helix( 0 ) << " phi0 " << tk->helix( 1 ) << " cpa " << tk->helix( 2 )
         << " z0 " << tk->helix( 3 ) << " tanl " << tk->helix( 4 ) << endl;
    bool printTrackDetail = true;
    if ( printTrackDetail )
    {
      std::cout << " q " << tk->charge() << " theta " << tk->theta() << " phi " << tk->phi()
                << " x0 " << tk->x() << " y0 " << tk->y() << " z0 " << tk->z() << " r0 "
                << tk->r() << endl;
      std::cout << " p " << tk->p() << " pt " << tk->pxy() << " pxyz(" << tk->px() << ","
                << tk->py() << "," << tk->pz() << ")" << endl;
      std::cout << " tkStat " << tk->stat() << " chi2 " << tk->chi2() << " ndof " << tk->ndof()
                << " nhit " << tk->getNhits() << " nst " << tk->nster() << endl;
      bool printErrMat = true;
      if ( printErrMat )
      {
        std::cout << "errmat   " << std::endl;
        for ( int i = 0; i < 15; i++ ) { std::cout << " " << tk->err( i ); }
        std::cout << "   " << std::endl;
      }
    }
 
    bool printHit = true;
    if ( printHit )
    {
      int  haveDigi[43][288];
      bool printMcTk = true;
      if ( printMcTk )
      {
        for ( int ii = 0; ii < 43; ii++ )
        {
          for ( int jj = 0; jj < 43; jj++ ) { haveDigi[ii][jj] = -9999; }
        }
        MdcDigiVec           mdcDigiVec = getRawDataProviderSvc()->getMdcDigiVec();
        MdcDigiCol::iterator iter       = mdcDigiVec.begin();
        std::cout << name() << "//==== " << name()
                  << " print MdcDigiVec, nDigi=" << mdcDigiVec.size() << " :" << std::endl;
        for ( int iDigi = 0; iter != mdcDigiVec.end(); iter++, iDigi++ )
        {
          long l         = MdcID::layer( ( *iter )->identify() );
          long w         = MdcID::wire( ( *iter )->identify() );
          haveDigi[l][w] = ( *iter )->getTrackIndex();
        }
      }
 
      int nhits = tk->getVecHits().size();
      std::cout << "nHits=" << nhits << std::endl;
      cout << "No. ";
      if ( printMcTk ) cout << "mcTkId ";
      cout << "(layer,wire,lr) stat z" << endl;
      for ( int ii = 0; ii < nhits; ii++ )
      {
        Identifier id( tk->getVecHits()[ii]->getMdcId() );
        int        layer = MdcID::layer( id );
        int        wire  = MdcID::wire( id );
        cout << ii << ":";
        if ( printMcTk ) { cout << haveDigi[layer][wire]; }
        cout << " (" << layer << "," << wire << "," << tk->getVecHits()[ii]->getFlagLR()
             << ")  " << tk->getVecHits()[ii]->getStat() << "  "
             << tk->getVecHits()[ii]->getZhit() << " " << std::endl;
      } // end of hit list
      std::cout << "  " << std::endl;
    }
    /*
       std::cout << "retrieved MDC track:"
       << "Track Id: " << (*iter_trk)->getTrkId()
       << "   Pivot: " << (*iter_trk)->getVX0() << " "
       << (*iter_trk)->getVY0() << " " << (*iter_trk)->getVZ0()
       << endl
       << "Phi0: " << (*iter_trk)->getFi0() << "   Error of Phi0 "
       << (*iter_trk)->getError()(2,2) << endmsg
       << "kappa: " << (*iter_trk)->getCpa() << "   Error of kappa "
       << (*iter_trk)->getError()(3,3) << endmsg
       << "Tanl: " << (*iter_trk)->getTanl() << "   Error of Tanl "
       << (*iter_trk)->getError()(5,5) << endmsg
       << "Chisq of fit: "<< (*iter_trk)->getChisq()
       << "   Phi terminal: "<< (*iter_trk)->getFiTerm()
       << endl
       << "Number of hits: "<< (*iter_trk)->getNhits()
       << "   Number of stereo hits " << (*iter_trk)->getNster()
       << endmsg;
 
       log << MSG::DEBUG << "hitList of this  track:" << endmsg;
       std::vector<MdcRecHit*> gothits = (*iter_trk)->getVecHits();
       std::vector<MdcRecHit*>::iterator it_gothit = gothits.begin();
       for( ; it_gothit != gothits.end(); it_gothit++){
       log << MSG::DEBUG << "hits Id: "<<(*it_gothit)->getId()
       << "   hits DDL&DDR " <<(*it_gothit)->getDriftDistLeft()
       << "   hits MDC IDentifier " <<(*it_gothit)->getMdcId()
       << endmsg
       << "   hits TDC " <<(*it_gothit)->getTdc()
       << "   hits ADC " <<(*it_gothit)->getAdc() << endmsg;
       }
       */
  }
  std::cout << std::defaultfloat;
#endif
  return StatusCode::SUCCESS;
}
 
void TTrackManager::saveTables( void ) {
#ifdef TRKRECO_DEBUG
  std::cout << "TTrackManager::saveTables ... # 3D tracks=" << _tracks.length()
            << ", # 2D tracks=" << _tracks2D.length() << ", all tracks=" << _tracksAll.length()
            << std::endl;
#endif
 
  //...For 3D tracks...
  AList<TTrack> badTracks;
  unsigned      n  = _tracks.length();
  unsigned*     id = (unsigned*)malloc( n * sizeof( unsigned ) );
  //    bzero((char *) id, n * sizeof(unsigned));
  memset( (char*)id, 0, n * sizeof( unsigned ) );
  for ( unsigned i = 0; i < n; i++ )
  {
    TTrack& t = *_tracks[i];
    if ( !t.nLinks() )
    {
      badTracks.append( (TTrack&)t );
      continue;
    }
 
    //...Copy track parameters...
    MdcRec_trk*     r   = 0;
    MdcRec_trk_add* a   = 0;
    int             err = copyTrack( t, &r, &a );
    if ( err )
    {
      badTracks.append( t );
      continue;
    }
    _tracksFinal.append( t );
 
    //...Type and quality...
    id[i]       = r->id;
    r->stat     = t.state();
    a->kind     = t.type();
    a->decision = t.finder();
    a->stat     = t.fitting();
    //      if (a->m_kind == TrackTypeCosmic) {
    //          a->m_quality = TrackQualityCosmic;
    //      }
    //      if (t.daughter() && (_tracks.index(t.daughter()) >= 0))
    //          a->m_daughter = _tracks.index(t.daughter()) + 1;
  }
 
  //...Daughter treatment...
  for ( unsigned i = 0; i < n; i++ )
  {
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "id[" << i << "]=" << id[i] << std::endl;
#endif
    if ( !( id[i] ) ) continue;
    if ( !( _tracks[i]->daughter() ) ) continue;
 
    int dId = _tracks.index( _tracks[i]->daughter() );
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    dId=" << dId;
    if ( dId >= 0 ) std::cout << ", id[dId]=" << id[dId];
    std::cout << std::endl;
#endif
 
    if ( dId >= 0 )
    {
      if ( id[dId] )
      {
        //      reccdc_trk_add * a = (reccdc_trk_add *)
        //          BsGetEnt(RECMDC_TRK_ADD, id[i], BBS_No_Index);
        //      a->m_daughter = id[dId];
        MdcRec_trk_add* a = &( *MdcRecTrkAddCol::getMdcRecTrkAddCol() )[i];
      }
    }
  }
  free( id );
 
  //...Remove bad tracks...
  _tracks.remove( badTracks );
  badTracks.removeAll();
 
  //...For 2D tracks...
  n = _tracks2D.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TTrack& t = *_tracks2D[i];
 
    //...Copy track parameters...
    MdcRec_trk*     r   = 0;
    MdcRec_trk_add* a   = 0;
    int             err = copyTrack( t, &r, &a );
    if ( err )
    {
#ifdef TRKRECO_DEBUG
      std::cout << "TTrackManager::saveTables !!! bad 2D tracks found"
                << " : err=" << err << std::endl
                << TrackDump( t ) << std::endl;
#endif
      badTracks.append( t );
      continue;
    }
    _tracksFinal.append( t );
 
    //...Reset helix parameter...
    //      r->m_helix[3] = 0.;
    //      r->m_helix[4] = 0.;
    //      r->m_nhits -= r->m_nster;
    //      r->m_nster = 0;
 
    //...Table filling...
    r->stat     = t.state();
    a->kind     = t.type();
    a->decision = t.finder();
    // a->m_quality = t.quality();
    a->quality = TrackQuality2D;
    a->stat    = t.fitting();
 
#ifdef TRKRECO_DEBUG
    if ( ( r->ndf == 0 ) && ( r->chiSq > 0. ) )
    {
      std::cout << "TTrackManager::saveTables !!! chisq>0 with ndf=0." << std::endl
                << "    Here is a track dump"
                << "    " << TrackDump( t ) << std::endl;
      t.dump( "detail" );
    }
    if ( ( r->ndf > 0 ) && ( r->chiSq == 0. ) )
    {
      std::cout << "TTrackManager::saveTables !!! chisq=0 with ndf>0." << std::endl
                << "    Here is a track dump"
                << "    " << TrackDump( t ) << std::endl;
      t.dump( "detail" );
    }
 
    if ( r->ndf == 0 )
      std::cout << "TTrackManager::saveTables ... ndf = 0" << std::endl
                << "    " << TrackDump( t ) << std::endl;
    if ( r->chiSq == 0. )
      std::cout << "TTrackManager::saveTables ... chisq = 0" << std::endl
                << "    " << TrackDump( t ) << std::endl;
#endif
  }
  _tracks2D.remove( badTracks );
 
  //...Statistics...
  ++_s->_nEvents;
  _s->_nTracks += _tracks.length();
  _s->_nTracksAll += _tracksAll.length();
  _s->_nTracks2D += _tracks2D.length();
  _s->_nTracksFinal += _tracksFinal.length();
}
 
void TTrackManager::saveMCTables( void ) const {
  unsigned n = _tracksFinal.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    const TTrack& t = *_tracksFinal[i];
 
    //  struct reccdc_trk * r;
    //  r = (struct reccdc_trk *) BsGetEnt(RECMDC_TRK, i + 1, BBS_No_Index);
    MdcRec_trk* r = &( *MdcRecTrkCol::getMdcRecTrkCol() )[i];
 
    //...Set type...
 
    //...Hit loop...
    const AList<TMLink>& hits  = t.finalHits();
    unsigned             nHits = hits.length();
    for ( unsigned j = 0; j < nHits; j++ )
    {
      TMLink*          l = hits[j];
      MdcRec_wirhit*   h = l->hit()->reccdc();
      MdcDat_mcwirhit* m = l->hit()->mc()->datcdc();
      m->trk             = r;
      //        struct reccdc_mctrk2hep * c;
      //        c = (struct reccdc_mctrk2hep *) BsNewEnt(RECMDC_MCTRK2HEP);
      MdcRec_mctrk2hep* c = new MdcRec_mctrk2hep;
      MdcRecMctrk2hepCol::getMdcRecMctrk2hepCol()->push_back( *c );
      c->wir = h;
      c->trk = r;
      c->hep = l->hit()->mc()->hep()->gen();
    }
 
    const TTrackMC* const mc = t.mc();
    //  struct reccdc_mctrk * m;
    //  m = (struct reccdc_mctrk *) BsNewEnt(RECMDC_MCTRK);
    // //   MdcRec_mctrk* m = &(*MdcRecMctrkCol::getMdcRecMctrkCol())[0];
    MdcRec_mctrk* m = new MdcRec_mctrk;
    MdcRecMctrkCol::getMdcRecMctrkCol()->push_back( *m );
    m->wirFrac    = mc->wireFraction();
    m->wirFracHep = mc->wireFractionHEP();
    m->charge     = mc->charge();
    m->ptFrac     = mc->ptFraction();
    m->pzFrac     = mc->pzFraction();
    m->quality    = mc->quality();
    if ( mc->hep() ) m->hep = mc->hep()->gen();
    else m->hep = 0;
  }
}
 
void TTrackManager::movePivot( void ) {
  unsigned n = _tracksAll.length();
  //  unsigned n = _tracks.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    if ( _tracksAll[i]->nLinks() ) _tracksAll[i]->movePivot();
    //  if (_tracks[i]->nLinks()) _tracks[i]->movePivot();
  }
  nameTracks();
}
 
void TTrackManager::clear( void ) {
  HepAListDeleteAll( _tracksAll );
  _tracks.removeAll();
  _tracks2D.removeAll();
  _tracksFinal.removeAll();
  HepAListDeleteAll( _associateHits );
  static bool first = true;
  if ( first )
  {
    first = false;
    int size;
    //  _s = (struct summary *)
    //      BASF_Sharedmem->get_pointer(BASF_Sharedmem->get_id(),
    //                  "TrkMgrSum",
    //                  & size);
  }
}
 
void TTrackManager::finish( void ) {
  refit();
  movePivot();
  if ( _debugLevel > 1 )
  {
    std::cout << name() << " ... finishing" << std::endl;
    //  unsigned n = _tracksAll.length();
    unsigned n = _tracks.length();
    for ( unsigned i = 0; i < n; i++ )
    {
      //        TTrack & t = * _tracksAll[i];
      TTrack& t = *_tracks[i];
      std::cout << "    " << t.name() << std::endl;
      t.dump( "hits mc track flag sort", "    " );
    }
  }
}
 
void TTrackManager::append( AList<TTrack>& list ) {
  _tracksAll.append( list );
  _tracks.append( selectGoodTracks( list ) );
  list.removeAll();
}
 
void TTrackManager::append2D( AList<TTrack>& list ) {
  _tracksAll.append( list );
  _tracks2D.append( selectGoodTracks( list, true ) );
  _tracks2D.sort( SortByPt );
  list.removeAll();
}
 
void TTrackManager::refit( void ) {
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << name() << " ... refitting" << std::endl;
#endif
  unsigned      n = _tracks.length();
  AList<TTrack> bads;
  for ( unsigned i = 0; i < n; i++ )
  {
    TTrack& t = *_tracks[i];
    int     err;
    err = _fitter.fit( t );
    if ( err < 0 )
    {
      bads.append( t );
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout << "    " << t.name();
      std::cout << " rejected because of fitting failure" << std::endl;
#endif
      continue;
    }
    t.refine( 30. * 10. );
    err = _fitter.fit( t );
    if ( err < 0 )
    {
      bads.append( t );
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout << "    " << t.name();
      std::cout << " rejected because of fitting failure" << std::endl;
#endif
      continue;
    }
    t.refine( 30. * 1. );
    err = _fitter.fit( t );
    if ( err < 0 )
    {
      bads.append( t );
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout << "    " << t.name();
      std::cout << " rejected because of fitting failure" << std::endl;
#endif
      continue;
    }
  }
  _tracks.remove( bads );
}
 
void TTrackManager::mask( void ) const {
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << name() << " ... masking" << std::endl;
#endif
 
  unsigned n = _tracks.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TTrack& t = *_tracks[i];
 
    //...Skip if no core...
    //   This should not be happend...
    if ( !t.cores().length() ) continue;
 
    //...Counts # of hits per layer...
    unsigned nHits[43];
    NHits( t.cores(), nHits );
 
    //...Check each layer...
    bool needMask = false;
    for ( unsigned j = 0; j < 43; j++ )
    {
      if ( nHits[j] > 1 )
      {
        AList<TMLink> linksInLayer = SameLayer( t.links(), j );
        if ( Width( linksInLayer ) > 2 )
        {
          needMask = true;
 
#ifdef TRKRECO_DEBUG_DETAIL
          Dump( linksInLayer, "sort", "    -->" );
#endif
          break;
        }
      }
    }
    if ( !needMask ) continue;
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    trk" << i << "(id is tmp) needs mask" << std::endl;
    std::cout << "        type = " << t.type() << std::endl;
#endif
 
    //...Switch by track type...
    switch ( t.type() )
    {
    case TrackTypeNormal:
      maskNormal( t );
      maskMultiHits( t );
      break;
    case TrackTypeCurl:
      maskCurl( t );
      maskMultiHits( t );
      break;
    default: break;
    }
 
    //...Refit...
    // refit() ???
    _fitter.fit( t );
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    masking result : ";
    t.dump( "detail sort", "    " );
#endif
  }
}
 
void TTrackManager::nameTracks( void ) {
  unsigned n = _tracks.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TTrack& t = *_tracks[i];
    t._name   = "trk" + std::to_string( i ) + "(" + t._name + ")";
  }
  AList<TTrack> tmp = _tracksAll;
  tmp.remove( _tracks );
  unsigned n1 = tmp.length();
  for ( unsigned i = 0; i < n1; i++ )
  {
    TTrack& t = *tmp[i];
    t._name   = "trk" + std::to_string( i + n ) + "(" + t._name + ")";
  }
}
 
TMLink& TTrackManager::divide( const TTrack& t, AList<TMLink>* l ) const {
  TMLink&           start  = *OuterMost( t.links() );
  const HepPoint3D& center = t.helix().center();
  const HepVector3D a      = start.positionOnTrack() - center;
  for ( unsigned j = 0; j < t.nLinks(); j++ )
  {
    if ( t[j] == &start ) continue;
    TMLink&           k = *t[j];
    const HepVector3D b = k.positionOnTrack() - center;
    if ( a.cross( b ).z() >= 0. ) l[0].append( k );
    else l[1].append( k );
  }
 
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    outer link = " << start.hit()->wire()->name() << std::endl;
  std::cout << "        nLinks of 0 = " << l[0].length() << std::endl;
  std::cout << "        nLinks of 1 = " << l[1].length() << std::endl;
#endif
 
  if ( l[0].length() == 0 || l[1].length() == 0 ) return divideByIp( t, l );
 
  return start;
}
 
TMLink& TTrackManager::divideByIp( const TTrack& t, AList<TMLink>* l ) const {
  l[0].removeAll();
  l[1].removeAll();
 
  const HepPoint3D& center = t.helix().center();
  const HepVector3D a      = ORIGIN - center;
  for ( unsigned j = 0; j < t.nLinks(); j++ )
  {
    TMLink&           k = *t[j];
    const HepVector3D b = k.positionOnTrack() - center;
    if ( a.cross( b ).z() >= 0. ) l[0].append( k );
    else l[1].append( k );
  }
 
  //...This is a dummy...
  TMLink& start = *OuterMost( t.links() );
  return start;
}
 
void TTrackManager::removeHitsAcrossOverIp( AList<TMLink>& l ) const {
 
  //...Calculate average phi...
  unsigned n      = l.length();
  float    phiSum = 0.;
  for ( unsigned i = 0; i < n; i++ )
  {
    const TMDCWire& w     = *l[i]->hit()->wire();
    unsigned        j     = w.localId();
    unsigned        nWire = w.layer()->nWires();
 
    float phi = (float)j / (float)nWire;
    phiSum += phi;
  }
  float average = phiSum / (float)n;
 
  AList<TMLink> cross;
  for ( unsigned i = 0; i < n; i++ )
  {
    const TMDCWire& w     = *l[i]->hit()->wire();
    unsigned        j     = w.localId();
    unsigned        nWire = w.layer()->nWires();
 
    float phi = (float)j / (float)nWire;
    float dif = fabs( phi - average );
    if ( dif > 0.5 ) dif = 1. - dif;
 
    if ( dif > 0.3 ) cross.append( l[i] );
  }
  l.remove( cross );
 
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    Cross over IP reduction : ";
  for ( unsigned i = 0; i < cross.length(); i++ )
  { std::cout << cross[i]->wire()->name() << ","; }
  std::cout << std::endl;
#endif
}
 
void TTrackManager::maskOut( TTrack& t, const AList<TMLink>& links ) const {
  unsigned n = links.length();
  if ( !n ) return;
  for ( unsigned i = 0; i < n; i++ )
  {
    const TMDCWireHit& hit = *links[i]->hit();
    hit.state( hit.state() | WireHitInvalidForFit );
  }
  t._fitted = false;
 
#ifdef TRKRECO_DEBUG_DETAIL
  Dump( links, "detail", "    TTrackManager::maskOut ... masking " );
#endif
}
 
void TTrackManager::maskMultiHits( TTrack& t ) const {
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "... masking multi-hits" << std::endl;
#endif
 
  if ( !t.cores().length() ) return;
  AList<TMLink> cores        = t.cores();
  unsigned      n            = cores.length();
  bool          layerLimited = false;
  AList<TMLink> bads;
 
  cores.sort( SortByWireId );
  for ( unsigned i = 0; i < n; i++ )
  {
    if ( layerLimited )
    {
      bads.append( cores[i] );
      continue;
    }
    AList<TMLink> linksInLayer = SameLayer( cores, cores[i]->wire()->layerId() );
    if ( linksInLayer.length() > 3 )
    {
      bads.append( cores[i] );
      layerLimited = true;
    }
  }
  maskOut( t, bads );
}
 
void TTrackManager::maskNormal( TTrack& t ) const {
 
  //...Divide into two tracks...
  AList<TMLink> l[2];
  TMLink&       start = divideByIp( t, l );
 
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    normal : divided by IP" << std::endl;
  std::cout << "    0=";
  for ( unsigned j = 0; j < l[0].length(); j++ )
  { std::cout << "," << l[0][j]->wire()->name(); }
  std::cout << std::endl;
  std::cout << "    1=";
  for ( unsigned j = 0; j < l[1].length(); j++ )
  { std::cout << "," << l[1][j]->wire()->name(); }
  std::cout << std::endl;
#endif
 
  //...Which should be masked out ?...
  unsigned maskSide = 2;
 
  //...1. Check by # of super layers...
  if ( NSuperLayers( l[0] ) < NSuperLayers( l[1] ) ) maskSide = 0;
  else if ( NSuperLayers( l[0] ) > NSuperLayers( l[1] ) ) maskSide = 1;
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    NSuperLayers 0, 1 = " << NSuperLayers( l[0] ) << ", ";
  std::cout << NSuperLayers( l[1] ) << std::endl;
#endif
  if ( maskSide != 2 )
  {
    maskOut( t, l[maskSide] );
    return;
  }
 
  //...2. Check by the inner-most layer...
  unsigned i0 = InnerMost( l[0] )->wire()->layerId();
  unsigned i1 = InnerMost( l[1] )->wire()->layerId();
  if ( i0 < i1 ) maskSide = 1;
  else if ( i0 > i1 ) maskSide = 0;
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    i0, i1 = " << i0 << ", " << i1 << std::endl;
#endif
  if ( maskSide != 2 )
  {
    maskOut( t, l[maskSide] );
    return;
  }
 
  //...3. Check by # of layers...
  if ( NLayers( l[0] ) < NLayers( l[1] ) ) maskSide = 0;
  else if ( NLayers( l[0] ) > NLayers( l[1] ) ) maskSide = 1;
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    NLayers 0, 1 = " << NLayers( l[0] ) << ", ";
  std::cout << NLayers( l[1] ) << std::endl;
#endif
  if ( maskSide != 2 )
  {
    maskOut( t, l[maskSide] );
    return;
  }
 
  //...4. Check by pt...
  if ( maskSide == 2 )
  {
    TTrack* tt[2];
    for ( unsigned j = 0; j < 2; j++ )
    {
      tt[j] = new TTrack( t );
      tt[j]->remove( l[j] );
      _fitter.fit( *tt[j] );
    }
    if ( tt[1]->pt() > tt[0]->pt() ) maskSide = 1;
    else maskSide = 0;
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    pt 0 = " << tt[1]->pt() << std::endl;
    std::cout << "    pt 1 = " << tt[0]->pt() << std::endl;
#endif
    delete tt[0];
    delete tt[1];
  }
  maskOut( t, l[maskSide] );
  return;
}
 
void TTrackManager::maskCurl( TTrack& t ) const {
 
  //...Divide into two tracks...
  AList<TMLink> l[2];
  TMLink&       start = divideByIp( t, l );
  if ( l[0].length() == 0 ) return;
  if ( l[1].length() == 0 ) return;
 
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    curl : divided by IP" << std::endl;
  std::cout << "    0:";
  Dump( l[0], "flag sort" );
  std::cout << "    1:";
  Dump( l[1], "flag sort" );
  std::cout << std::endl;
#endif
 
  //...Which should be masked out ?...
  unsigned maskSide = 2;
 
  //...1. Check by # of super layers...
  if ( NSuperLayers( l[0] ) < NSuperLayers( l[1] ) ) maskSide = 0;
  else if ( NSuperLayers( l[0] ) > NSuperLayers( l[1] ) ) maskSide = 1;
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    NSuperLayers 0, 1 = " << NSuperLayers( l[0] ) << ", ";
  std::cout << NSuperLayers( l[1] ) << std::endl;
#endif
  if ( maskSide != 2 )
  {
    maskOut( t, l[maskSide] );
    return;
  }
 
  //...Make two tracks...
  TTrack* tt[2];
  tt[0] = new TTrack( t );
  tt[1] = new TTrack( t );
  tt[0]->remove( l[1] );
  tt[1]->remove( l[0] );
  _fitter.fit( *tt[0] );
  _fitter.fit( *tt[1] );
  Helix h0 = Helix( tt[0]->helix() );
  Helix h1 = Helix( tt[1]->helix() );
 
  //...Check by z...
  h0.pivot( ORIGIN );
  h1.pivot( ORIGIN );
  if ( fabs( h0.dz() ) < fabs( h1.dz() ) ) maskSide = 1;
  else maskSide = 0;
 
  delete tt[0];
  delete tt[1];
  maskOut( t, l[maskSide] );
  return;
}
 
StatusCode TTrackManager::determineT0( unsigned level, unsigned nMax ) {
#ifdef TRKRECO_DEBUG_DETAIL
  if ( level == 0 )
  {
    std::cout << "TTrackManager::determineT0 !!! called with level = 0";
    std::cout << std::endl;
  }
#endif
 
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "TTrackManager" );
 
  IDataProviderSvc* eventSvc = NULL;
  Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
 
  static bool     first   = true;
  static unsigned methode = 0;
  if ( first )
  {
    first = false;
 
    if ( level == 1 ) { _cFitter.fit2D( true ); }
    else if ( level == 2 )
    {
      // default setting
    }
    else if ( level == 3 )
    {
      //        _cFitter.sag(true);  //Liuqg
    }
    else if ( level == 4 )
    {
      //        _cFitter.sag(true);  //Liuqg
      _cFitter.propagation( true );
    }
    else if ( level == 5 )
    {
      //        _cFitter.sag(true);  //Liuqg
      _cFitter.propagation( true );
      _cFitter.tof( true );
    }
    else if ( level == 6 )
    {
      methode = 1;
      //        _cFitter.sag(true); //Liuqg
      _cFitter.propagation( true );
      _cFitter.tof( true );
    }
    else if ( level == 7 )
    {
      methode = 2;
      //        _cFitter.sag(true);  //Liuqg
      _cFitter.propagation( true );
      _cFitter.tof( true );
    }
  }
 
  unsigned n = _tracks.length();
  if ( !n ) return StatusCode::SUCCESS;
 
  if ( nMax == 0 ) nMax = n;
  if ( n > nMax ) n = nMax;
 
  //    float t0 = 0.;
  _t0 = 999.;
 
  // read t0 from TDS
  float                      t0_sta = -1;
  float                      tev    = 0;
  SmartDataPtr<RecEsTimeCol> aevtimeCol( eventSvc, "/Event/Recon/RecEsTimeCol" );
  if ( aevtimeCol )
  {
    RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
    t0_sta                          = ( *iter_evt )->getStat();
    tev                             = ( *iter_evt )->getTest();
    //  cout<<"t0_sta: "<<t0_sta<<"  tev: "<<tev<<endl;
  }
  else { log << MSG::WARNING << "Could not find RecEsTimeCol" << endmsg; }
 
  if ( methode == 0 ) _t0 = T0( n );
 
  else if ( methode == 1 ) _t0 = T0Fit( n );
 
  else if ( methode == 2 )
  { // revise method == 2 to BESIII environment. Liuqg
    if ( t0_sta != 1 )
    { // 1: tof  11:tof after reco  //no tof
      _t0 = T0Fit( n );
      // cout<<"t0: "<<_t0<<endl;
    }
  }
 
  // std::cout << "reccdc_timing=" << BsCouTab(RECMDC_TIMING) << std::endl;
  /*    else if (methode == 2 && BsCouTab(RECMDC_TIMING) != 0) {
        struct reccdc_timing * r0 = (struct reccdc_timing *)
        BsGetEnt(RECMDC_TIMING, BsCouTab(RECMDC_TIMING), BBS_No_Index);
        if (r0->m_quality == 102) {
        if (BsCouTab(BELLE_EVENT)) {
        struct belle_event * b0 = (struct belle_event *)
        BsGetEnt(BELLE_EVENT, 1, BBS_No_Index);
        if(1==b0->m_ExpMC)                    t0 = T0Fit(n);
        if(2==b0->m_ExpMC && r0->m_time !=0.) t0 = T0Fit(n);
        }
        }
        else if (r0->m_quality == 100) t0 = T0Fit(n);
  //      std::cout << "quality=" << r0->m_quality << std::endl;
  }
  */
 
  //...For debug...
  if ( _debugLevel )
  {
    std::cout << "TTrackManager::determineT0 ... methode=" << methode;
    std::cout << ", T0 offset=" << -_t0;
    std::cout << ", # of tracks used=" << n << std::endl;
  }
 
  //...store them... Liuqg
  float t0_rec    = 999.;
  int   t0_recSta = 8;
  if ( fabs( _t0 + tev ) < 4 ) t0_rec = 0;
  if ( fabs( _t0 + tev - 8 ) < 4 ) t0_rec = 8;
  if ( fabs( _t0 + tev - 16 ) < 4 ) t0_rec = 16;
  log << MSG::INFO << "beginning to make RecEsTimeCol" << endmsg;
  IDataManagerSvc* dataManSvc = dynamic_cast<IDataManagerSvc*>( eventSvc );
  DataObject*      aEvTime;
  eventSvc->findObject( "/Event/Recon/RecEsTimeCol", aEvTime );
  if ( aEvTime != NULL && t0_rec < 25 )
  {
    dataManSvc->clearSubTree( "/Event/Recon/RecEsTimeCol" );
    eventSvc->unregisterObject( "/Event/Recon/RecEsTimeCol" );
  }
  if ( t0_rec < 25 )
  {
 
    RecEsTimeCol* aEvTimeCol = new RecEsTimeCol;
    RecEsTime*    aevtime    = new RecEsTime;
    aevtime->setTest( t0_rec );
    aevtime->setStat( t0_recSta );
    aEvTimeCol->push_back( aevtime );
 
    // register event time  to TDS
    // check whether the t0  has been already registered
    StatusCode evtime = eventSvc->registerObject( "/Event/Recon/RecEsTimeCol", aEvTimeCol );
 
    if ( evtime != StatusCode::SUCCESS )
    {
      log << MSG::FATAL << "Could not register Event Start Time" << endmsg;
      return ( StatusCode::FAILURE );
    }
  }
  return ( StatusCode::SUCCESS );
}
 
float TTrackManager::T0( unsigned n ) {
 
#define X0 -10.
#define X1 0.
#define X2 10.
#define STEP 10.
 
  //...Determine T0 for each track...
  float t0Sum = 0.;
  for ( unsigned i = 0; i < n; i++ )
  {
    TTrack& t = *_tracks[i];
    float   y[3];
    for ( unsigned j = 0; j < 3; j++ )
    {
      float offset = X0 + j * STEP;
      _cFitter.fit( t, offset );
      y[j] = t.chi2();
    }
    t0Sum += minimum( y[0], y[1], y[2] );
  }
  float t0 = t0Sum / (float)n;
  if ( isnan( t0 ) ) t0 = 0.;
 
  //...Fit with T0 correction...
  n = _tracks.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TTrack& t = *_tracks[i];
    _cFitter.fit( t, t0 );
  }
 
  //...Store it...
  //    reccdc_timing * t = (reccdc_timing *) BsNewEnt(RECMDC_TIMING);
  /*    MdcRec_timing* t = new MdcRec_timing;
        MdcRecTimingCol::getMdcRecTimingCol()->push_back(*t);
        t->time = - t0;
        t->quality = 11;
        */
  return -t0;
}
 
float TTrackManager::T0Fit( unsigned n ) {
 
  float tev_err;
  float tev_sum0 = 0.;
  float tev_sum  = 0.;
  float tev_sum2 = 0.;
  float w_sum    = 0.;
 
  // sort in order of pt
  //     std::cout << "length=" << _tracks.length() << std::endl;
  const int cn        = _tracks.length();
  float*    sort      = new float[cn];
  float     ptmax_pre = 1.e10;
  for ( unsigned i = 0; i < cn; i++ )
  {
    float ptmax = -999.;
    int   jmax = 0; // add initialize 25-05-15
    for ( unsigned j = 0; j < cn; j++ )
    {
      TTrack& tj = *_tracks[j];
      float   pt = fabs( 1. / tj.helix().a()[2] );
      if ( pt < ptmax_pre && pt > ptmax )
      {
        ptmax = pt;
        jmax  = j;
      }
      sort[i] = jmax;
    }
    ptmax_pre = ptmax;
  }
 
  for ( unsigned i = 0; i < n; i++ )
  {
    // srtbypt  TTrack & t = * _tracks[i];
    TTrack& t = *_tracks[sort[i]];
    // float pt = fabs(1./t.helix().a()[2]);
    // std::cout << "pt=" << pt << endl;
    float tev = 0.;
    _cFitter.fit( t, tev, tev_err );
    //  std::cout << "tev,tev_err=" <<tev<< " "<<tev_err<<endl;
    float w = 1. / tev_err / tev_err;
    tev_sum += w * tev;
    w_sum += w;
    // tev_sum0 += tev;
    // tev_sum2 += tev * tev;
  }
 
  delete[] sort;
 
  float tev_mean = tev_sum / w_sum;
  // float tev_err_a = 1. / sqrt(w_sum);
  // float tev_err_b = (tev_sum2 - tev_sum0 * tev_sum0 / (n + 1)) / n;
  // tev_err_b = sqrt(tev_err_b);
  // std::cout << "comp,t0,mean tev,err ="<<t0<<" "<<tev_mean<<" "<<tev_err_a<<"
  // "<<tev_err_b<<std::endl;
  if ( isnan( tev_mean ) ) tev_mean = 0.;
 
  //...Store it...
  //    reccdc_timing * tt = (reccdc_timing *) BsNewEnt(RECMDC_TIMING);
  //    MdcRec_timing* tt = new MdcRec_timing;
  //    MdcRecTimingCol::getMdcRecTimingCol()->push_back(*tt);
  //    tt->time = tev_mean;
  //    tt->quality = 151;
 
  return -tev_mean;
}
 
float TTrackManager::minimum( float y0, float y1, float y2 ) const {
  float xMin = X1 + 0.5 * STEP * ( y0 - y2 ) / ( y0 + y2 - 2. * y1 );
  return xMin;
}
 
// added by matsu ( 1999/05/24 )
void TTrackManager::merge( void ) {
 
  //...Merging...
  unsigned n = _tracks.length();
  // cout<<"tracks: "<<n<<endl;
  AList<TTrack> bads;
  unsigned*     flagTrk = (unsigned*)malloc( n * sizeof( unsigned ) );
  for ( unsigned i = 0; i < n; i++ ) flagTrk[i] = 0;
 
  //...Search a track to be merged...
  for ( unsigned i0 = 0; i0 < n; i0++ )
  {
 
    if ( flagTrk[i0] != 0 ) continue;
    TTrack& t0 = *_tracks[i0];
    if ( !( t0.pt() < 0.13 ) ) continue;
 
    unsigned Noverlap( 0 ), Nall( 0 );
    float    OverlapRatioMax( -999. );
    unsigned MaxID( 0 );
 
    for ( unsigned i1 = 0; i1 < n; i1++ )
    {
 
      if ( i0 == i1 || flagTrk[i1] != 0 ) continue;
      TTrack& t1 = *_tracks[i1];
      if ( !( t1.pt() < 0.13 ) ) continue;
      Nall = t1.links().length();
      if ( !Nall ) continue;
 
      Noverlap = 0;
      for ( unsigned j = 0; j < Nall; j++ )
      {
        TMLink&            l    = *t1.links()[j];
        const TMDCWireHit& whit = *l.hit();
        double             load( 3. );                 // jialk original is 2
        if ( whit.state() & WireHitStereo ) load = 4.; // jialk original is 3
 
        double x = t0.helix().center().x() - l.positionOnTrack().x();
        double y = t0.helix().center().y() - l.positionOnTrack().y();
        double r = sqrt( x * x + y * y );
        double R = fabs( t0.helix().radius() );
 
        if ( ( R - load ) < r && r < ( R + load ) ) Noverlap++;
      }
 
      if ( !Noverlap ) continue;
      float tmpRatio = float( Noverlap ) / float( Nall );
 
      if ( tmpRatio > OverlapRatioMax )
      {
        OverlapRatioMax = tmpRatio;
        MaxID           = i1;
      }
    }
    // jialk caution original is 0.8
    if ( OverlapRatioMax < 0.7 ) continue;
 
    //...Mask should be done...
    unsigned      MaskID[2] = { MaxID, i0 };
    AList<TMLink> l[2];
 
    for ( unsigned j0 = 0; j0 < 2; j0++ )
    {
      for ( unsigned j1 = 0; j1 < _tracks[MaskID[j0]]->nLinks(); j1++ )
      {
        TMLink& k = *_tracks[MaskID[j0]]->links()[j1];
        l[j0].append( k );
      }
    }
    // _tracks[i0]->links().append( _tracks[MaxID]->links() );
    // _tracks[MaxID]->links().append( _tracks[i0]->links ());
    _tracks[i0]->append( _tracks[MaxID]->links() );
    _tracks[MaxID]->append( _tracks[i0]->links() );
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    mask & merge " << std::endl;
    std::cout << "    0:";
    Dump( l[0], "flag sort" );
    std::cout << "    1:";
    Dump( l[1], "flag sort" );
    std::cout << std::endl;
#endif
 
    //...Which should be masked out ?...
    unsigned maskSide = 2;
 
#if 0
    //...0. Check by # of super layers... ( not applied now )
    unsigned super0 = NSuperLayers(l[0]);
    unsigned super1 = NSuperLayers(l[1]);
 
    if( super0 < super1 ) maskSide = 0;
    else if ( super0 > super1 ) maskSide = 1;
 
#  ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    NSuperLayers 0, 1 = " << NSuperLayers(l[0]) << ", ";
    std::cout << NSuperLayers(l[1]) << std::endl;
#  endif
 
    if (maskSide == 2) {
#endif
 
    //...1. Check by the inner-most layer...
    unsigned inner0 = InnerMost( l[0] )->wire()->layerId();
    unsigned inner1 = InnerMost( l[1] )->wire()->layerId();
    if ( inner0 < inner1 ) maskSide = 1;
    else if ( inner0 > inner1 ) maskSide = 0;
 
    if ( maskSide == 2 )
    {
 
      //...2. Check by dz
 
      //...Make two tracks...
      TTrack* tt[2];
      tt[0] = new TTrack( *( _tracks[MaskID[0]] ) );
      tt[1] = new TTrack( *( _tracks[MaskID[1]] ) );
      _fitter.fit( *tt[0] );
      _fitter.fit( *tt[1] );
      Helix h0 = Helix( tt[0]->helix() );
      Helix h1 = Helix( tt[1]->helix() );
 
      //...Check dz...
      h0.pivot( ORIGIN );
      h1.pivot( ORIGIN );
      if ( fabs( h0.dz() ) < fabs( h1.dz() ) ) maskSide = 1;
      else maskSide = 0;
 
      delete tt[0];
      delete tt[1];
    }
#if 0
    }
#endif
    bads.append( _tracks[MaskID[maskSide]] );
    flagTrk[MaskID[maskSide]] = 1;
  }
 
  // cout<<"bads: "<<bads.length()<<endl;
  _tracks.remove( bads );
 
  //*****  Masking *****
  n = _tracks.length();
 
  for ( unsigned i = 0; i < n; i++ )
  {
    TTrack& t = *_tracks[i];
    for ( unsigned j = 0; j < t.links().length(); j++ )
    {
      TMLink&            l    = *t.links()[j];
      const TMDCWireHit& whit = *l.hit();
 
      if ( !( whit.state() & WireHitFittingValid ) ) continue;
 
      // within half circle or not?
      double q = t.helix().center().x() * l.positionOnTrack().y() -
                 t.helix().center().y() * l.positionOnTrack().x();
      double qq = q * t.charge();
 
      if ( qq > 0 ) whit.state( whit.state() & ~WireHitInvalidForFit );
      else whit.state( whit.state() | WireHitInvalidForFit );
    }
  }
 
  free( flagTrk );
}
// end of addition
 
int TTrackManager::copyTrack( TTrack& t, MdcRec_trk** pr, MdcRec_trk_add** pra ) const {
 
  static const unsigned GoodHitMask =
      ( WireHitTimeValid | WireHitChargeValid | WireHitFindingValid | WireHitFittingValid );
  int err = 0;
 
  //...Hit loop...
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "    checking hits ... " << t.name() << " quality = " << t.quality();
  std::cout << " : " << t.cores().length() << ", " << t.ndf() << " : ";
  unsigned nnn = 0;
#endif
  unsigned      j         = 0;
  unsigned      nClst     = 0;
  unsigned      nStereos  = 0;
  unsigned      nOccupied = 0;
  AList<TMLink> hits;
  AList<TMLink> badHits;
  unsigned      n = t.links().length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TMLink*        l = t.links()[i];
    MdcRec_wirhit* h = l->hit()->reccdc();
 
#ifdef TRKRECO_DEBUG_DETAIL
    if ( h->trk ) std::cout << l->wire()->name() << "(n/a),";
#endif
 
    if ( h->trk )
    {
      ++nOccupied;
      if ( !( h->stat & WireHitInvalidForFit ) ) continue;
    }
    if ( ( l->hit()->state() & GoodHitMask ) == GoodHitMask )
    {
      if ( l->hit()->state() & WireHitInvalidForFit )
      {
        if ( !( h->stat & WireHitInvalidForFit ) ) badHits.append( l );
      }
      else
      {
        hits.append( l );
        if ( l->wire()->stereo() ) ++nStereos;
      }
    }
  }
  t.finalHits( hits );
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << std::endl;
#endif
 
  //...Check # of hits...
  if ( t.quality() & TrackQuality2D )
  {
    if ( hits.length() < 3 ) err = 3;
    if ( nOccupied > 2 ) err = 4;
  }
  else
  {
    if ( hits.length() < 5 ) err = 1;
    if ( nStereos < 2 ) err = 2;
  }
  if ( err ) return err;
 
  //...Create new tables...
  //    * pr = (reccdc_trk *) BsNewEnt(RECMDC_TRK);
  //    * pra = (reccdc_trk_add *) BsNewEnt(RECMDC_TRK_ADD);
  //    reccdc_trk * r = * pr;
  //    reccdc_trk_add * ra = * pra;
  *pr                = new MdcRec_trk;
  *pra               = new MdcRec_trk_add;
  MdcRec_trk*     r  = *pr;
  MdcRec_trk_add* ra = *pra;
 
  //...Copy hit information...
  // const AList<TMLink> & cores = t.cores();
  // const AList<TMLink> & links = t.links();
  // unsigned allHits = cores.length();
  // unsigned stereoHits = NStereoHits(cores);
  // r.m_chiSq = t.chi2();
  // r.m_confl = t.confidenceLevel();
  // r.m_ndf = t.ndf();
  // r.m_nhits = allHits;
  // r.m_nster = stereoHits;
  float    chisq = 0.;
  unsigned nHits = hits.length();
  for ( unsigned i = 0; i < nHits; i++ )
  {
    TMLink*        l = hits[i];
    MdcRec_wirhit* h = hits[i]->hit()->reccdc();
    h->trk           = r;
    h->pChiSq        = l->pull();
    h->lr            = l->leftRight();
    // zsl  if (l->usecathode() == 4) ++nClst;
    chisq += h->pChiSq;
  }
  r->chiSq = chisq;
  r->nhits = nHits;
  r->nster = nStereos;
  r->ndf   = nHits - 5;
  if ( t.quality() & TrackQuality2D ) r->ndf = nHits - 3;
 
  //...Bad hits...
  n = badHits.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    MdcRec_wirhit* h = badHits[i]->hit()->reccdc();
    h->trk           = r;
    h->stat |= WireHitInvalidForFit;
  }
 
  //...Cathode...
  r->nclus = nClst;
 
  //...Helix parameter...
  const Vector&     a  = t.helix().a();
  const SymMatrix&  ea = t.helix().Ea();
  const HepPoint3D& x  = t.helix().pivot();
  r->helix[0]          = a[0];
  r->helix[1]          = a[1];
  r->helix[2]          = a[2];
  r->helix[3]          = a[3];
  r->helix[4]          = a[4];
 
  r->pivot[0] = x.x();
  r->pivot[1] = x.y();
  r->pivot[2] = x.z();
 
  r->error[0]  = ea[0][0];
  r->error[1]  = ea[1][0];
  r->error[2]  = ea[1][1];
  r->error[3]  = ea[2][0];
  r->error[4]  = ea[2][1];
  r->error[5]  = ea[2][2];
  r->error[6]  = ea[3][0];
  r->error[7]  = ea[3][1];
  r->error[8]  = ea[3][2];
  r->error[9]  = ea[3][3];
  r->error[10] = ea[4][0];
  r->error[11] = ea[4][1];
  r->error[12] = ea[4][2];
  r->error[13] = ea[4][3];
  r->error[14] = ea[4][4];
 
  //...Get outer most hit(=termination point)...
  TMLink* last = OuterMost( hits );
 
  //...Calculate phi of the termination point...
  t.approach( *last );
  r->fiTerm = last->dPhi();
 
  return err;
}
 
void TTrackManager::sortTracksByQuality( void ) {
  unsigned n = _tracks.length();
  if ( n < 2 ) return;
 
  for ( unsigned i = 0; i < n - 1; i++ )
  {
    TTrack& t0         = *_tracks[i];
    float   bestRChisq = HUGE_VAL;
    if ( t0.ndf() > 0 ) bestRChisq = t0.chi2() / t0.ndf();
    for ( unsigned j = i + 1; j < n; j++ )
    {
      TTrack& t1     = *_tracks[j];
      float   rChisq = HUGE_VAL;
      if ( t1.ndf() > 0 ) rChisq = t1.chi2() / t1.ndf();
      if ( rChisq < bestRChisq )
      {
        bestRChisq = rChisq;
        _tracks.swap( i, j );
      }
    }
  }
}
 
void TTrackManager::sortTracksByPt( void ) {
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "trkmgr::sortTracksByPt : # of tracks=" << _tracks.length() << std::endl;
#endif
 
  unsigned n = _tracks.length();
  if ( n < 2 ) return;
 
  for ( unsigned i = 0; i < n - 1; i++ )
  {
    TTrack& t0     = *_tracks[i];
    float   bestPt = t0.pt();
    for ( unsigned j = i + 1; j < n; j++ )
    {
      TTrack& t1 = *_tracks[j];
      float   pt = t1.pt();
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout << "i,j=" << i << "," << j << " : pt i,j=" << bestPt << "," << pt << std::endl;
#endif
      if ( pt > bestPt )
      {
        bestPt = pt;
        _tracks.swap( i, j );
      }
    }
  }
}
 
void TTrackManager::treatCurler( MdcTrk& trk1, MdcRec_trk_add& cdc1, unsigned flag ) const {
  //...Originally coded by j.tanaka...
 
  //...Check inputs...
  if ( trk1.zero[2] == 0 ) return;
  if ( cdc1.daughter == 0 ) return;
 
  //...The other side...
  //    reccdc_trk_add & cdc2 = * (reccdc_trk_add *) BsGetEnt(RECMDC_TRK_ADD,
  //                              cdc1.m_daughter,
  //                              BBS_No_Index);
  //    MdcRec_trk_add & cdc2 = (*MdcRecTrkAddCol::getMdcRecTrkAddCol())[cdc1.daughter.id];
  MdcRec_trk_add& cdc2 = *cdc1.daughter->add;
 
  if ( cdc2.daughter == 0 ) return;
  if ( cdc2.rectrk == 0 ) return;
  //    rectrk & trk2 = * (rectrk *) BsGetEnt(RECTRK, cdc2.m_rectrk, BBS_No_Index);
  MdcTrk& trk2 = *cdc2.rectrk;
 
  if ( trk2.zero[2] == 0 ) return;
 
  //...Obtain RECTRK_LOCALZ...
  //    rectrk_localz & z1 = * (rectrk_localz *) BsGetEnt(RECTRK_LOCALZ,
  //                              trk1.m_zero[2],
  //                              BBS_No_Index);
  //    rectrk_localz & z2 = * (rectrk_localz *) BsGetEnt(RECTRK_LOCALZ,
  //                              trk2.m_zero[2],
  //                              BBS_No_Index);
  MdcTrk_localz& z1 = *trk1.zero[2];
  MdcTrk_localz& z2 = *trk2.zero[2];
 
  //...Pointer to mother and daughter...
  MdcRec_trk_add* mother   = &cdc1;
  MdcRec_trk_add* daughter = &cdc2;
 
  //      //...By dr...
  //      if (flag == 1) {
  //    float dr1 = fabs(z1.m_helix[0]);
  //    float dr2 = fabs(z2.m_helix[0]);
  //    if (dr1 > dr2) {
  //        mother = & cdc2;
  //        daughter = & cdc1;
  //    }
  //      }
 
  //      //...By dz...
  //      else {
  //    float dz1 = fabs(z1.m_helix[3]);
  //    float dz2 = fabs(z2.m_helix[3]);
  //    if (dz1 > dz2) {
  //        mother = & cdc2;
  //        daughter = & cdc1;
  //    }
  //      }
 
  //...By dz + dr...
  if ( flag == 3 )
  {
    float dz1 = fabs( z1.helix[3] );
    float dz2 = fabs( z2.helix[3] );
    if ( fabs( dz1 - dz2 ) < 2. ) flag = 1;
    else flag = 2;
  }
 
  //...By dr...
  if ( flag == 1 )
  {
    float dr1 = fabs( z1.helix[0] );
    float dr2 = fabs( z2.helix[0] );
    if ( dr1 > dr2 )
    {
      mother   = &cdc2;
      daughter = &cdc1;
    }
  }
 
  //...By dz...
  else if ( flag == 2 )
  {
    float dz1 = fabs( z1.helix[3] );
    float dz2 = fabs( z2.helix[3] );
    if ( dz1 > dz2 )
    {
      mother   = &cdc2;
      daughter = &cdc1;
    }
  }
 
  //...Update information...
  mother->quality &= ( ~TrackQualityOutsideCurler );
  mother->likelihood[0] = 1.;
  mother->decision |= TrackTrackManager;
  // zsl
  mother->daughter = daughter->body;
  mother->mother   = 0;
  // zsl end;
  daughter->quality |= TrackQualityOutsideCurler;
  daughter->likelihood[0] = 0.;
  daughter->mother        = mother->body;
  daughter->daughter      = 0;
  daughter->decision |= TrackTrackManager;
}
 
void TTrackManager::sortBanksByPt( void ) const {
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "trkmgr::sortBanksByPt : # of tracks="
            //   << BsCouTab(RECMDC_TRK_ADD) << std::endl;
            << MdcRecTrkAddCol::getMdcRecTrkAddCol()->size() << std::endl;
#endif
 
  //    unsigned n = BsCouTab(RECMDC_TRK_ADD);
  unsigned n = MdcRecTrkAddCol::getMdcRecTrkAddCol()->size();
  if ( n < 2 ) return;
 
  //...Sort RECMDC...
  unsigned* id = (unsigned*)malloc( n * sizeof( unsigned ) );
  for ( unsigned i = 0; i < n; i++ ) id[i] = i;
  for ( unsigned i = 0; i < n - 1; i++ )
  {
    //  reccdc_trk & cdc0 =
    //      * (reccdc_trk *) BsGetEnt(RECMDC_TRK,
    //                    i + 1,
    //                    BBS_No_Index);
    //  reccdc_trk_add & add0 =
    //      * (reccdc_trk_add *) BsGetEnt(RECMDC_TRK_ADD,
    //                    i + 1,
    //                    BBS_No_Index);
    //  reccdc_mctrk & mc0 =
    //      * (reccdc_mctrk *) BsGetEnt(RECMDC_MCTRK,
    //                  i + 1,
    //                  BBS_No_Index);
    MdcRec_trk&     cdc0 = ( *MdcRecTrkCol::getMdcRecTrkCol() )[i];
    MdcRec_trk_add& add0 = ( *MdcRecTrkAddCol::getMdcRecTrkAddCol() )[i];
    MdcRec_mctrk&   mc0  = ( *MdcRecMctrkCol::getMdcRecMctrkCol() )[i];
 
    float    bestPt      = 1. / fabs( cdc0.helix[2] );
    unsigned bestQuality = add0.quality;
    for ( unsigned j = i + 1; j < n; j++ )
    {
      //        reccdc_trk & cdc1 =
      //        * (reccdc_trk *) BsGetEnt(RECMDC_TRK,
      //                      j + 1,
      //                      BBS_No_Index);
      //        reccdc_trk_add & add1 =
      //        * (reccdc_trk_add *) BsGetEnt(RECMDC_TRK_ADD,
      //                          j + 1,
      //                          BBS_No_Index);
      //        reccdc_mctrk & mc1 =
      //        * (reccdc_mctrk *) BsGetEnt(RECMDC_MCTRK,
      //                        j + 1,
      //                        BBS_No_Index);
      MdcRec_trk&     cdc1 = ( *MdcRecTrkCol::getMdcRecTrkCol() )[j];
      MdcRec_trk_add& add1 = ( *MdcRecTrkAddCol::getMdcRecTrkAddCol() )[j];
      MdcRec_mctrk&   mc1  = ( *MdcRecMctrkCol::getMdcRecMctrkCol() )[j];
 
      float pt = 1. / fabs( cdc1.helix[2] );
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout << "i,j=" << i << "," << j << " : quality i,j=" << bestQuality << ","
                << add1.quality << std::endl;
#endif
      unsigned quality = add1.quality;
      if ( quality > bestQuality ) continue;
      else if ( quality < bestQuality )
      {
        bestQuality = quality;
        bestPt      = pt;
        swapReccdc( cdc0, add0, mc0, cdc1, add1, mc1 );
        unsigned tmp = id[i];
        id[i]        = id[j];
        id[j]        = tmp;
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "swapped" << std::endl;
#endif
        continue;
      }
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout << "i,j=" << i << "," << j << " : pt i,j=" << bestPt << "," << pt << std::endl;
#endif
      if ( pt > bestPt )
      {
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "swapping ... " << &cdc0 << "," << &add0 << "," << &mc0 << " <-> "
                  << &cdc1 << "," << &add1 << "," << &mc1 << std::endl;
#endif
        bestQuality = quality;
        bestPt      = pt;
        swapReccdc( cdc0, add0, mc0, cdc1, add1, mc1 );
        unsigned tmp = id[i];
        id[i]        = id[j];
        id[j]        = tmp;
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << "swapped" << std::endl;
#endif
      }
    }
  }
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "trkmgr::sortBanksByPt : first phase finished" << std::endl;
#endif
 
  tagReccdc( id, n );
  free( id );
 
#ifdef TRKRECO_DEBUG_DETAIL
  std::cout << "trkmgr::sortBanksByPt : second phase finished" << std::endl;
#endif
 
#if 0
  //...Sort RECTRK...
  n = BsCouTab(RECTRK);
  id = (unsigned *) malloc(n * sizeof(unsigned));
  for (unsigned i = 0; i < n; i++) id[i] = i;
  if (n > 1) {
    unsigned i = 0;
    while (i < n - 1) {
      rectrk & t = * (rectrk *) BsGetEnt(RECTRK, i + 1, BBS_No_Index);
      if (t.m_prekal == (i + 1)) {
    ++i;
    continue;
      }
 
      rectrk & s = * (rectrk *) BsGetEnt(RECTRK,
      t.m_prekal,
      BBS_No_Index);
 
      swapRectrk(t, s);
      unsigned tmp = id[i];
      id[i] = id[s.m_ID - 1];
      id[s.m_ID - 1] = tmp;
 
      //std::cout << "swap " << i + 1 << " and " << s.m_ID << std::endl;
 
      reccdc_trk_add & a =
    * (reccdc_trk_add *) BsGetEnt(RECMDC_TRK_ADD,
        t.m_prekal,
        BBS_No_Index);
      reccdc_trk_add & b =
    * (reccdc_trk_add *) BsGetEnt(RECMDC_TRK_ADD,
        s.m_prekal,
        BBS_No_Index);
      a.m_rectrk = t.m_ID;
      b.m_rectrk = s.m_ID;
    }
  }
#else
  /*
  // jtanaka 000925 -->
  n = BsCouTab(RECTRK);
  id = (unsigned *) malloc(n * sizeof(unsigned));
  for (unsigned i = 0; i < n; i++) id[i] = i;
  int foundId = 0;
  while(foundId != n){
  rectrk & t = * (rectrk *) BsGetEnt(RECTRK, foundId + 1, BBS_No_Index);
  int minPrekal = t.m_prekal;
  int exchangeId = foundId;
  for(int i=foundId+1;i<n;++i){
  rectrk & s = * (rectrk *) BsGetEnt(RECTRK, i + 1, BBS_No_Index);
  if(s.m_prekal < minPrekal){
  minPrekal = s.m_prekal;
  exchangeId = i;
  }
  }
  if(exchangeId != foundId){
  rectrk & s = * (rectrk *) BsGetEnt(RECTRK,
  exchangeId + 1,
  BBS_No_Index);
 
  swapRectrk(t, s);
  unsigned tmp = id[t.m_ID - 1];
  id[t.m_ID - 1] = id[s.m_ID - 1];
  id[s.m_ID - 1] = tmp;
 
  reccdc_trk_add & a =
   * (reccdc_trk_add *) BsGetEnt(RECMDC_TRK_ADD,
   t.m_prekal,
   BBS_No_Index);
   reccdc_trk_add & b =
   * (reccdc_trk_add *) BsGetEnt(RECMDC_TRK_ADD,
   s.m_prekal,
   BBS_No_Index);
   a.m_rectrk = t.m_ID;
   b.m_rectrk = s.m_ID;
   }
   ++foundId;
   }
  // <-- jtanaka 000925
  */
#endif
 
  tagRectrk( id, n );
  free( id );
}
 
void TTrackManager::swapReccdc( MdcRec_trk& cdc0, MdcRec_trk_add& add0, MdcRec_mctrk& mc0,
                                MdcRec_trk& cdc1, MdcRec_trk_add& add1,
                                MdcRec_mctrk& mc1 ) const {
#define RECMDC_ACTUAL_SIZE 124
#define RECMDCADD_ACTUAL_SIZE 40
#define RECMDCMC_ACTUAL_SIZE 28
 
  static bool  first = true;
  static void* swapRegion;
  if ( first )
  {
    first      = false;
    swapRegion = malloc( RECMDC_ACTUAL_SIZE );
  }
 
  void* s0 = &cdc0.helix[0];
  void* s1 = &cdc1.helix[0];
  memcpy( swapRegion, s0, RECMDC_ACTUAL_SIZE );
  memcpy( s0, s1, RECMDC_ACTUAL_SIZE );
  memcpy( s1, swapRegion, RECMDC_ACTUAL_SIZE );
 
  s0 = &add0.quality;
  s1 = &add1.quality;
  memcpy( swapRegion, s0, RECMDCADD_ACTUAL_SIZE );
  memcpy( s0, s1, RECMDCADD_ACTUAL_SIZE );
  memcpy( s1, swapRegion, RECMDCADD_ACTUAL_SIZE );
 
  if ( ( &mc0 ) && ( &mc1 ) )
  {
    s0 = &mc0.hep;
    s1 = &mc1.hep;
    memcpy( swapRegion, s0, RECMDCMC_ACTUAL_SIZE );
    memcpy( s0, s1, RECMDCMC_ACTUAL_SIZE );
    memcpy( s1, swapRegion, RECMDCMC_ACTUAL_SIZE );
  }
}
 
void TTrackManager::swapRectrk( MdcTrk& rec0, MdcTrk& rec1 ) const {
#define RECTRK_ACTUAL_SIZE 84
 
  static bool  first = true;
  static void* swapRegion;
  if ( first )
  {
    first      = false;
    swapRegion = malloc( RECTRK_ACTUAL_SIZE );
  }
 
  void* s0 = &rec0.glob[0];
  void* s1 = &rec1.glob[0];
  memcpy( swapRegion, s0, RECTRK_ACTUAL_SIZE );
  memcpy( s0, s1, RECTRK_ACTUAL_SIZE );
  memcpy( s1, swapRegion, RECTRK_ACTUAL_SIZE );
}
 
void TTrackManager::tagReccdc( unsigned* id0, unsigned nTrk ) const {
  /*
     unsigned * id = (unsigned *) malloc(nTrk * sizeof(unsigned));
     for (unsigned i = 0; i < nTrk; i++)
     id[id0[i]] = i;
 
#ifdef TRKRECO_DEBUG_DETAIL
for (unsigned i = 0; i < nTrk; i++)
std::cout << "id0 " << i << " ... " << id0[i] << std::endl;
for (unsigned i = 0; i < nTrk; i++)
std::cout << "id  " << i << " ... " << id[i] << std::endl;
#endif
  //    unsigned n = BsCouTab(RECMDC_TRK_ADD);
  //    unsigned n = MdcRecTrkAddCol::getMdcRecTrkAddCol()->size();
 
  //    for (unsigned i = 0; i < n; i++) {
  //    reccdc_trk_add & w = * (reccdc_trk_add *) BsGetEnt(RECMDC_TRK_ADD,
  //                               i + 1,
  //                               BBS_No_Index);
  //    MdcRec_trk_add & w = (* MdcRecTrkAddCol::getMdcRecTrkAddCol())[i];
  //    if (w.mother) w.mother = id[w.mother - 1] + 1;
  //    if (w.daughter) w.daughter = id[w.daughter - 1] + 1;
  //    }
 
#ifdef TRKRECO_DEBUG_DETAIL
std::cout << "TTrackManager::tagReccdc ... RECMDC_TRK_ADD done" << std::endl;
#endif
 
n = BsCouTab(RECMDC_WIRHIT);
for (unsigned i = 0; i < n; i++) {
reccdc_wirhit & w = * (reccdc_wirhit *) BsGetEnt(RECMDC_WIRHIT,
i + 1,
BBS_No_Index);
if (w.m_trk == 0) continue;
w.m_trk = id[w.m_trk - 1] + 1;
}
 
#ifdef TRKRECO_DEBUG_DETAIL
std::cout << "TTrackManager::tagReccdc ... RECMDC_WIRHIT done" << std::endl;
#endif
 
n = BsCouTab(DATMDC_MCWIRHIT);
for (unsigned i = 0; i < n; i++) {
datcdc_mcwirhit & m =
   * (datcdc_mcwirhit *) BsGetEnt(DATMDC_MCWIRHIT,i + 1,BBS_No_Index);
   if (m.m_trk == 0) continue;
   m.m_trk = id[m.m_trk - 1] + 1;
   }
 
#ifdef TRKRECO_DEBUG_DETAIL
std::cout << "TTrackManager::tagReccdc ... DATMDC_MCWIRHIT done" << std::endl;
#endif
 
n = BsCouTab(RECTRK);
for (unsigned i = 0; i < n; i++) {
rectrk & r = * (rectrk *) BsGetEnt(RECTRK, i + 1, BBS_No_Index);
if (r.m_prekal == 0) continue;
r.m_prekal = id[r.m_prekal - 1] + 1;
}
 
#ifdef TRKRECO_DEBUG_DETAIL
std::cout << "TTrackManager::tagReccdc ... RECTRK done" << std::endl;
#endif
 
// jtanaka
n = BsCouTab(RECMDC_SVD_TRK);
for (unsigned i = 0; i < n; i++) {
reccdc_svd_trk & r = * (reccdc_svd_trk *) BsGetEnt(RECMDC_SVD_TRK, i + 1, BBS_No_Index);
if (r.m_cdc_trk == 0) continue;
r.m_cdc_trk = id[r.m_cdc_trk - 1] + 1;
}
 
#ifdef TRKRECO_DEBUG_DETAIL
std::cout << "TTrackManager::tagReccdc ... RECMDC_SVD_TRK done" << std::endl;
#endif
 
free(id);
*/
}
 
void TTrackManager::setCurlerFlags( void ) {
  unsigned n = _tracks.length();
  if ( n < 2 ) return;
 
  for ( unsigned i = 0; i < n - 1; i++ )
  {
    TTrack& t0 = *_tracks[i];
    if ( t0.type() != TrackTypeCurl ) continue;
    float c0 = t0.charge();
 
    for ( unsigned j = i + 1; j < n; j++ )
    {
      TTrack& t1 = *_tracks[j];
      float   c1 = t1.charge();
      if ( c0 * c1 > 0. ) continue;
      if ( t1.type() != TrackTypeCurl ) continue;
 
      bool     toBeMerged = false;
      unsigned n0         = t0.testByApproach( t1.cores(), _sigmaCurlerMergeTest );
      if ( n0 > _nCurlerMergeTest ) toBeMerged = true;
      if ( !toBeMerged )
      {
        unsigned n1 = t1.testByApproach( t0.cores(), _sigmaCurlerMergeTest );
        if ( n1 > _nCurlerMergeTest ) toBeMerged = true;
      }
 
      if ( toBeMerged )
      {
        //      ++_s->_nToBeMerged;
        if ( ( t0.daughter() ) || ( t1.daughter() ) ) ++_s->_nToBeMergedMoreThanTwo;
        t0.daughter( &t1 );
        t1.daughter( &t0 );
      }
    }
  }
}
 
void TTrackManager::salvageAssociateHits( const AList<TMDCWireHit>& hits, float maxSigma2 ) {
  // #define TRKRECO_DEBUG
  // #define TRKRECO_DEBUG_DETAIL
 
#ifdef TRKRECO_DEBUG
  std::cout << "... trkmgr::salvaging associate hits" << std::endl;
  std::cout << "    # of given hits=" << hits.length() << std::endl;
#endif
 
  //...Check arguments...
  unsigned nTracks = _tracks.length();
  if ( nTracks == 0 ) return;
  unsigned nHits = hits.length();
  if ( nHits == 0 ) return;
 
  static const TPoint2D o( 0., 0. );
 
  //...Hit loop...
  for ( unsigned i = 0; i < nHits; i++ )
  {
    TMDCWireHit& h = *hits[i];
 
    //...Already used?...
    if ( h.state() & WireHitUsed ) continue;
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << "    checking " << h.wire()->name();
#endif
 
    //...Track loop...
    AList<TMLink> toBeDeleted;
    TMLink*       best      = NULL;
    TTrack*       bestTrack = NULL;
    for ( unsigned j = 0; j < nTracks; j++ )
    {
      TTrack& t = *_tracks[j];
 
      //...Pre-selection...
      TPoint2D c  = t.center();
      TPoint2D co = -c;
      TPoint2D x  = h.wire()->xyPosition();
 
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout << " : c= " << co.cross( x - c ) * t.charge();
      std::cout << ", d=" << fabs( ( x - c ).mag() - fabs( t.radius() ) );
#endif
 
      if ( co.cross( x - c ) * t.charge() > 0. ) continue;
      if ( fabs( ( x - c ).mag() - fabs( t.radius() ) ) > 5. ) continue;
 
      //...Try to append this hit...
      TMLink& link = *new TMLink( 0, &h );
      int     err  = t.approach( link );
      if ( err < 0 )
      {
#ifdef TRKRECO_DEBUG_DETAIL
        std::cout << " : " << t.name() << " approach failure";
#endif
        toBeDeleted.append( link );
        continue;
      }
 
      //...Calculate sigma...
      float distance = link.distance();
      float diff     = fabs( distance - link.hit()->drift() );
      float sigma    = diff / link.hit()->dDrift();
      link.pull( sigma * sigma );
 
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout << " : " << t.name() << " pull = " << link.pull();
#endif
      if ( link.pull() > maxSigma2 )
      {
        toBeDeleted.append( link );
        continue;
      }
 
      if ( best )
      {
        if ( best->pull() > link.pull() )
        {
          toBeDeleted.append( best );
          best      = &link;
          bestTrack = &t;
        }
        else { toBeDeleted.append( link ); }
      }
      else
      {
        best      = &link;
        bestTrack = &t;
      }
    }
 
    if ( best )
    {
      bestTrack->append( *best );
      best->hit()->state( best->hit()->state() | WireHitInvalidForFit );
      _associateHits.append( best );
#ifdef TRKRECO_DEBUG
      std::cout << "    " << best->hit()->wire()->name();
      std::cout << " -> " << bestTrack->name() << std::endl;
#endif
    }
    HepAListDeleteAll( toBeDeleted );
 
#ifdef TRKRECO_DEBUG_DETAIL
    std::cout << std::endl;
#endif
  }
}
 
void TTrackManager::maskBadHits( const AList<TTrack>& tracks, float maxSigma2 ) {
#ifdef TRKRECO_DEBUG
  std::cout << "... trkmgr::maskBadHits" << std::endl;
#endif
 
  unsigned n = tracks.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    TTrack&             t           = *tracks[i];
    bool                toBeUpdated = false;
    const AList<TMLink> links       = t.links();
    unsigned            nHits       = links.length();
    for ( unsigned j = 0; j < nHits; j++ )
    {
      if ( links[j]->pull() > maxSigma2 )
      {
        links[j]->hit()->state( links[j]->hit()->state() | WireHitInvalidForFit );
        toBeUpdated = true;
#ifdef TRKRECO_DEBUG
        std::cout << "    " << t.name() << " : ";
        std::cout << links[j]->wire()->name() << "(pull=";
        std::cout << links[j]->pull() << ") is masked" << std::endl;
#endif
      }
    }
    if ( toBeUpdated ) t.update();
  }
}
 
void TTrackManager::clearTables( void ) const {
  //    BsDelEnt(RECMDC_TRK, BBS_ID_ALL);
  //    BsDelEnt(RECMDC_TRK_ADD, BBS_ID_ALL);
  //    BsDelEnt(RECMDC_MCTRK, BBS_ID_ALL);
  //    BsDelEnt(RECMDC_MCTRK2HEP, BBS_ID_ALL);
  unsigned n = MdcRecTrkCol::getMdcRecTrkCol()->size();
  for ( unsigned i = 0; i < n; i++ ) delete &( *MdcRecTrkCol::getMdcRecTrkCol() )[i];
 
  n = MdcRecTrkAddCol::getMdcRecTrkAddCol()->size();
  for ( unsigned i = 0; i < n; i++ ) delete &( *MdcRecTrkAddCol::getMdcRecTrkAddCol() )[i];
 
  n = MdcRecMctrkCol::getMdcRecMctrkCol()->size();
  for ( unsigned i = 0; i < n; i++ ) delete &( *MdcRecMctrkCol::getMdcRecMctrkCol() )[i];
 
  n = MdcRecMctrk2hepCol::getMdcRecMctrk2hepCol()->size();
  for ( unsigned i = 0; i < n; i++ )
    delete &( *MdcRecMctrk2hepCol::getMdcRecMctrk2hepCol() )[i];
 
  //...Clear track association...
  //    unsigned n = BsCouTab(RECMDC_WIRHIT);
  n = MdcRecWirhitCol::getMdcRecWirhitCol()->size();
  for ( unsigned i = 0; i < n; i++ )
  {
    //  reccdc_wirhit & h = * (reccdc_wirhit *)
    //      BsGetEnt(RECMDC_WIRHIT, i + 1, BBS_No_Index);
    //  h.m_trk = 0;
    MdcRec_wirhit& h = ( *MdcRecWirhitCol::getMdcRecWirhitCol() )[i];
    h.trk            = 0;
  }
  //    n = BsCouTab(DATMDC_MCWIRHIT);
  n = MdcDatMcwirhitCol::getMdcDatMcwirhitCol()->size();
  for ( unsigned i = 0; i < n; i++ )
  {
    //  datcdc_mcwirhit & h = * (datcdc_mcwirhit *)
    //      BsGetEnt(DATMDC_MCWIRHIT, i + 1, BBS_No_Index);
    //  h.m_trk = 0;
    MdcDat_mcwirhit& h = ( *MdcDatMcwirhitCol::getMdcDatMcwirhitCol() )[i];
    h.trk              = 0;
  }
}
 
AList<TTrack> TTrackManager::selectGoodTracks( const AList<TTrack>& list,
                                               bool                 track2D ) const {
  AList<TTrack> goodTracks;
  unsigned      n = list.length();
  for ( unsigned i = 0; i < n; i++ )
  {
    const TTrack& t = *list[i];
    if ( !goodTrack( t, track2D ) ) continue;
 
    //...Remove super momentum...
    if ( _maxMomentum > 0. )
    {
      if ( t.ptot() > _maxMomentum )
      {
        //      ++_s->_nSuperMoms;
        continue;
      }
    }
 
    goodTracks.append( (TTrack&)t );
  }
 
  if ( _debugLevel )
  {
    if ( list.length() != goodTracks.length() )
    {
      std::cout << "TTrackManager::selectGoodTracks ... bad tracks found" << std::endl
                << "                                    # of bad tracks = "
                << list.length() - goodTracks.length() << " : 2D flag = " << track2D
                << std::endl;
      AList<TTrack> tmp;
      tmp.append( list );
      tmp.remove( goodTracks );
      std::cout << "    Track dump" << std::endl;
      for ( unsigned i = 0; i < tmp.length(); i++ )
      { std::cout << "    " << TrackDump( *tmp[i] ) << std::endl; }
    }
  }
 
  return goodTracks;
}
 
bool TTrackManager::checkNumberOfHits( const TTrack& t, bool track2D ) {
  const AList<TMLink>& cores = t.cores();
 
  if ( track2D )
  {
    unsigned axialHits = NAxialHits( cores );
    if ( axialHits < 3 ) return false;
  }
  else
  {
    unsigned allHits = cores.length();
    if ( allHits < 5 ) return false;
    unsigned stereoHits = NStereoHits( cores );
    if ( stereoHits < 2 ) return false;
    unsigned axialHits = allHits - stereoHits;
    if ( axialHits < 3 ) return false;
  }
  return true;
}
 
void TTrackManager::determineIP( void ) {
  static const HepVector3D InitialVertex( 0., 0., 0. );
 
  //...Track selection...
  //    unsigned n = BsCouTab(RECTRK);
  unsigned             n = MdcTrkCol::getMdcTrkCol()->size();
  AList<MdcTrk_localz> zList;
  for ( unsigned i = 0; i < n; i++ )
  {
    //  const rectrk & t = * (rectrk *) BsGetEnt(RECTRK, i + 1, BBS_No_Index);
    const MdcTrk& t = ( *MdcTrkCol::getMdcTrkCol() )[i];
 
    if ( t.prekal == 0 ) continue;
    //  const reccdc_trk_add & c = * (reccdc_trk_add *)
    //      BsGetEnt(RECMDC_TRK_ADD, t.m_prekal, BBS_No_Index);
    const MdcRec_trk_add& c = *t.prekal->add;
 
    //...Only good tracks...
    if ( c.quality ) continue;
 
    //...Require SVD hits...
    //  const rectrk_global & g = * (rectrk_global *) BsGetEnt(RECTRK_GLOBAL,
    //                                 t.m_glob[2],
    //                                 BBS_No_Index);
    const MdcTrk_global& g = *t.glob[2];
 
    if ( !&g ) continue;
    if ( g.nhits[3] < 2 ) continue;
    if ( g.nhits[4] < 2 ) continue;
 
    //...OK...
    //  const rectrk_localz & z = * (rectrk_localz *) BsGetEnt(RECTRK_LOCALZ,
    //                                 t.m_zero[2],
    //                                 BBS_No_Index);
    MdcTrk_localz& z = *t.zero[2];
 
    if ( !&z ) continue;
    //  zList.append((rectrk_localz &) z);
    zList.append( z );
  }
  unsigned nZ = zList.length();
  if ( nZ < 2 ) return;
 
  //...Fitting...
  //      kvertexfitter kvf;
  //      kvf.initialVertex(initialVertex);
  //      for (unsigned i = 0; i < nZ; i++) {
  //    kvf.addTrack();
  //      }
}
 
void TTrackManager::tagRectrk( unsigned* id0, unsigned nTrk ) const {
  /*
     unsigned * id = (unsigned *) malloc(nTrk * sizeof(unsigned));
     for (unsigned i = 0; i < nTrk; i++)
     id[id0[i]] = i;
 
  //      for (unsigned i = 0; i < nTrk; i++)
  //        std::cout << "id0 " << i << " ... " << id0[i] << std::endl;
  //      for (unsigned i = 0; i < nTrk; i++)
  //        std::cout << "id  " << i << " ... " << id[i] << std::endl;
  //      BsShwDat(RECTRK_TOF);
 
  unsigned n = BsCouTab(RECTRK_TOF);
  for (unsigned i = 0; i < n; i++) {
  rectrk_tof & t = * (rectrk_tof *) BsGetEnt(RECTRK_TOF,
  i + 1,
  BBS_No_Index);
  if (t.m_rectrk) t.m_rectrk = id[t.m_rectrk - 1] + 1;
  }
 
  //      BsShwDat(RECTRK_TOF);
 
  // jtanaka
  n = BsCouTab(RECSVD_HIT);
  for (unsigned i = 0; i < n; i++) {
  recsvd_hit & t = * (recsvd_hit *) BsGetEnt(RECSVD_HIT,
  i + 1,
  BBS_No_Index);
  if (t.m_trk) t.m_trk = id[t.m_trk - 1] + 1;
  }
 
  free(id);
  */
}
 
/*
// jtanaka 000925 -->
#define TRKRECO_REPLACE_TABLE 1
#if !(TRKRECO_REPLACE_TABLE)
void
copyRecMDC_trk_Table(const Reccdc_trk & org,
Reccdc_trk & copied)
{
copied.helix(0, org.helix(0));
copied.helix(1, org.helix(1));
copied.helix(2, org.helix(2));
copied.helix(3, org.helix(3));
copied.helix(4, org.helix(4));
copied.pivot(0, org.pivot(0));
copied.pivot(1, org.pivot(1));
copied.pivot(2, org.pivot(2));
copied.error(0, org.error(0));
copied.error(1, org.error(1));
copied.error(2, org.error(2));
copied.error(3, org.error(3));
copied.error(4, org.error(4));
copied.error(5, org.error(5));
copied.error(6, org.error(6));
copied.error(7, org.error(7));
copied.error(8, org.error(8));
copied.error(9, org.error(9));
copied.error(10, org.error(10));
copied.error(11, org.error(11));
copied.error(12, org.error(12));
copied.error(13, org.error(13));
copied.error(14, org.error(14));
copied.chiSq(org.chiSq());
copied.ndf(org.ndf());
copied.fiTerm(org.fiTerm());
copied.nhits(org.nhits());
copied.nster(org.nster());
copied.nclus(org.nclus());
copied.stat(org.stat());
copied.mass(org.mass());
}
 
void
copyRecMDC_trk_add_Table(const Reccdc_trk_add & org,
Reccdc_trk_add & copied)
{
copied.quality(org.quality());
copied.kind(org.kind());
copied.mother(org.mother());
copied.daughter(org.daughter());
copied.decision(org.decision());
copied.likelihood(0, org.likelihood(0));
copied.likelihood(1, org.likelihood(1));
copied.likelihood(2, org.likelihood(2));
copied.stat(org.stat());
copied.rectrk(org.rectrk());
}
 
void
copyRecMDC_MCtrk_Table(const Reccdc_mctrk & org,
Reccdc_mctrk & copied)
{
copied.hep(org.hep());
copied.wirFrac(org.wirFrac());
copied.wirFracHep(org.wirFracHep());
copied.charge(org.charge());
copied.ptFrac(org.ptFrac());
copied.pzFrac(org.pzFrac());
copied.quality(org.quality());
}
 
void
copyRecMDC_MCtrk2hep_Table(const Reccdc_mctrk2hep & org,
    Reccdc_mctrk2hep & copied)
{
  copied.wir(org.wir());
  copied.clust(org.clust());
  copied.trk(org.trk());
  copied.hep(org.hep());
}
#endif
 
void
TTrackManager::addSvd(const int mcFlag) const {
  TSvdAssociator svdA(-20000.,20000.);
  svdA.fillClusters();
 
  //BsShwDat(RECTRK);
  //BsShwDat(RECMDC_TRK);
  //BsShwDat(RECMDC_MCTRK);
 
  Reccdc_trk_Manager& trkMgr =
    Reccdc_trk_Manager::get_manager();
  Reccdc_trk_add_Manager& trkMgr2 =
    Reccdc_trk_add_Manager::get_manager();
  Reccdc_svd_trk_Manager& svdMgr =
    Reccdc_svd_trk_Manager::get_manager();
 
#if !(TRKRECO_REPLACE_TABLE)
  Reccdc_wirhit_Manager& wirMgr =
    Reccdc_wirhit_Manager::get_manager();
  Reccdc_mctrk_Manager& mcMgr =
    Reccdc_mctrk_Manager::get_manager();
  Reccdc_mctrk2hep_Manager& mcMgr2 =
    Reccdc_mctrk2hep_Manager::get_manager();
  Datcdc_mcwirhit_Manager& mcMgr3 =
    Datcdc_mcwirhit_Manager::get_manager();
#endif
 
  int nSize = trkMgr.count();
  for(int i=0;i<nSize;++i){
    //std::cout << "trk " << i << ": " << trkMgr[i].helix(0) << std::endl;
#if 1
    // Reconstruction Info --> SVD Recon.
    if(trkMgr2[i].decision() != TrackPMCurlFinder &&
        (trkMgr2[i].quality() & TrackQuality2D) != TrackQuality2D &&
        trkMgr[i].helix(2) != 0. && fabs(1./trkMgr[i].helix(2))<0.2){
      HepVector a(5);
      a[0] = trkMgr[i].helix(0);
      a[1] = trkMgr[i].helix(1);
      a[2] = trkMgr[i].helix(2);
      a[3] = trkMgr[i].helix(3);
      a[4] = trkMgr[i].helix(4);
      HepPoint3D p(trkMgr[i].pivot(0),
          trkMgr[i].pivot(1),
          trkMgr[i].pivot(2));
      Helix th(p,a);
      th.pivot(HepPoint3D(0.,0.,0.)); // pivot = (0,0,0)
      AList<TSvdHit> cand;
      double tz,tt;
      if(svdA.recTrk(th,tz,tt,0.5,50.0,cand,0.5)){
        //std::cout << "SVD in " << i << std::endl;
#if TRKRECO_REPLACE_TABLE
        Reccdc_svd_trk & newSvd  = svdMgr.add();
#else
        Reccdc_trk     & newTrk  = trkMgr.add();
        Reccdc_trk_add & newTrk2 = trkMgr2.add();
        Reccdc_svd_trk & newSvd  = svdMgr.add();
        // copy all information
        copyRecMDC_trk_Table(trkMgr[i],newTrk);
        copyRecMDC_trk_add_Table(trkMgr2[i],newTrk2);
#endif
        if(mcFlag){
#if TRKRECO_REPLACE_TABLE
          ;
#else
          Reccdc_mctrk & newMcTrk = mcMgr.add();
          copyRecMDC_MCtrk_Table(mcMgr[i],mcMgr[mcMgr.count()-1]);
          int nMCt2h = mcMgr2.count();
          for(int j=0;j<nMCt2h;++j){
            if(mcMgr2[j].trk() &&
                mcMgr2[j].trk().get_ID() == trkMgr[i].get_ID()){
              Reccdc_mctrk2hep & newMcTrk2Hep = mcMgr2.add();
              copyRecMDC_MCtrk2hep_Table(mcMgr2[j],newMcTrk2Hep);
              newMcTrk2Hep.trk(newTrk);
            }
          }
          int nMCwire = mcMgr3.count();
          for(int j=0;j<nMCwire;++j){
            if(mcMgr3[j].trk().get_ID() == trkMgr[i].get_ID()){
              mcMgr3[j].trk(newTrk);
            }
          }
#endif
        }
        HepVector ta = th.a(); // pivot = (0,0,0)
        ta[3] = tz;
        ta[4] = tt;
        th.a(ta);
        th.pivot(p); // pivot, (0,0,0) --> p
#if TRKRECO_REPLACE_TABLE
        trkMgr[i].helix(3, th.a()[3]);
        trkMgr[i].helix(4, th.a()[4]);
#else
        newTrk.helix(3, th.a()[3]);
        newTrk.helix(4, th.a()[4]);
#endif
 
        newSvd.Helix(0, ta[0]);
        newSvd.Helix(1, ta[1]);
        newSvd.Helix(2, ta[2]);
        newSvd.Helix(3, ta[3]);
        newSvd.Helix(4, ta[4]);
#if TRKRECO_REPLACE_TABLE
        newSvd.cdc_trk(trkMgr2[i]);
#else
        newSvd.cdc_trk(newTrk2);
#endif
        newSvd.Status(0); // 0 is normal.
        int indexSvd  = 0;
        for(int j=0;j<cand.length();++j){
          if(indexSvd >= 16)break;
          if((cand[j])->rphi() && (cand[j])->z()){
            newSvd.svd_cluster(indexSvd, *(cand[j]->rphi()));
            ++indexSvd;
            newSvd.svd_cluster(indexSvd, *(cand[j]->z()));
            ++indexSvd;
          }else{
            std::cerr << "[TTrackManager of TrkReco] Why ? no associated SVDhit ?" <<
std::endl;
          }
        }
#if TRKRECO_REPLACE_TABLE
        trkMgr2[i].quality(0); // set to 0 --> GOOD!
        trkMgr2[i].decision((trkMgr2[i].decision() | TrackSVDAssociator));
#else
        newTrk2.quality(1); // temporary
        newTrk2.decision((newTrk2.decision() | TrackSVDAssociator));
#endif
#if !(TRKRECO_REPLACE_TABLE)
        // MDC Wire information
        for(int j=0;j<wirMgr.count();++j){
          if(wirMgr[j].trk() &&
              wirMgr[j].trk().get_ID() == trkMgr[i].get_ID()){
            wirMgr[j].trk(newTrk);
          }
        }
#endif
      }else if(fabs(th.a()[3]) > 30.){
        if(svdA.recTrk(th,tz,tt,0.5,-1.0,cand,0.5)){
          //std::cout << "SVD in " << i << std::endl;
#if TRKRECO_REPLACE_TABLE
          Reccdc_svd_trk & newSvd  = svdMgr.add();
#else
          Reccdc_trk     & newTrk  = trkMgr.add();
          Reccdc_trk_add & newTrk2 = trkMgr2.add();
          Reccdc_svd_trk & newSvd  = svdMgr.add();
          // copy all information
          copyRecMDC_trk_Table(trkMgr[i],newTrk);
          copyRecMDC_trk_add_Table(trkMgr2[i],newTrk2);
#endif
          if(mcFlag){
#if TRKRECO_REPLACE_TABLE
            ;
#else
            Reccdc_mctrk & newMcTrk = mcMgr.add();
            copyRecMDC_MCtrk_Table(mcMgr[i],mcMgr[mcMgr.count()-1]);
            int nMCt2h = mcMgr2.count();
            for(int j=0;j<nMCt2h;++j){
              if(mcMgr2[j].trk() &&
                  mcMgr2[j].trk().get_ID() == trkMgr[i].get_ID()){
                Reccdc_mctrk2hep & newMcTrk2Hep = mcMgr2.add();
                copyRecMDC_MCtrk2hep_Table(mcMgr2[j],newMcTrk2Hep);
                newMcTrk2Hep.trk(newTrk);
              }
            }
            int nMCwire = mcMgr3.count();
            for(int j=0;j<nMCwire;++j){
              if(mcMgr3[j].trk().get_ID() == trkMgr[i].get_ID()){
                mcMgr3[j].trk(newTrk);
              }
            }
#endif
          }
          HepVector ta = th.a(); // pivot = (0,0,0)
          ta[3] = tz;
          ta[4] = tt;
          th.a(ta);
          th.pivot(p); // pivot, (0,0,0) --> p
#if TRKRECO_REPLACE_TABLE
          trkMgr[i].helix(3, th.a()[3]);
          trkMgr[i].helix(4, th.a()[4]);
#else
          newTrk.helix(3, th.a()[3]);
          newTrk.helix(4, th.a()[4]);
#endif
 
          newSvd.Helix(0, ta[0]);
          newSvd.Helix(1, ta[1]);
          newSvd.Helix(2, ta[2]);
          newSvd.Helix(3, ta[3]);
          newSvd.Helix(4, ta[4]);
#if TRKRECO_REPLACE_TABLE
          newSvd.cdc_trk(trkMgr2[i]);
#else
          newSvd.cdc_trk(newTrk2);
#endif
          newSvd.Status(0); // 0 is normal.
          int indexSvd  = 0;
          for(int j=0;j<cand.length();++j){
            if(indexSvd >= 16)break;
            if((cand[j])->rphi() && (cand[j])->z()){
              newSvd.svd_cluster(indexSvd, *(cand[j]->rphi()));
              ++indexSvd;
              newSvd.svd_cluster(indexSvd, *(cand[j]->z()));
              ++indexSvd;
            }else{
              std::cerr << "[TTrackManager of TrkReco] Why ? no associated SVDhit ?" <<
std::endl;
            }
          }
#if TRKRECO_REPLACE_TABLE
          trkMgr2[i].quality(0); // set to 0 --> GOOD!
          trkMgr2[i].decision((trkMgr2[i].decision() | TrackSVDAssociator));
#else
          newTrk2.quality(1); // temporary
          newTrk2.decision((newTrk2.decision() | TrackSVDAssociator));
#endif
#if !(TRKRECO_REPLACE_TABLE)
          // MDC Wire information
          for(int j=0;j<wirMgr.count();++j){
            if(wirMgr[j].trk() &&
                wirMgr[j].trk().get_ID() == trkMgr[i].get_ID()){
              wirMgr[j].trk(newTrk);
            }
          }
#endif
        }
      }
    }
  }
#endif
}
// <-- jtanaka 000925
*/
 
bool TTrackManager::goodTrack( const TTrack& t, bool track2D ) {
 
  //...Check number of hits...
  if ( !checkNumberOfHits( t, track2D ) ) return false;
 
  //...Check helix parameter...
  if ( HelixHasNan( t.helix() ) ) return false;
 
  return true;
}
 
IBesMagFieldSvc* TTrackManager::getPmgnIMF() const {
  if ( nullptr == m_pmgnIMF )
  {
    StatusCode sc = Gaudi::svcLocator()->service( "MagneticFieldSvc", m_pmgnIMF );
    if ( sc.isFailure() )
    { std::cout << "Unable to open Magnetic field service" << std::endl; }
  }
  if ( nullptr == m_pmgnIMF ) std::cout << " Unable to get IBesMagFieldSvc " << std::endl;
  return m_pmgnIMF;
}
 
IRawDataProviderSvc* TTrackManager::getRawDataProviderSvc() const {
  if ( !m_rawDataProviderSvc )
  {
    StatusCode sc = Gaudi::svcLocator()->service( "RawDataProviderSvc", m_rawDataProviderSvc );
    if ( sc.isFailure() ) { std::cout << "Unable to open RawDataProviderSvc" << std::endl; }
  }
  return m_rawDataProviderSvc;
}