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    : [email protected]
//-----------------------------------------------------------------------------
// 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;}  
    /// 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*(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());
           
           ///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" << 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;
   
 ///check the result:MdcTrackCol   
#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;
}