TrkReco.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TrkReco.cxx,v 1.100 2012/05/28 05:16:29 maoh Exp $
//-----------------------------------------------------------------------------
// Filename : TrkReco.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : [email protected]
//-----------------------------------------------------------------------------
// Description : A tracking module.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//               http://belle.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
#include <math.h> 
#include "CLHEP/String/Strings.h"
 
#include "MdcGeomSvc/MdcGeomSvc.h" 
#include "MdcTables/MdcTables.h"
 
#include "TrkReco/TrkReco.h"
#include "TrkReco/TMDC.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TPerfectFinder.h"
#include "TrkReco/TConformalFinder0.h"
#include "TrkReco/TConformalFinder.h"
#include "TrkReco/TCurlFinder.h"
#include "TrkReco/TFastFinder.h"
#include "TrkReco/TTrack.h"
#include "TrkReco/TTrackMC.h"
#include "TrkReco/TTrackHEP.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TTrackBase.h"
 
#include "EvTimeEvent/RecEsTime.h"
#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
#include "MdcCalibFunSvc/MdcCalibFunSvc.h"
 
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/IDataManagerSvc.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/IJobOptionsSvc.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/StatusCode.h"
 
#include "Identifier/Identifier.h"
#include "Identifier/MdcID.h"
#include "ReconEvent/ReconEvent.h" 
#include "EventModel/EventHeader.h"
#include "MdcRawEvent/MdcDigi.h"
#include "McTruth/McParticle.h"
#include "McTruth/MdcMcHit.h"
#include "RawEvent/RawDataUtil.h"
#include "RawDataProviderSvc/IRawDataProviderSvc.h"
 
#include "RawDataProviderSvc/MdcRawDataProvider.h"//jialk 20090624
 
 
#include <vector>
#include <iostream>
#include <fstream>
 
#include "BesTimerSvc/IBesTimerSvc.h"
#include "BesTimerSvc/BesTimerSvc.h"
#include "BesTimerSvc/BesTimer.h"
 
 
using namespace std;
using namespace Event;
 
int TrkRecoTest = 0;
float TrkRecoHelixFitterChisqMax = 0;
 
TrkReco::TrkReco(const std::string& name, ISvcLocator* pSvcLocator) :
  Algorithm(name, pSvcLocator),
  _cdc(0),
  _perfectFinder(0),
  _rkfitter("range fitter",
      false,
      0,
      true),
  _confFinder(0),
  _curlFinder(0),
  _nEvents(0) {
 
    /// Initiate Parameters for TrkReco
    initPara();
    
    TrkRecoHelixFitterChisqMax = b_helixFitterChisqMax;
}
 
//************************************************************************Initialize
StatusCode TrkReco::initialize(){
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "in initialize()" << endreq;
 
    StatusCode sc;
 
    /// Get RawDataProviderSvc
    IRawDataProviderSvc* irawDataProviderSvc;
    sc = service ("RawDataProviderSvc", irawDataProviderSvc);
    _rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
    if ( sc.isFailure() ){
      log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq;
      return StatusCode::FAILURE;
    }
 
    /// Get MdcGeomSvc
 //   IMdcGeomSvc* imdcGeomSvc;
 //   sc = service("MdcGeomSvc", imdcGeomSvc);
 //   _mdcGeomSvc =  dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
 //   if (sc.isFailure()) {
 //     return( StatusCode::FAILURE);  
 //   }
 
    /// Get MdcCalibFunSvc
    IMdcCalibFunSvc* imdcCalibSvc;
    sc = service ("MdcCalibFunSvc",imdcCalibSvc);
    _mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*> (imdcCalibSvc);
    if ( sc.isFailure() ){
      log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
    }
 
    /// Initialize NTuple
    if (b_tuple) InitTuple();
 
#ifdef TRKRECO_DEBUG
    b_debugLevel = 2;
#endif
    
    //...Create TMDC...
//    _cdc = cdcInit();
 
    //...TrkReco...
    TrkRecoHelixFitterChisqMax = b_helixFitterChisqMax;
 
    //...Create rphi finder...
/*    _perfectFinder = new TPerfectFinder(b_perfectFitting,
                    b_conformalMaxSigma,
                    b_conformalStereoMaxSigma,
                    b_conformalFittingCorrections);
*/
    if ((! _confFinder) && (b_conformalFinder == 0)) {
       _confFinder = new TConformalFinder0(b_conformalMaxSigma,
                        b_conformalFraction,
                        b_conformalStereoZ3,
                        b_conformalStereoZ4,
                        b_conformalStereoChisq3,
                        b_conformalStereoChisq4,
                        b_conformalStereoMaxSigma,
                        b_conformalFittingCorrections,
                        b_salvageLevel,
                        b_doConformalFinderCosmic);
    }
    else if (! _confFinder) {
    _confFinder = new TConformalFinder(b_doConformalFastFinder,
                       b_doConformalSlowFinder,
                       b_conformalPerfectSegmentFinding,
                       b_conformalMaxSigma,
                       b_conformalStereoMaxSigma,
                       b_salvageLevel,
                       b_conformalMinNLinksForSegment,
                       b_conformalMinNCores,
                       b_conformalMinNSegments,
                       b_conformalSalvageLoadWidth,
                       b_conformalStereoMode,
                       b_conformalStereoLoadWidth,
                       b_conformalStereoSzSegmentDistance,
                       b_conformalStereoSzLinkDistance,
                       b_conformalFittingFlag);
    }
    _confFinder->debugLevel(b_debugLevel);
    _confFinder->doStereo(b_doConformalFinderStereo);
 
    //...Salvage flag is ignored in new conf. finder...
    if (b_doSalvage == 2) _confFinder->doSalvage(true);
 
    //...Create curl finder...
    if (! _curlFinder)
        _curlFinder = new TCurlFinder((unsigned)min_segment,
                      (unsigned)min_salvage,
                      bad_distance_for_salvage,
                      good_distance_for_salvage,
                      (unsigned)min_sequence,
                      (unsigned)min_fullwire,
                      range_for_axial_search,
                      range_for_stereo_search,
                      (unsigned)superlayer_for_stereo_search,
                      range_for_axial_last2d_search,
                      range_for_stereo_last2d_search,
                      trace2d_distance,
                      trace2d_first_distance,
                      trace3d_distance,
                      (unsigned)determine_one_track,
                      selector_max_impact,
                      selector_max_sigma,
                      selector_strange_pz,
                      selector_replace_dz,
                      (unsigned)stereo_2dfind,
                      (unsigned)merge_exe,
                      merge_ratio,
                      merge_z_diff,
                      mask_distance,
                      ratio_used_wire,
                      range_for_stereo1,
                      range_for_stereo2,
                      range_for_stereo3,
                      range_for_stereo4,
                      range_for_stereo5,
                      range_for_stereo6,
                      z_cut,
                      z_diff_for_last_attend,
                      (unsigned)svd_reconstruction,
                      min_svd_electrons,
                      (unsigned)on_correction,
                      (unsigned)output_2dtracks,
                      (unsigned)curl_version,
                      //jialk
                      minimum_seedLength, 
                      minimum_2DTrackLength,
                      minimum_3DTrackLength,
                      minimum_closeHitsLength,
                      MIN_RADIUS_OF_STRANGE_TRACK, 
                      ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK);
    _curlFinder->debugLevel(b_debugLevel);
 
    //...Track manager setup...
    _trackManager.maxMomentum(b_momentumCut);
    _trackManager.debugLevel(b_debugLevel);
    _trackManager.fittingFlag(b_fittingFlag);
 
    //...Initialize...
//zsl    TUpdater::initialize();
 
    //...For debug...
//    dump("parameter");
 
    TrkRecoTest = b_test;
 
    if (b_timeTest && b_tuple) {
      StatusCode sc = service( "BesTimerSvc", m_timersvc);
      if( sc.isFailure() ) {
    log << MSG::WARNING << " Unable to locate BesTimerSvc" << endreq;
    return StatusCode::FAILURE;
      } 
      m_timer[1] = m_timersvc->addItem("Execution");
      m_timer[1]->propName("nExecution");
    }
    
    return StatusCode::SUCCESS;   
}
//************************************************************************MDCInit
TMDC *
TrkReco::cdcInit(void) {
    _cdcVersion = dtostring(b_cdcVersion);
    if (! _cdc) _cdc = TMDC::getTMDC(_cdcVersion);
    _cdc->debugLevel(b_debugLevel);
 
    //...Obtain fudge factor...
    float fudgeFactor = b_fudgeFactor;
 
/*   if (fudgeFactor == 0.) {
    const calcdc_const4 & c = * (calcdc_const4 *) BsGetEnt(CALMDC_CONST4,
                                   1,
                                   BBS_No_Index);
    if (!(& c)) fudgeFactor = 1.;
    else        fudgeFactor = c.m_fudge;
    }*/
   
    //...Turn off fudge factor before checking performance...
    fudgeFactor = 1.0;
 
    _cdc->fudgeFactor(fudgeFactor);
    return _cdc;
}
 
//************************************************************************Finalize
StatusCode TrkReco::finalize() {
 
    if (b_timeTest && b_tuple) m_timersvc->print();
 
    //...Clear TrkReco objects...
    if (_cdc) _cdc->fastClear();
    if (b_doConformalFinder && _confFinder) _confFinder->clear();
    if (b_doCurlFinder && _curlFinder) _curlFinder->clear();
 
    //...Summary...
//    dump("summary");
 
    //...Delete TrkReco objects...
    if (_confFinder) delete _confFinder;
    if (_curlFinder) delete _curlFinder;
 
    return StatusCode::SUCCESS;
}
StatusCode TrkReco::beginRun(){
  /// Get MdcGeomSvc
    IMdcGeomSvc* imdcGeomSvc;
    StatusCode sc1 =Gaudi::svcLocator()->service("MdcGeomSvc", imdcGeomSvc);
    _mdcGeomSvc =  dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
    if (sc1.isFailure()) {
      return( StatusCode::FAILURE);
    }
     _cdc = cdcInit();
 
    return StatusCode::SUCCESS;  
}
void
TrkReco::disp_stat(const char * m) {
    std::string msg = m;
//    dump(msg);
}
 
//************************************************************************Execute
StatusCode TrkReco::execute() {
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "in execute()" << endreq;
  StatusCode sc;   
  if (b_timeTest && b_tuple) m_timer[1]->start();
 
  // Initiate state of all sense wire
  if (b_goodTrk && b_tuple)
    for (int ii=0;ii<43;ii++){
      for (int jj=0;jj<288;jj++){
    havedigi[ii][jj]= -99;//no hit/noise
      }
    }
  
 //------------------------------------
  // Initialize track collection in TDS
  //------------------------------------
  //yzhang add 2010-04-30 
  //Clear TDS tracks and hits
  if(!b_RungeKuttaCorrection){
  IDataManagerSvc *dataManSvc;
  DataObject *aTrackCol;
  DataObject *aRecHitCol;
  if(!m_combineTracking){
    dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc().get());
    eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
    if(aTrackCol != NULL) {
      dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
      eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
    }
    eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
    if(aRecHitCol != NULL) {
      dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
      eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
    }
  }
  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);
    }
  }
  RecMdcTrackCol* trkcol;
  eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
  if (aTrackCol) {
    trkcol = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
  }else{
    trkcol= new RecMdcTrackCol;
    sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trkcol);
    if(!sc.isSuccess()) {
      log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
      return StatusCode::FAILURE;
    }
  }
  RecMdcHitCol* hitcol;
  eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
  if (aRecHitCol) {
    hitcol = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
  }else{
    hitcol= new RecMdcHitCol;
    sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitcol);
    if(!sc.isSuccess()) {
     log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
      return StatusCode::FAILURE;
    }
  }  
  // Part 1: Get Tev
  t0_bes = 0.;
  t0Sta  = -1;
  if (useESTime) { 
    SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
//bugcheck jialk    
    if (aevtimeCol  && aevtimeCol->size() ) {
      RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
      t0_bes =  (*iter_evt)->getTest();
      t0Sta  =  (*iter_evt)->getStat();
    }else{
      log << MSG::WARNING << "Could not find EsTimeCol" << endreq;
      return StatusCode::SUCCESS;
    }
  }
 
  // Part 2: Get the event header
  SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
  if (!eventHeader) {
    log << MSG::FATAL << "Could not find Event Header" << endreq;
    return( StatusCode::FAILURE);
  }
  _nEvents = eventHeader->eventNumber();
  if (b_tuple) std::cout << "TrkReco ... processing ev# " << _nEvents << std::endl;
 
  // Part 3: Retrieve MDC digi 
  int digiId;
  uint32_t getDigiFlag = 0;
  if(m_dropHot || m_combineTracking)getDigiFlag |= MdcRawDataProvider::b_dropHot;
  if(m_keepBadTdc)  getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
  if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
  MdcDigiVec mdcDigiVec = _rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
  if (0 == mdcDigiVec.size()){
    log << MSG::WARNING << " No hits in MdcDigiVec" << endreq;
    return StatusCode::SUCCESS;
  }
/*  //jialk in order to reject events with exceeding # of total hits   2009/03/31   
if (mdcDigiVec.size() > 2000){   
     log << MSG::WARNING << " Too many hits in MdcDigiVec" << endreq;   
     return StatusCode::SUCCESS;   
   }*/ 
  MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
  MC_DIGI_SIZE = mdcDigiVec.size();
  
  digiId = 0; 
  Identifier mdcId;
  int layerId;
  int wireId;
 
  //Clear the old MdcRec_wirhit tables and create the hits' info for the new event.
  unsigned nt = MdcRecWirhitCol::getMdcRecWirhitCol()->size();
  //cout<<"Col size of last Event's WirHit = "<<nt<<endl;
  MdcRecWirhitCol::getMdcRecWirhitCol()->clear();
 
  for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
    //log << MSG::INFO << "MDC digit No: " << digiId << endreq;
    mdcId = (*iter1)->identify();
    layerId = MdcID::layer(mdcId);
    wireId = MdcID::wire(mdcId);
    /*log << MSG::INFO
      << " time_channel = " << (*iter1)->getTimeChannel()
      <<" time = "<<(*iter1)->getTimeChannel()* 40./10000/1000000
      << " charge_channel = " << (*iter1)->getChargeChannel()
      << " layerId = " << layerId 
      << " wireId = " << wireId
      << endreq;*/
 
    if (b_goodTrk && b_tuple) havedigi[layerId][wireId] = (*iter1)->getTrackIndex();  //-1:noise  0-n:tracks
 
        MdcRec_wirhit mhit;
 
        mhit.id = digiId;
    mhit.geo = _mdcGeomSvc->Wire(layerId,wireId);
    //Apply crude TOF, Belle: tof=1.074*radius/c, here we use 1.28 instead of 1.074. 
    double tof;
    tof = 1.28 * mhit.geo->Lyr()->Radius()/299.8;  //unit of Radius is mm.  
    mhit.tdc = RawDataUtil::MdcTime((*iter1)->getTimeChannel());
 
    //cout<<" ..  mhit.tdc = "<<mhit.tdc<<"    t0_bes = "<<t0_bes<<endl;
    mhit.tdc -= t0_bes;
 
    //jialk 
    _trackManager.sett0bes(t0_bes);
    //mhit.adc = RawDataUtil::MdcCharge((*iter1)->getChargeChannel());
    //mhit.tdc = (*iter1)->getTimeChannel();
    mhit.adc = (*iter1)->getChargeChannel();
    //cout<<"raw tdc(ns): "<<mhit.tdc<<"; tof(ns): "<<tof<<endl;
//  double dist2 = _mdcCalibFunSvc->rawTimeNoTOFToDist(mhit.tdc-tof, layerId, wireId, 2, 0.0);
 //       mhit.erddl = _mdcCalibFunSvc->getSigma(layerId, 2, dist2, 0.0);
        double timewalk=_mdcCalibFunSvc->getTimeWalk(layerId,mhit.adc); 
        double T0 = _mdcCalibFunSvc->getT0(layerId,wireId);
        double drifttime = mhit.tdc-tof-timewalk-T0;
        double dist2 = _mdcCalibFunSvc->driftTimeToDist(drifttime,layerId, wireId, 2, 0);//by liucy 2010/05/12
        mhit.erddl = _mdcCalibFunSvc->getSigma(layerId, 2 , dist2, 0.0,0.0,0.0,mhit.adc);
//cout<<"getSigma: "<<mhit.erddl<<endl;
    mhit.ddl = dist2/10.;  //mm->cm
        mhit.ddr = mhit.ddl;
    mhit.erddl = mhit.erddl/10.;  //mm->cm
        mhit.erddr = mhit.erddl;
 
        mhit.lr = 2;
        mhit.stat = 0;
        mhit.stat = mhit.stat |= 1048576; //bit20  WireHitTimeValid
        mhit.stat = mhit.stat |= 2097152; //bit21  WireHitChargeValid
        mhit.stat = mhit.stat |= 4194304; //bit22  WireHitFindingValid
        mhit.stat = mhit.stat |= 1073741824; //bit30 
    //cout<<"layerNo = "<<mhit.geo->Layer()<<"; "<<mdigi.digi[i].layerNo<<endl;
    //cout<<"cellNo = "<<mhit.geo->Cell()<<"; "<<mdigi.digi[i].cellNo<<endl;
    //cout<<"NCell of this layer = "<<mhit.geo->Lyr()->NCell()<<endl;
    //cout<<"NCell of this layer = "<<mhit.geo->Lyr()->NCell()<<endl;
        MdcRecWirhitCol::getMdcRecWirhitCol()->push_back(mhit);
 
    if (b_tuple) {
    t10_tdc = mhit.tdc;
    t10_adc = mhit.adc;
    t10_drift = mhit.ddl;
    t10_dDrift = mhit.erddl;
    t10_lyrId = layerId;
    t10_localId = wireId;
    m_tuple10->write();
    }
   }
 
    unsigned nT0Reset = 0;
 
    //...Starting point...
TrkReco_start:
 
    //...Clear myself...
    clear();
    
    //jialk in order to reject events with exceeding # of total hits   2009/03/31   
    if (mdcDigiVec.size() > 2000){   
     log << MSG::WARNING << " Too many hits in MdcDigiVec" << endreq;   
     return StatusCode::SUCCESS;   
    }
 
    //...Update TMDC...
    _cdc->update(b_doMCAnalysis);
 
    //...Get lists of hits...
    unsigned mask = 0;
    if (! b_useAllHits) mask = WireHitFindingValid;
    const AList<TMDCWireHit> & axialHits = _cdc->axialHits(mask);
    const AList<TMDCWireHit> & stereoHits = _cdc->stereoHits(mask);
    const AList<TMDCWireHit> & allHits = _cdc->hits(mask);
    //cout<<"axial: "<<axialHits.length()<<"  stereo: "<<stereoHits.length()<<endl; 
    
    //...Storage for tracks...
    AList<TTrack> tracks;
    AList<TTrack> tracks2D;
 
    //...Perfect finder...
    if (b_doPerfectFinder) {
      _perfectFinder->doit(axialHits, stereoHits, tracks, tracks2D);
      _trackManager.append(tracks);
    }
 
    else {
    //...Conformal finder...
    if (b_doConformalFinder) {
 
        //...T0 reset option...
        if (b_doT0Reset) {
        if (b_nT0ResetMax > nT0Reset)
            ((TConformalFinder *) _confFinder)->doT0Reset(true);
        else
            ((TConformalFinder *) _confFinder)->doT0Reset(false);
        }           
 
        _confFinder->doit(axialHits, stereoHits, tracks, tracks2D);
 
        //...T0 reset...
        if (b_doT0Reset) {
        ++nT0Reset;
        if (((TConformalFinder *) _confFinder)->T0ResetDone())
            goto TrkReco_start;
        }
 
        //cout<<"tracks: "<<tracks.length()<<endl;
        
        //...Stores tracks...
        _trackManager.append(tracks);
        _trackManager.append2D(tracks2D);
        if (b_conformalFinder == 0) {
        if (b_doSalvage == 1) _trackManager.salvage(allHits);
        if (b_doSalvage) _trackManager.mask();
        }
    }
    
    //...Curl finder...
        if (b_doCurlFinder) {
        if ((! b_doSalvage) && (b_conformalFinder == 0))
        _trackManager.maskCurlHits(axialHits,
                       stereoHits,
                       _trackManager.tracks());
        AList<TTrack> confTracks = _trackManager.tracks();
        tracks.append(confTracks);
        _curlFinder->doit(axialHits, stereoHits, tracks, tracks2D);
            tracks.remove(confTracks);
        //_trackManager.append(tracks);
    }
 
    //...Finishes tracks...
//  if ((b_doSalvage) && (b_conformalFinder == 0)) _trackManager.refit();   
 
        //...Appends tracks which are reconstructed by CurlFinder...
    if (b_doCurlFinder) {
      _trackManager.append(tracks);
      _trackManager.append2D(tracks2D);
    }
       
    //...Merge & Mask ...
    //if ((b_doMerge) && (b_conformalFinder == 0)) _trackManager.merge();
    _trackManager.merge();  //Liuqg
    //if (b_conformalFinder != 0) _trackManager.setCurlerFlags();
 
    //...Salvage for dE/dx...
    //if (b_doAssociation)
    //  _trackManager.salvageAssociateHits(allHits, b_associateSigma);
    }
 
    //...Move a pivot... //move to the innermost hit, remove this step now 2005/10/17 zang shilei
    //_trackManager.movePivot();
 
    //...Save Panther tables...
    // _trackManager.checkNumberOfHits();
    //if (b_sortMode == 0) _trackManager.sortTracksByQuality();
    //else                 _trackManager.sortTracksByPt();
    /*if (b_doMCAnalysis) {
      if (mcEvent()) {
    mcInformation();
    _trackManager.saveMCTables();
      }
    }*/
 
    if (b_tuple) FillTuple();
 
    //write Tds 
    int t_tkStat = 2;
    if (!b_doConformalFinder && b_doCurlFinder) t_tkStat = 3;//FIXME
    StatusCode scTds = _trackManager.makeTds(trkcol, hitcol, t_tkStat,b_RungeKuttaCorrection,m_CalibFlag);
    if (scTds !=  StatusCode::SUCCESS) return( StatusCode::FAILURE);
 
    //...T0 correction...
    if (b_doT0Determination) {
      _trackManager.determineT0(b_doT0Determination, b_nTracksForT0);
      if (b_tuple) t3_t0Rec = _trackManager.paraT0();
    }
  }
  else if(b_RungeKuttaCorrection==1){
    if (useESTime) { 
      SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
      if (aevtimeCol  && aevtimeCol->size() ) {
    RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
    t0_bes =  (*iter_evt)->getTest();
    t0Sta  =  (*iter_evt)->getStat();
      }else{
    log << MSG::WARNING << "Could not find EsTimeCol" << endreq;
    return StatusCode::SUCCESS;
      }
    }
 
    SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
    if (!eventHeader) {
      log << MSG::FATAL << "Could not find Event Header" << endreq;
      return( StatusCode::FAILURE);
    }
    _nEvents = eventHeader->eventNumber();
    if (b_tuple) std::cout << "TrkReco ... processing ev# " << _nEvents << std::endl;
    AList<TTrack> rktracks;
    int digiId;
    uint32_t getDigiFlag = 0;
    if(m_dropHot || m_combineTracking)getDigiFlag |= MdcRawDataProvider::b_dropHot;
    if(m_keepBadTdc)  getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
    if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
    MdcDigiVec mdcDigiVec = _rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
    if (0 == mdcDigiVec.size()){
      log << MSG::WARNING << " No hits in MdcDigiVec" << endreq;
      return StatusCode::SUCCESS;
    }
    MdcDigiVec::iterator iter1 = mdcDigiVec.begin();
    MC_DIGI_SIZE = mdcDigiVec.size();
 
    digiId = 0;
    Identifier mdcId;
    int layerId;
    int wireId;
 
    unsigned nt = MdcRecWirhitCol::getMdcRecWirhitCol()->size();
    MdcRecWirhitCol::getMdcRecWirhitCol()->clear();
 
    for (;iter1 != mdcDigiVec.end(); iter1++, digiId++) {
      mdcId = (*iter1)->identify();
      layerId = MdcID::layer(mdcId);
      wireId = MdcID::wire(mdcId);
      if (b_goodTrk && b_tuple) havedigi[layerId][wireId] = (*iter1)->getTrackIndex();
 
      MdcRec_wirhit mhit;
 
      mhit.id = digiId;
      mhit.geo = _mdcGeomSvc->Wire(layerId,wireId);
      double tof;
      tof = 1.28 * mhit.geo->Lyr()->Radius()/299.8;
      mhit.tdc = RawDataUtil::MdcTime((*iter1)->getTimeChannel());
      mhit.timechannel=(*iter1)->getTimeChannel();
      mhit.tdc -= t0_bes;
 
      _trackManager.sett0bes(t0_bes);
      mhit.adc = (*iter1)->getChargeChannel();
      double timewalk=_mdcCalibFunSvc->getTimeWalk(layerId,mhit.adc);
      double T0 = _mdcCalibFunSvc->getT0(layerId,wireId);
      double drifttime = mhit.tdc-tof-timewalk-T0;
      double dist2 = _mdcCalibFunSvc->driftTimeToDist(drifttime,layerId, wireId, 2, 0);
      mhit.erddl = _mdcCalibFunSvc->getSigma(layerId, 2 , dist2, 0.0,0.0,0.0,mhit.adc);
      mhit.ddl = dist2/10.;
      mhit.ddr = mhit.ddl;
      mhit.erddl = mhit.erddl/10.;
      mhit.erddr = mhit.erddl;
 
      mhit.lr = 2;
      mhit.stat = 0;
      mhit.stat = mhit.stat |= 1048576;
      mhit.stat = mhit.stat |= 2097152;
      mhit.stat = mhit.stat |= 4194304;
      mhit.stat = mhit.stat |= 1073741824;
      MdcRecWirhitCol::getMdcRecWirhitCol()->push_back(mhit);
 
      if (b_tuple) {
    t10_tdc = mhit.tdc;
    t10_adc = mhit.adc;
    t10_drift = mhit.ddl;
    t10_dDrift = mhit.erddl;
    t10_lyrId = layerId;
    t10_localId = wireId;
    m_tuple10->write();
      }
    }
 
    unsigned nT0Reset = 0;
 
    clear();
    if (mdcDigiVec.size() > 2000){
      log << MSG::WARNING << " Too many hits in MdcDigiVec" << endreq;
      return StatusCode::SUCCESS;
    }
    int counter=0;
    _cdc->update(b_doMCAnalysis);
 
    unsigned mask = 0;
    if (! b_useAllHits) mask = WireHitFindingValid;
    const AList<TMDCWireHit> & axialHits = _cdc->axialHits(mask);
    const AList<TMDCWireHit> & stereoHits = _cdc->stereoHits(mask);
    const AList<TMDCWireHit> & allHits = _cdc->hits(mask);
    AList<TMLink> _allHits[3];
    TConformalFinder0::conformalTransformationDriftCircle(ORIGIN, axialHits, _allHits[0]);
    TConformalFinder0::conformalTransformationDriftCircle(ORIGIN, stereoHits, _allHits[1]);
    _allHits[2].append(_allHits[0]);
    _allHits[2].append(_allHits[1]);
 
 
   SmartDataPtr<RecMdcTrackCol> mdcTracks(eventSvc(),EventModel::Recon::RecMdcTrackCol);
   if( ! mdcTracks )
   {
     log << MSG::ERROR << "Unable to retrieve RecMdcTrackCol" << endreq;
     return  StatusCode::FAILURE;
   } else {
     log << MSG::DEBUG << "RecMdcTrackCol retrieved of size "<< mdcTracks->size() << endreq;
     int ntrk=0;
     for(RecMdcTrackCol::iterator it=mdcTracks->begin(); it!=mdcTracks->end(); it++)
     {ntrk++;
       TRunge *r = new TRunge(**it);
       TTrack *t1 = new TTrack(*r);
       HitRefVec gothits = (*it)->getVecHits();
       HitRefVec::iterator it_gothit = gothits.begin();
       unsigned stat=(*it)->stat();
       int nhit=0;
       for( ; it_gothit != gothits.end(); it_gothit++){
         for(int i=0;i<_allHits[2].length();i++){
           int lyrraw=_allHits[2][i]->wire()->layerId();
           int wireraw=_allHits[2][i]->wire()->localId();
           int g_layer       = MdcID::layer((*it_gothit)->getMdcId());
           int g_wire        = MdcID::wire((*it_gothit)->getMdcId());
           if(lyrraw==g_layer&&g_wire==wireraw){
             nhit++;
             t1->append(*_allHits[2][i]);
           }
         }
       }
       TRunge *rr = new TRunge(*t1);
       TRunge *tt = new TRunge(*t1);
       int err=_rkfitter.fit(*rr);
       int nhits=rr->links().length();
       int ndrop=0;
       int nmax=0;
       if(nhits<=10)nmax=0;
       if(nhits>10)nmax=(int)nhit*0.3;
       for(int ii=0;ii<nmax;ii++){ 
    
        ndrop=rr->DropWorst();
        if(ndrop)err=_rkfitter.fit(*rr);
        if(err)break;
       }  
       if(err==-2) counter++;
       if(m_CalibFlag==1) {//This is for Calibration 
         tt->removeLinks();
         tt->append(rr->links());
         for(int i=0;i<43;i++){
            err= _rkfitter.fit(*tt,0,i);
         }
        rr->removeLinks();
        rr->append(tt->links());
        TTrack *t =new TTrack(*rr);
        if(err==0){
            rktracks.append(*t); 
            t->setFinderType(100+stat);
        }
       }
       if(m_CalibFlag==0){
       TTrack *t =new TTrack(*rr);
       if(err==0){
           
        t->setFinderType(100+stat);
         rktracks.append(*t);
       }
       }
       delete r ;
       delete t1;
       delete rr;
       delete tt;
     }
   }
 
  _trackManager.append(rktracks);
  IDataManagerSvc *dataManSvc;
  DataObject *aTrackCol;
  DataObject *aRecHitCol;
  if(!m_combineTracking){
    dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc().get());
    eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
    if(aTrackCol != NULL) {
      dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
      eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
    }
    eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
    if(aRecHitCol != NULL) {
      dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
      eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
    }
  }
  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);
    }
  }
  RecMdcTrackCol* trkcol;
  eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
  if (aTrackCol) {
    trkcol = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
  }else{
    trkcol= new RecMdcTrackCol;
    sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trkcol);
    if(!sc.isSuccess()) {
      log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
      return StatusCode::FAILURE;
    }
  }
  RecMdcHitCol* hitcol;
  eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
  if (aRecHitCol) {
    hitcol = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
  }else{
    hitcol= new RecMdcHitCol;
    sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitcol);
    if(!sc.isSuccess()) {
      log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
      return StatusCode::FAILURE;
    }
  }
  if (b_tuple) FillTuple();
 
  int t_tkStat = 2;
  if (!b_doConformalFinder && b_doCurlFinder) t_tkStat = 3;//FIXME
  StatusCode scTds = _trackManager.makeTds(trkcol, hitcol, t_tkStat,b_RungeKuttaCorrection,m_CalibFlag);
  if (b_doT0Determination) {
    _trackManager.determineT0(b_doT0Determination, b_nTracksForT0);
    if (b_tuple) t3_t0Rec = _trackManager.paraT0();
  }
   for (unsigned i = 0; i < 3; i++) {
        if (i == 2)
            HepAListDeleteAll(_allHits[i]);
        else
            _allHits[i].removeAll();
    }
  
  rktracks.removeAll(); 
  }
  //...For debug...
  /*if (b_debugLevel) {
    std::cout << "TrkReco ... ev# " << _nEvents << " processed,"
    << " #tracks found=" << _trackManager.allTracks().length()
    << ", #good tracks=" << _trackManager.tracks().length()
    << ", #2D tracks=" << _trackManager.tracks2D().length()
    << std::endl;
    if (b_debugLevel > 1) _trackManager.dump("eventSummary hits");
    else _trackManager.dump("eventSummary");
    }*/
 
  //TUpdater::update();
  return StatusCode::SUCCESS;
}
 
void
TrkReco::mcInformation(void) {
 
  //...Preparation...
  //  const AList<TTrack> & allTracks = _trackManager.allTracks();
  const AList<TTrack> & allTracks = _trackManager.tracksFinal();
 
  unsigned nHep = TTrackHEP::list().length();
  unsigned nTrk = allTracks.length();
  unsigned * N = (unsigned *) malloc(nHep * sizeof(unsigned));
  for (unsigned i = 0; i < nHep; i++) N[i] = 0;
 
  //...Loop over all tracks...
  for (unsigned i = 0; i < nTrk; i++) {
    TTrackMC * mc = allTracks[i]->_mc;
    if (! mc) {
      mc = new TTrackMC(* allTracks[i]);
      _mcTracks.append(mc);
      allTracks[i]->_mc = mc;
    }
 
    mc->update();
    if (mc->hepId() != -1)
      if (mc->charge())
    ++N[mc->hepId()];
  }
 
  //...Check uniqueness...
  for (unsigned i = 0; i < nHep; i++) {
    if (N[i] < 2) {
      for (unsigned j = 0; j < nTrk; j++)
    if (allTracks[j]->_mc->hepId() == i)
      allTracks[j]->_mc->_quality += TTrackUnique;
    }
  }
 
  //...Good tracks...
  for (unsigned i = 0; i < nTrk; i++) {
    unsigned & quality = allTracks[i]->_mc->_quality;
    if ((quality & TTrackGhost) && (quality & TTrackUnique))
      quality += TTrackGood;
  }
 
  //...Termination...
  free(N);
}
 
void
TrkReco::clear(void) {
 
  //...Clear track candidates of the last event...
  HepAListDeleteAll(_mcTracks);
 
  //...Clear finders...
  if (b_doPerfectFinder) _perfectFinder->clear();
  if (b_doConformalFinder) _confFinder->clear();
  if (b_doCurlFinder) _curlFinder->clear();
  _trackManager.clear();
  _trackManager.clearTables();
}
 
void
TrkReco::fastClear(void) {
  clear();
}
 
bool
TrkReco::mcEvent(void) const {
  //    struct belle_event * ev = 
  //    (struct belle_event *) BsGetEnt(BELLE_EVENT, 1, BBS_No_Index);
 
  //...No BELLE_EVENT ???...
  //    if (! ev) return false;
  //    if (ev->m_ExpMC == 2) return true;
  //    return false;
  return true;
}
 
//std::string
//TrkReco::version(void) const {
//    return "TrkReco-00-00-04";
//}
 
void
TrkReco::initPara(void){
  _rkfitter.setBesFromGdml();
  declareProperty("mcPar", b_mcPar = 0);
  declareProperty("mcHit", b_mcHit = 0);
  declareProperty("tuple", b_tuple = 0);
  declareProperty("goodTrk", b_goodTrk = 0);
  declareProperty("timeTest", b_timeTest = 0);
 
  declareProperty("cdcVersion", b_cdcVersion = 1.0);
  declareProperty("fudgeFactor", b_fudgeFactor = 1.0);
 
  declareProperty("useESTime", useESTime = 1.0);
 
  declareProperty("debugLevel", b_debugLevel = 0);
  declareProperty("useAllHits", b_useAllHits = 0);
  declareProperty("doT0Reset", b_doT0Reset = 0);
  declareProperty("nT0ResetMax", b_nT0ResetMax = 1);
  declareProperty("doMCAnalysis", b_doMCAnalysis = 1);
  declareProperty("helixFitterChisqMax", b_helixFitterChisqMax = 0.);
 
  declareProperty("RungeKuttaCorrection", b_RungeKuttaCorrection = 0);
  declareProperty("doPerfectFinder", b_doPerfectFinder = 0);
  declareProperty("perfectFitting", b_perfectFitting = 0);
 
  declareProperty("conformalFinder", b_conformalFinder = 1);
  declareProperty("doConformalFinder", b_doConformalFinder = 0);
  declareProperty("doConformalFastFinder", b_doConformalFastFinder = 1);
  declareProperty("doConformalSlowFinder", b_doConformalSlowFinder = 1);
  declareProperty("conformalPerfectSegmentFinding", b_conformalPerfectSegmentFinding = 0);
  declareProperty("conformalFittingFlag", b_conformalFittingFlag = 4);  //1: sag  2: propagation  4: tof  8: freeT0
  declareProperty("conformalMaxSigma", b_conformalMaxSigma = 30.);
  declareProperty("conformalMinNLinksForSegment", b_conformalMinNLinksForSegment = 2);
  declareProperty("conformalMinNCores", b_conformalMinNCores = 2);    //min core for a seg
  declareProperty("conformalMinNSegments", b_conformalMinNSegments = 2); //min seg for trk, default: 3
  declareProperty("salvageLevel", b_salvageLevel = 30.);
  declareProperty("conformalSalvageLoadWidth", b_conformalSalvageLoadWidth = 1);
  declareProperty("conformalStereoMode", b_conformalStereoMode = 1);
  declareProperty("conformalStereoLoadWidth", b_conformalStereoLoadWidth = 5);
  declareProperty("conformalStereoMaxSigma", b_conformalStereoMaxSigma = 30.);
  declareProperty("conformalStereoSzSegmentDistance", b_conformalStereoSzSegmentDistance = 10.);
  declareProperty("conformalStereoSzLinkDistance", b_conformalStereoSzLinkDistance = 15.);
 
  declareProperty("doConformalFinderStereo", b_doConformalFinderStereo = 1);  //cosmic, use stereo
  declareProperty("doConformalFinderCosmic", b_doConformalFinderCosmic = 0);
  declareProperty("conformalFraction", b_conformalFraction = 0.7);      //cores fraction in trk
  declareProperty("conformalStereoZ3", b_conformalStereoZ3 = 20.);
  declareProperty("conformalStereoZ4", b_conformalStereoZ4 = 20.);
  declareProperty("conformalStereoChisq3", b_conformalStereoChisq3 = 30.);
  declareProperty("conformalStereoChisq4", b_conformalStereoChisq4 = 30.);
  declareProperty("conformalFittingCorrections", b_conformalFittingCorrections = 0);
 
  declareProperty("doCurlFinder", b_doCurlFinder = 1);  //CurlFinder on: 1  off: 0
  declareProperty("doSalvage", b_doSalvage = 2);
  declareProperty("doMerge", b_doMerge = 0);
  declareProperty("momentumCut", b_momentumCut = 0.0);
  declareProperty("doT0Determination", b_doT0Determination = 7);
  declareProperty("nTracksForT0", b_nTracksForT0 = 2);
  declareProperty("sortMode", b_sortMode = 0);
  declareProperty("test", b_test = 0);
  declareProperty("fittingFlag", b_fittingFlag = 0);
  declareProperty("doAssociation", b_doAssociation = 1);
  declareProperty("associateSigma", b_associateSigma = 60);
 
  declareProperty("dropHot",      m_dropHot= 0);//jialk 20090624
  declareProperty("combineTracking",m_combineTracking=0);
  declareProperty("keepBadTdc",m_keepBadTdc=0);
  declareProperty("keepUnmatch",m_keepUnmatch=0);
 
  declareProperty("CalibFlag",m_CalibFlag=0);
  //   For Curl Finder -->
  declareProperty("min_segment", min_segment = 5);   //min links for a segment
  declareProperty("min_salvage", min_salvage = 10);  //salvage while trk.links < min
  declareProperty("bad_distance_for_salvage", bad_distance_for_salvage = 1.0);
  declareProperty("good_distance_for_salvage", good_distance_for_salvage = 0.2);
  declareProperty("min_sequence", min_sequence = 6);
  declareProperty("min_fullwire", min_fullwire = 5);  //Belle: 7
  declareProperty("range_for_axial_search", range_for_axial_search = 1.5);
  declareProperty("range_for_stereo_search", range_for_stereo_search = 2.5);
  declareProperty("superlayer_for_stereo_search", superlayer_for_stereo_search = 3);
  declareProperty("range_for_axial_last2d_search", range_for_axial_last2d_search = 1.5);
  declareProperty("range_for_stereo_last2d_search", range_for_stereo_last2d_search = 2.0);
  declareProperty("trace2d_distance", trace2d_distance = 35.0);
  declareProperty("trace2d_first_distance", trace2d_first_distance = 0.5);
  declareProperty("trace3d_distance", trace3d_distance = 30.0);
  declareProperty("determine_one_track", determine_one_track = 0);
  declareProperty("selector_max_impact", selector_max_impact = 4.0);
  declareProperty("selector_max_sigma", selector_max_sigma = 36.0);
  declareProperty("selector_strange_pz", selector_strange_pz = 10.0);
  declareProperty("selector_replace_dz", selector_replace_dz = 15.0);
  declareProperty("stereo_2dfind", stereo_2dfind = 0);
  declareProperty("merge_exe", merge_exe = 1);
  declareProperty("merge_ratio", merge_ratio = 0.1);
  declareProperty("merge_z_diff", merge_z_diff = 10.0);
  declareProperty("mask_distance", mask_distance = 0.5);
  declareProperty("ratio_used_wire", ratio_used_wire = 0.3);
  declareProperty("range_for_stereo1", range_for_stereo1 = 2.4);
  declareProperty("range_for_stereo2", range_for_stereo2 = 2.7);
  declareProperty("range_for_stereo3", range_for_stereo3 = 2.9);
  declareProperty("range_for_stereo4", range_for_stereo4 = 3.4);
  declareProperty("range_for_stereo5", range_for_stereo5 = 4.1);
  declareProperty("range_for_stereo6", range_for_stereo6 = 5.0);  //Liuqg
  declareProperty("z_cut", z_cut = 50.0);
  declareProperty("z_diff_for_last_attend", z_diff_for_last_attend = 1.5);
  declareProperty("on_correction", on_correction = 4);   //1 sag, 2 propagation, 4 tof. definition changed in BES
  declareProperty("output_2dtracks", output_2dtracks = 1); //fill 2D track
  declareProperty("curl_version", curl_version = 0);
  //jialk
  declareProperty("minimum_seedLength", minimum_seedLength = 5);
  declareProperty("minimum_2DTrackLength", minimum_2DTrackLength = 7);
  declareProperty("minimum_3DTrackLength", minimum_3DTrackLength = 10);
  declareProperty("minimum_closeHitsLength", minimum_closeHitsLength = 5);
  declareProperty("MIN_RADIUS_OF_STRANGE_TRACK",MIN_RADIUS_OF_STRANGE_TRACK = 50);
  declareProperty("ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK",ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK = 5);
  //   <-- For Curl Finder 
}
 
void
TrkReco::InitTuple(void){
  m_tuple=ntupleSvc()->book("FILE101/rec3D",CLID_ColumnWiseTuple,"MdcRecEvent");
  m_tuple->addItem ("mc_phi",   t_mcphi);
  m_tuple->addItem ("mc_theta", t_mctheta);
  m_tuple->addItem ("mc_ptot",  t_mcptot);
  m_tuple->addItem ("mc_pt",    t_mcpt);
  m_tuple->addItem ("mc_pz",    t_mcpz);
  m_tuple->addItem ("mc_t0",    t_mct0);
  m_tuple->addItem ("nDigi",    t_nDigi);
 
  m_tuple->addItem ("dr",    t_dr);    
  m_tuple->addItem ("dz",    t_dz);   
  m_tuple->addItem ("pt",    t_pt);    
  m_tuple->addItem ("ptot",  t_ptot);        
  m_tuple->addItem ("tanlmd",t_tanlmd);
  m_tuple->addItem ("phi",   t_phi);
  m_tuple->addItem ("radius",t_radius);
  m_tuple->addItem ("chi2",  t_chi2);
  m_tuple->addItem ("nHits", t_nHits);
  m_tuple->addItem ("nCores",t_nCores);
  m_tuple->addItem ("nSegs", t_nSegs);
  m_tuple->addItem ("ndf",   t_ndf);
 
  m_tuple->addItem ("dpt",   t_dpt);
  m_tuple->addItem ("dptot", t_dptot);
  m_tuple->addItem ("dlmd",  t_dlmd);
  m_tuple->addItem ("dphi",  t_dphi);
 
  m_tuple->addItem ("t0",    t_t0);
  m_tuple->addItem ("t0_sta",t0_sta);  //mdc, tof
 
  m_tuple->addItem ("evtNo", t_evtNo);
  m_tuple->addItem ("length", t_length);
  m_tuple->addItem ("length2", t_length2);
 
  m_tuple->addItem ("gd_theta", t_good_theta);
  m_tuple->addItem ("gd_nLayers", t_gdNLayers);
  m_tuple->addItem ("mc_nLayers", t_mcNLayers);
  m_tuple->addItem ("best_nLayers",t_bestNLayers);
  m_tuple->addItem ("best_mc_nLayers",t_bestMcNLayers);
 
  m_tuple3=ntupleSvc()->book("FILE101/raw",CLID_ColumnWiseTuple,"Raw Data");
  m_tuple3->addItem ("mc_t0",    t3_mct0);
  m_tuple3->addItem ("mc_theta", t3_mctheta);
  m_tuple3->addItem ("mc_phi",   t3_mcphi);
  m_tuple3->addItem ("mc_ptot",  t3_mcptot);
  m_tuple3->addItem ("mc_pt",    t3_mcpt);
  m_tuple3->addItem ("mc_pid",   t3_mcpid);
  m_tuple3->addItem ("evtNo",    t3_evtNo);
 
  m_tuple31=ntupleSvc()->book("FILE101/rawEvt",CLID_ColumnWiseTuple,"Raw Data");
  m_tuple31->addItem ("nDigi",  t3_nDigi);
  m_tuple31->addItem ("goodLength", t3_goodLength);
  m_tuple31->addItem ("length", t3_length);
  m_tuple31->addItem ("t0_rec", t3_t0Rec);
  m_tuple31->addItem ("t0", t3_t0);
  m_tuple31->addItem ("t0_sta", t3_t0Sta);
  m_tuple31->addItem ("finalLength", t3_finalLength);
 
  m_tuple31->addItem ("mc_theta", t3_mctheta);
  m_tuple31->addItem ("mc_ptot",  t3_mcptot);
  m_tuple31->addItem ("mc_pt",    t3_mcpt);
  m_tuple31->addItem ("evtNo",    t3_evtNo);
 
  m_tuple2=ntupleSvc()->book("FILE101/rec2D",CLID_ColumnWiseTuple,"MdcRecEvent 2D");
  m_tuple2->addItem ("mc_theta", t2_mctheta);
  m_tuple2->addItem ("mc_pt",    t2_mcpt);
  m_tuple2->addItem ("nDigi",   t2_nDigi);
  m_tuple2->addItem ("nHits",   t2_nHits);
  m_tuple2->addItem ("nSegs",   t2_nSegs);
  m_tuple2->addItem ("chi2",    t2_chi2);
  m_tuple2->addItem ("evtNo",   t2_evtNo);
  m_tuple2->addItem ("ndf",     t2_ndf);
  m_tuple2->addItem ("length",  t2_length);
  m_tuple2->addItem ("length2", t2_length2);
  m_tuple2->addItem ("radius",  t2_radius);
 
  m_tuple4=ntupleSvc()->book("FILE101/hit",CLID_ColumnWiseTuple,"MdcRecEvent Hits");
  m_tuple4->addItem ("d_cal", t4_Dist);
  m_tuple4->addItem ("d_meas", t4_drift);   //d_cal-d_meas
  m_tuple4->addItem ("e_meas", t4_dDrift);  //error measure
  m_tuple4->addItem ("mc_drift", t4_mcDrift);
  m_tuple4->addItem ("mcLR", t4_mcLR);
  m_tuple4->addItem ("pull", t4_pull);
  m_tuple4->addItem ("lyrId", t4_lyrId);
  m_tuple4->addItem ("localId", t4_localId);
  m_tuple4->addItem ("LR", t4_LR);
  m_tuple4->addItem ("tdc", t4_tdc);
  m_tuple4->addItem ("z", t4_z);
  m_tuple4->addItem ("bz", t4_bz);  //backward
  m_tuple4->addItem ("fz", t4_fz);  //forward
  m_tuple4->addItem ("fy", t4_fy);  //forward
  m_tuple4->addItem ("phi", t4_phi);
  m_tuple4->addItem ("nHits", t4_nHits);
 
  m_tuple5=ntupleSvc()->book("FILE101/char",CLID_ColumnWiseTuple,"MdcRecEvent Charge");
  m_tuple5->addItem ("ptotPos", t5_ptotPos);
  m_tuple5->addItem ("ptotNeg", t5_ptotNeg);
  m_tuple5->addItem ("drPos", t5_drPos);
  m_tuple5->addItem ("drNeg", t5_drNeg);
  m_tuple5->addItem ("dzPos", t5_dzPos);
  m_tuple5->addItem ("dzNeg", t5_dzNeg);
 
  m_tuple6=ntupleSvc()->book("FILE101/urec",CLID_ColumnWiseTuple,"MdcRecEvent urec");
  m_tuple6->addItem ("length2",  u_length2);
  m_tuple6->addItem ("mc_ptot",   u_mcptot);
  m_tuple6->addItem ("mc_pt",     u_mcpt);
  m_tuple6->addItem ("mc_theta",  u_mctheta);
  m_tuple6->addItem ("nDigi",    u_nDigi);
  m_tuple6->addItem ("evtNo",    u_evtNo);
  m_tuple6->addItem ("mc_t0",     u_mct0);
  m_tuple6->addItem ("t0",       ut_t0);
  m_tuple6->addItem ("t0_sta",   ut0_sta);
 
  if (b_timeTest && b_tuple) {
    m_tuple7=ntupleSvc()->book("FILE101/time",CLID_ColumnWiseTuple,"MdcRecEvent time");
    m_tuple7->addItem ("time", ti_eventTime);   //total time in a event.
    m_tuple7->addItem ("recNum", ti_recTrkNum); //total raw-tracks, include decayed tracks!
    m_tuple7->addItem ("evtNo", ti_evtNo);
    m_tuple7->addItem ("nHits", ti_nHits);  //total hits of rec-tracks
    m_tuple7->addItem ("nDigi", ti_nDigi);  //total hits in rawdata
  }
 
  m_tuple9=ntupleSvc()->book("FILE101/seg",CLID_ColumnWiseTuple,"MdcRecEvent segments");
  m_tuple9->addItem ("times", t9_times);
  m_tuple9->addItem ("nLinks", t9_nLinks);
  m_tuple9->addItem ("nUsed", t9_nUsed);
  m_tuple9->addItem ("nSL", t9_nSL);
  m_tuple9->addItem ("mctheta", t9_mctheta);
 
  m_tuple10=ntupleSvc()->book("FILE101/rawHit",CLID_ColumnWiseTuple,"MdcRecEvent mchitCOL");
  m_tuple10->addItem ("tdc", t10_tdc);
  m_tuple10->addItem ("adc", t10_adc);
  m_tuple10->addItem ("Drift", t10_drift);
  m_tuple10->addItem ("dDrift", t10_dDrift);
  m_tuple10->addItem ("lyrId", t10_lyrId);
  m_tuple10->addItem ("localId", t10_localId);
}
 
void
TrkReco::FillTuple(void) {
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "fill Tuple()" << endreq;
 
  /// Get Tracks
  AList<TTrack> tracks;
  AList<TTrack> tracks2D;
  tracks = _trackManager.tracks();
  tracks2D = _trackManager.tracks2D();
  //if(t3_t0Rec!=999. && fabs(t3_t0Rec + t0_bes - MC_EVENT_TIME) > 4) 
  //cout<<"3Dtrack in manager: "<<tracks.length()<<endl;
  if(tracks.length()==0) cout<<"zslength: 3D length=0, and the 2D length is"<<tracks2D.length()<<endl;
 
  /// Run Time Calculation  
  if (b_timeTest && b_tuple) { 
    m_timer[1]->stop();
    ti_eventTime = m_timer[1]->elapsed();
    ti_nDigi = MC_DIGI_SIZE;
    ti_recTrkNum = tracks.length();
    ti_evtNo = _nEvents;
    for (unsigned i = 0; i < ti_recTrkNum; ++i)
      ti_nHits += tracks[i]->nLinks();
    m_tuple7->write();
  }
 
  /// Retrieve MC
  CLHEP::Hep3Vector MC_TRACK_VEC;
  CLHEP::HepLorentzVector MC_TRACK_LRZVEC;
  float  MC_TRACK_NUM;
  double MC_EVENT_TIME;
 
  int digiId;
  int num_par = 0;
  MC_EVENT_TIME = -1;
  if (b_mcPar) {
    SmartDataPtr<Event::McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
    if (mcParticleCol) {
      Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
      digiId = 0;
      for (;iter_mc != mcParticleCol->end(); iter_mc++, digiId++) {
    if((*iter_mc)->statusFlags()==8320||(*iter_mc)->statusFlags()==128) {
      MC_EVENT_TIME = (*iter_mc)->initialFourMomentum().t();
      break;
    }
      }
      //jialk bugcheck
      iter_mc = mcParticleCol->begin();
      MC_TRACK_LRZVEC = (*iter_mc)->initialFourMomentum();  //particularly for 1 track events.
      MC_TRACK_VEC = (*iter_mc)->initialFourMomentum().vect();
      MC_TRACK_NUM = mcParticleCol->size();
 
#ifdef TRKRECO_DEBUG
      iter_mc = mcParticleCol->begin();
      digiId = 0;
      for (;iter_mc != mcParticleCol->end(); iter_mc++, digiId++) {
    log << MSG::INFO << "MDC digit No: " << digiId << endreq;
    log << MSG::INFO
      << " particleId = " << (*iter_mc)->particleProperty()
      << " theta = " << (*iter_mc)->initialFourMomentum().theta()
      <<"  phi= "<< (*iter_mc)->initialFourMomentum().phi()
      <<"  px= "<< (*iter_mc)->initialFourMomentum().px()
      <<"  py= "<< (*iter_mc)->initialFourMomentum().py()
      <<"  pz= "<< (*iter_mc)->initialFourMomentum().pz()
      << endreq;
      }
#endif
      digiId = 0;
      iter_mc = mcParticleCol->begin();
      for (;iter_mc != mcParticleCol->end(); iter_mc++, digiId++) {
    CLHEP::HepLorentzVector LRZVEC = (*iter_mc)->initialFourMomentum();
    CLHEP::Hep3Vector  VEC = (*iter_mc)->initialFourMomentum().vect();
    //staFlag = (*iter_mc)->statusFlags();
 
    t3_mctheta = LRZVEC.theta();
    t3_mcphi   = VEC.phi();
    t3_mcptot  = VEC.mag()/1000.;
    t3_mcpt    = VEC.perp()/1000.;
    t3_mct0    = (*iter_mc)->initialFourMomentum().t();
    t3_mcpid   = (*iter_mc)->particleProperty();
    t3_evtNo   = _nEvents;
    if((t3_mcpid==211||t3_mcpid==-211||t3_mcpid==321)&&fabs(cos(t3_mctheta))<0.93) ++num_par;
    m_tuple3->write();
      }
    }
  }
 
  /// Fill Tuple
  if (tracks.length()==0) {
    u_length2 = tracks2D.length();
    u_mct0    = MC_EVENT_TIME;
    u_mcptot  = MC_TRACK_VEC.mag()/1000.;
    u_mcpt    = MC_TRACK_VEC.perp()/1000.;
    u_mctheta = MC_TRACK_LRZVEC.theta();
    u_nDigi   = MC_DIGI_SIZE;
    u_evtNo   = _nEvents;
    ut0_sta = -1;
    ut_t0 = t0_bes;
    ut0_sta = t0Sta;
    m_tuple6->write();  //unrec...tracks
  }
 
  //---Pos and Neg
  float dDot = 2;
  int flagi = -1, flagj = -1;
  HepVector3D ppp1, ppp2;
  float pppdDot = 2;
  for (unsigned i = 0; i < tracks.length(); i++){
    if(tracks.length()<2) break;
    if(tracks[i]->charge()>0) {
      ppp1 = tracks[i]->p().unit();
      for (unsigned j = 0; j < tracks.length(); j++){
    if (tracks[j]->charge()<0) {
      ppp2 = tracks[j]->p().unit();
      pppdDot = ppp1.dot(ppp2);
      if (pppdDot<dDot) {
        flagi = i;
        flagj = j;
        dDot  = pppdDot;
      }
    }
      }
    }
  }
  if (flagi != -1 && flagj != -1 && dDot < 0) {
    t5_ptotPos = tracks[flagi]->ptot();
    t5_ptotNeg = tracks[flagj]->ptot();
    t5_drPos   = tracks[flagi]->helix().dr();
    t5_drNeg   = tracks[flagj]->helix().dr();
    t5_dzPos   = tracks[flagi]->helix().dz();
    t5_dzNeg   = tracks[flagj]->helix().dz();
    m_tuple5->write();
  }  //Pos and Neg...back on back
 
  unsigned nGood = 0;   
  for (unsigned i = 0; i < tracks.length(); i++){
    for(unsigned k = 0; k < tracks[i]->links().length(); k++){ 
      t4_Dist  = tracks[i]->links()[k]->distance();     //_onTrack - _onWire
      t4_drift = tracks[i]->links()[k]->drift();    //drift updated in THelixFitter.  
      t4_dDrift= tracks[i]->links()[k]->dDrift();       //dDrift updated in THelixFitter.
      t4_pull  = tracks[i]->links()[k]->pull();
      t4_LR    = 2;  //initial
      t4_LR    = tracks[i]->links()[k]->leftRight();
      if (t4_LR == 0) t4_LR = -1;
      t4_lyrId = tracks[i]->links()[k]->wire()->layerId();
      t4_localId = tracks[i]->links()[k]->wire()->localId();
      t4_tdc = tracks[i]->links()[k]->hit()->reccdc()->tdc;
      t4_z = tracks[i]->links()[k]->positionOnTrack().z();
      t4_bz = tracks[i]->links()[k]->wire()->backwardPosition().z();
      t4_fz = tracks[i]->links()[k]->wire()->forwardPosition().z();
      t4_fy = tracks[i]->links()[k]->wire()->forwardPosition().y();
      t4_nHits = tracks[i]->links().length();
      unsigned lyrID = tracks[i]->links()[k]->wire()->layerId(); 
      float phi0 = _cdc->layer(lyrID)->offset();
      int nWir = (int) _cdc->layer(lyrID)->nWires();
      t4_phi = phi0 + t4_localId*2*pi/nWir;
      if(t4_phi<0) t4_phi+=2*pi;
 
      if (b_mcHit) { 
    SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
    if (mcMdcMcHitCol) {
      Event::MdcMcHitCol::iterator iter_mchit1 = mcMdcMcHitCol->begin();
      digiId=0;
      for (;iter_mchit1 != mcMdcMcHitCol->end(); iter_mchit1++, digiId++) {
        const Identifier ident = (*iter_mchit1)->identify();
        if(MdcID::layer(ident) != t4_lyrId) continue;
        if(MdcID::wire(ident)==t4_localId) {
          t4_mcDrift = (*iter_mchit1)->getDriftDistance();  //drift in MC.
          t4_mcLR = (*iter_mchit1)->getPositionFlag();
          if (t4_mcLR == 0) t4_mcLR = -1;
          break;
        }
      }
    }
      }
      //cout<<" lyrId, localId: "<<t4_lyrId<<"  "<<t4_localId<<endl;
      //cout<<" tdc: "<<t4_tdc<<"  mcDrift: " << t4_mcDrift<<"  Drift: "<<t4_drift<<"  t4_Dist: "<<t4_Dist<<endl;
      m_tuple4->write();  //rec Hits
    }
 
    unsigned segSize = tracks[i]->segments().length();
    for(unsigned kk = 0; kk < segSize; ++kk) {
      unsigned segLength = tracks[i]->segments()[kk]->links().length();
      AList<TMLink> ll = tracks[i]->segments()[kk]->links();
      int nSmall = 0;
      for(unsigned nn = 0; nn < segLength; ++nn) {
    double tmpX = ll[nn]->positionD().x();
    double tmpY = ll[nn]->positionD().y();
    double tmpA = tracks[i]->segments()[kk]->lineTsf().x();
    double tmpB = tracks[i]->segments()[kk]->lineTsf().y();
    double tmpC = tracks[i]->segments()[kk]->lineTsf().z();
    double dis = fabs(tmpA * tmpX + tmpB * tmpY + tmpC) / sqrt(tmpA * tmpA + tmpB * tmpB);
    double idealDis = maxdDistance(ll[nn]);
    if(fabs(dis-ll[nn]->cDrift())/idealDis<0.001 && nSmall<2) {
      ++nSmall;
      continue;
    }
    t9_times += fabs(dis-ll[nn]->cDrift())/idealDis;
      }
      t9_nSL   = ll[0]->wire()->superLayerId();
      t9_nLinks = segLength;
      t9_mctheta = MC_TRACK_LRZVEC.theta();
      t9_times = t9_times / (t9_nLinks - nSmall);
      m_tuple9->write();  //seg info
    }
 
    t_mcphi   = MC_TRACK_VEC.phi();
    t_mctheta = MC_TRACK_LRZVEC.theta();
    t_mcptot  = MC_TRACK_VEC.mag()/1000.;
    t_mcpt    = MC_TRACK_VEC.perp()/1000.;
    t_mcpz    = MC_TRACK_LRZVEC.pz()/1000.;
    t_mct0    = MC_EVENT_TIME;
    t_nDigi   = MC_DIGI_SIZE;
 
    t_dr      = tracks[i]->helix().dr();
    t_dz      = tracks[i]->helix().dz();
    t_pt      = tracks[i]->pt();
    t_ptot    = tracks[i]->ptot();
    t_tanlmd  = tracks[i]->helix().tanl();
    t_phi     = tracks[i]->helix().phi0();
    t_chi2    = tracks[i]->chi2();
    t_nHits   = tracks[i]->nLinks();
    t_nCores  = tracks[i]->nCores();
    t_nSegs   = tracks[i]->segments().length();
    t_ndf     = tracks[i]->ndf();
    t_radius  = tracks[i]->radius();
    t_evtNo   = _nEvents;
    t_length  = tracks.length();
    t_length2 = tracks2D.length();
 
    t_dpt     = tracks[i]->pt() - t_mcpt;
    t_dptot   = tracks[i]->ptot() - t_mcptot;
    t_dlmd    = atan(tracks[i]->helix().tanl()) - (pi/2 - t_mctheta);
    t_dphi    = tracks[i]->helix().phi0() - t_mcphi;
 
    t_t0 = t0_bes;
    t0_sta = t0Sta;
 
    if (b_mcPar) {
      if (b_goodTrk && b_tuple) {
    unsigned mcTrackSize = MC_TRACK_NUM;
    unsigned recTrackSize = tracks.length();
    const unsigned matchSize = 20;
    //  if(recTrackSize>40 || mcTrackSize>40) matchSize = 50;
 
    int mLayers[matchSize];
    for (int ii=0; ii<matchSize; ii++)
      mLayers[ii] = 0;
    for(int jj = 0; jj < 43; ++jj) {
      int tmp[matchSize];
      for(unsigned kk = 0; kk < matchSize; ++kk) 
        tmp[kk] = 0;
      for(int kk = 0; kk < 288; ++kk) {
        if (havedigi[jj][kk]<0) continue;
        tmp[havedigi[jj][kk]] = 1;
      }
      for(int kk = 0; kk < matchSize; ++kk)
        if(tmp[kk] == 1) ++mLayers[kk];
    }
 
    unsigned trackSize = tracks[i]->nLinks();
    int trkIndex[43];
    int nLayers[matchSize];
    for (int j = 0; j < 43; ++j)
      trkIndex[j] = -99;
    for (int j = 0; j < matchSize; ++j)
      nLayers[j] = 0;
 
    for (int j = 0; j < trackSize; ++j) {
      unsigned lId = tracks[i]->links()[j]->wire()->layerId();
      unsigned loId = tracks[i]->links()[j]->wire()->localId();
      if (havedigi[lId][loId] >= 0) trkIndex[lId] = havedigi[lId][loId];
    }
    for (int j = 0; j < 43; ++j)
      if (trkIndex[j] >= 0) ++nLayers[trkIndex[j]];
 
    for (int j = 0; j < matchSize; ++j) {
      //cout<<"nLayers: "<<nLayers[j]<<"  mLayers:"<<mLayers[j]<<endl;
      if (nLayers[j]==0) continue;
      if ((float)nLayers[j]/(float)mLayers[j] > 0.51) {
        t_gdNLayers = nLayers[j];
        t_mcNLayers = mLayers[j];
        t_good_theta = MC_TRACK_LRZVEC.theta();
        ++nGood;
        break;
      }
    }
    if (t_good_theta == 0.) {
      int tmpLayers = 0;
      int tmpi = -1;
      for (int j = 0; j < matchSize; ++j) {
        if (nLayers[j]==0) continue;
        if (nLayers[j] > tmpLayers) {
          tmpLayers = nLayers[j];
          tmpi = j;
        }
      }
      if (tmpi != -1) {
        t_bestNLayers = nLayers[tmpi];
        t_bestMcNLayers = mLayers[tmpi];
      }
      else {
        t_bestNLayers = -1;
        t_bestMcNLayers = -1;
      }
    }
    else {
      t_bestNLayers = -2;
      t_bestMcNLayers = -2;
    }
      } //if b_goodTrk
    }
    m_tuple->write();  //rec track
  }
 
  t3_mctheta = MC_TRACK_LRZVEC.theta();
  t3_mcptot  = MC_TRACK_VEC.mag()/1000.;
  t3_mcpt    = MC_TRACK_VEC.perp()/1000.;
 
  t3_nDigi   = MC_DIGI_SIZE;
  t3_t0 = t0_bes;
  t3_t0Sta = t0Sta;
  t3_goodLength = nGood;
  t3_length = tracks.length();
  t3_finalLength = num_par;
  m_tuple31->write();    //raw...
 
  for (unsigned i = 0; i < tracks2D.length(); i++) {
    t2_mctheta = MC_TRACK_LRZVEC.theta();
    t2_mcpt    = MC_TRACK_VEC.perp()/1000.;
 
    t2_nDigi   = MC_DIGI_SIZE;
    t2_nHits   = tracks2D[i]->nLinks();
    t2_nSegs   = tracks2D[i]->segments().length();
    t2_chi2    = tracks2D[i]->chi2();
    t2_ndf     = tracks2D[i]->ndf();
    t2_radius  = tracks2D[i]->radius();
    t2_evtNo   = _nEvents;
    t2_length  = tracks.length();
    t2_length2 = tracks2D.length();
    m_tuple2->write();   //unused 2D tracks
  }
}
 
double
TrkReco::maxdDistance(TMLink *l) const{
  unsigned sl = l->wire()->superLayerId();
  unsigned lId = l->wire()->layerId();
  double _averR[11] = {9.7, 14.5, 22.14, 28.62, 35.1, 42.39, 48.87, 55.35, 61.83, 69.12, 74.79};
  double _averR2[43] = { 7.885, 9.07, 10.29, 11.525,
    12.72, 13.875, 15.01, 16.16,
    19.7, 21.3, 22.955, 24.555,
    26.215, 27.82, 29.465, 31.06,
    32.69, 34.265, 35.875, 37.455,
    39.975, 41.52, 43.12, 44.76,
    46.415, 48.02, 49.685, 51.37,
    53.035, 54.595, 56.215, 57.855,
    59.475, 60.995, 62.565, 64.165,
    66.68, 68.285, 69.915, 71.57,
    73.21, 74.76, 76.345};
  double radius = _averR2[lId];
  const double singleSigma = l->dDrift();
  return 2 * singleSigma / (radius * radius);
}
 
/*void
  TrkReco::dump(const std::string & msg, const std::string & pre) const {
  bool def = (msg == "") ? true : false;
 
  if (msg.find("summary") != -1 || msg.find("detail") != -1) {
  std::cout << "TrkReco Summary (" << version() << ")" << std::endl;
  std::cout << "    Track Manager Summary" << std::endl;
  _trackManager.dump("summary", "        ");
  if (_confFinder) {
  std::cout << "    Conformal Finder Summary" << std::endl;
  _confFinder->dump("summary", "        ");
  }
  }
  if (def || msg.find("parameter") != std::string::npos || msg.find("detail") != std::string::npos) {
  std::string t0 = pre;
  std::string t1 = pre + "    ";
  std::string t2 = pre + "        ";
 
  std::cout << t0 << name() << "(" << version() << ")";
  std::cout << " : debug level = " << b_debugLevel << std::endl;
  if (b_doPerfectFinder != 0) {
  std::cout << t1 << "Perfect finder is active. ";
  std::cout << "Other finders will not called." << std::endl;
  }
  if (b_conformalFinder == 1) {
  std::cout << t1 << "New TrkReco(new conf. finder) is active.";
  std::cout << std::endl;
  }
 
  std::cout << t1 << _cdc->name() << "(" << _cdc->version() << ")"
  << std::endl;
  std::cout << t2 << "cdc version           : " << _cdc->cdcVersion()
  << std::endl;
  std::cout << t2 << "fudge factor          : " << _cdc->fudgeFactor()
  << std::endl;
 
//  if (_cdccat)
//      std::cout << t1 << _cdccat->name() << "("
//            << _cdccat->version() << ")" << std::endl;
 
std::cout << t1 << "ignore hit quality      : " << b_useAllHits
<< std::endl;
std::cout << t1 << "do T0 reset             : " << b_doT0Reset
<< std::endl;
std::cout << t1 << "    max number of reset : " << b_nT0ResetMax
<< std::endl;
std::cout << t1 << "MC analysis             : " << b_doMCAnalysis
<< std::endl;
std::cout << t1 << "helix fitter chisq max  : "
<< b_helixFitterChisqMax << std::endl;
std::cout << t1 << "test flag               : " << b_test << std::endl;
 
 
std::cout << t1 << _trackManager.name();
std::cout << "(" << _trackManager.version() << ") options"
<< std::endl;
std::cout << t2 << "T0 determination      : " << b_doT0Determination
<<std::endl;
std::cout << t2 << "    # of tracks       : " << b_nTracksForT0
<< std::endl;
std::cout << t2 << "bank sort             : " << b_sortMode
<< std::endl;
std::cout << t2 << "max. momentum allowed : " << b_momentumCut
<< " GeV/c";
std::cout << std::endl;
std::cout << t2 << "salvage               : " << b_doSalvage
<< std::endl;
std::cout << t2 << "salvage level         : " << b_salvageLevel
<< std::endl;
 
std::cout << t1 << _perfectFinder->name();
std::cout << "(" << _perfectFinder->version() << ") options";
std::cout << std::endl;
std::cout << t2 << "do finding            : " << b_doPerfectFinder;
std::cout << std::endl;
std::cout << t2 << "perfect fitting       : " << b_perfectFitting;
std::cout << std::endl;
 
std::cout << t1 << _confFinder->name();
std::cout << "(" << _confFinder->version() << ") options" << std::endl;
std::cout << t2 << "do finding            : " << b_doConformalFinder;
std::cout << std::endl;
if (b_conformalFinder == 1) {
  std::cout << t2 << "fast finder           : ";
  std::cout << b_doConformalFastFinder << std::endl;
  std::cout << t2 << "slow finder           : ";
  std::cout << b_doConformalSlowFinder << std::endl;
  std::cout << t2 << "perfect segment find  : ";
  std::cout << b_conformalPerfectSegmentFinding << std::endl;
  std::cout << t2 << "max sigma             : ";
  std::cout << b_conformalMaxSigma << std::endl;
  std::cout << t2 << "min # hits for segment: ";
  std::cout << b_conformalMinNLinksForSegment << std::endl;
  std::cout << t2 << "min # cores in segment: ";
  std::cout << b_conformalMinNCores << std::endl;
  std::cout << t2 << "min # segments        : ";
  std::cout << b_conformalMinNSegments << std::endl;
  std::cout << t2 << "salvage level         : " << b_salvageLevel
    << std::endl;
  std::cout << t2 << "salvage load width    : ";
  std::cout << b_conformalSalvageLoadWidth << std::endl;
  std::cout << t2 << "stereo mode           : "
    << b_conformalStereoMode << std::endl;
  std::cout << t2 << "stereo load width     : ";
  std::cout << b_conformalStereoLoadWidth << std::endl;
  std::cout << t2 << "stereo max sigma      : ";
  std::cout << b_conformalStereoMaxSigma << std::endl;
  std::cout << t2 << "stereo segment dist.  : ";
  std::cout << b_conformalStereoSzSegmentDistance << std::endl;
  std::cout << t2 << "stereo link distance  : ";
  std::cout << b_conformalStereoSzLinkDistance << std::endl;
}
else {
  std::cout << t2 << "Old TrkReco(old conf. finder) is active."
    << std::endl;
  std::cout << t2 << "do stereo finding     : ";
  std::cout << b_doConformalFinderStereo;
  std::cout << std::endl;
  std::cout << t2 << "use cosmic builder    : ";
  std::cout << b_doConformalFinderCosmic;
  std::cout << std::endl;
  std::cout << t2 << "max sigma             : ";
  std::cout << b_conformalMaxSigma <<std::endl;
  std::cout << t2 << "fraction              : ";
  std::cout << b_conformalFraction <<std::endl;
  std::cout << t2 << "fitting corrections   : ";
  std::cout << b_conformalFittingCorrections << std::endl;
  std::cout << t2 << "stereo max sigma      : ";
  std::cout << b_conformalStereoMaxSigma;
  std::cout << std::endl;
  std::cout << t2 << "stereo z3             : ";
  std::cout << b_conformalStereoZ3 <<std::endl;
  std::cout << t2 << "stereo z4             : ";
  std::cout << b_conformalStereoZ4 <<std::endl;
  std::cout << t2 << "stereo chisq 3        : ";
  std::cout << b_conformalStereoChisq3;
  std::cout << std::endl;
  std::cout << t2 << "stereo chisq 4        : ";
  std::cout << b_conformalStereoChisq4;
  std::cout << std::endl;
}
std::cout << t1 << _curlFinder->name();
std::cout << "(" << _curlFinder->version() << ") options" << std::endl;
std::cout << t2 << "do finding            : " << b_doCurlFinder
<< std::endl;
 
if (_clustFinder) {
  std::cout << t1 << _clustFinder->name();
  std::cout << "(" << _clustFinder->version() << ") options"
    << std::endl;
  std::cout << t2 << "do finding            : " << b_doClustFinder
    << std::endl;
  std::cout << t2 << "cathode window        : " << b_cathodeWindow
    << std::endl;
  std::cout << t2 << "systematic correction : "
    << b_cathodeSystematics<<std::endl;
  std::cout << t2 << "cathode cosmic switch : " << b_cathodeCosmic
    << std::endl;
}
 
std::cout << t1 << "Patten Matching CurlFinder"
<< "(0.2beta) options" << std::endl;
std::cout << t2 << "do finding                    : "
<< b_doPMCurlFinder << std::endl;
std::cout << t2 << "mode(1=rec, 2=map, 3=plot)    : "
  << b_pmCurlFinder << std::endl;
  std::cout << t2 << "min electrons of svd clusters : "
  << min_svd_electrons_in_pmc << std::endl;
  std::cout << t1 << "SVD Associator" << "(0.2beta) options"
  << std::endl;
  std::cout << t2 << "doing : " << b_doSvdAssociator << std::endl;
 
  }
if (def || msg.find("tracks") != std::string::npos || msg.find("detail") != std::string::npos) {
  _trackManager.dump("eventSummary");
}
 
if (msg.find("cathode") != std::string::npos || msg.find("detail") != std::string::npos) {
  if (b_doClustFinder && _cdccat) {
    _cdccat->dump("hit");
    _clustFinder->dump("");
  }
}
}*/