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

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TTrackManager.cxx,v 1.68.14.1 2013/03/04 07:22:52 xuqn 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 <values.h>
#include <cstring>
#include <cstdlib>

//#include "belle.h"
#include "CLHEP/String/Strings.h"
//#include "basf/basfshm.h"
#include "TrkReco/TrkReco.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TTrack.h"
#include "TrkReco/TTrackMC.h"
#include "TrkReco/TTrackHEP.h"
#include "TrkReco/TTrackManager.h"
#include "MdcTables/MdcTables.h"
#include "MdcTables/TrkTables.h"

#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/IDataManagerSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/IJobOptionsSvc.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/StatusCode.h"

#include "GaudiKernel/ContainedObject.h"
#include "GaudiKernel/SmartRef.h"
#include "GaudiKernel/SmartRefVector.h"
#include "GaudiKernel/ObjectVector.h"
#include "EventModel/EventModel.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) {
//    BASF_Sharedmem->allocate("TrkMgrSum", sizeof(struct summary));

  StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF); 
  if(scmgn!=StatusCode::SUCCESS) { 
    std::cout<< "ERROR: Unable to open Magnetic field service"<<std::endl;
  }
  // Get RawDataProviderSvc
  IRawDataProviderSvc* irawDataProviderSvc;
  StatusCode sc = Gaudi::svcLocator()->service("RawDataProviderSvc",irawDataProviderSvc); 
  m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
  if(sc!=StatusCode::SUCCESS) { 
    std::cout<< "ERROR: Unable to load RawDataProviderSvc"<<std::endl;
  }
}

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" << endreq;
  } else {
     log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
    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" <<endreq;
//       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;}  
    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*(m_pmgnIMF->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());
           
           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" << endreq;
//   return trksc; 
// }
//
// StatusCode hitsc = eventSvc->registerObject("/Event/Recon/RecMdcHitCol", hitcol);
// if ( hitsc.isFailure() ) {
//   log << MSG::FATAL << "Could not register MdcRecHit" << endreq;
//   return hitsc;
// }
//log << MSG::INFO << "MdcTrackCol registered successfully!" <<endreq;
   
#ifdef TRKRECO_DEBUG
  SmartDataPtr<RecMdcTrackCol> newtrkCol(eventSvc,"/Event/Recon/RecMdcTrackCol");
  if (!newtrkCol) { 
     log << MSG::FATAL << "Could not find RecMdcTrackCol" << endreq;
     return( StatusCode::FAILURE);
  }
  log << MSG::INFO << "Begin to check RecMdcTrackCol"<<endreq; 
   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 = m_rawDataProviderSvc->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) << endreq
          << "kappa: " << (*iter_trk)->getCpa() << "   Error of kappa "
          << (*iter_trk)->getError()(3,3) << endreq
          << "Tanl: " << (*iter_trk)->getTanl() << "   Error of Tanl "
          << (*iter_trk)->getError()(5,5) << endreq      
          << "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()
          << endreq;

          log << MSG::DEBUG << "hitList of this  track:" << endreq;
          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()
          << endreq      
          << "   hits TDC " <<(*it_gothit)->getTdc()  
          << "   hits ADC " <<(*it_gothit)->getAdc() << endreq;
          } 
          */
   }
#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" + itostring(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" + itostring(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" << endreq;
  }

  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" <<endreq;
  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" << endreq;
      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;
    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;
}

---