TrkReco.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TrkReco.cxx,v 1.103 2022/01/28 22:14:13 maqm Exp $
//-----------------------------------------------------------------------------
// Filename : TrkReco.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A tracking module.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//               http://belle.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------
 
// #include "CLHEP/String/Strings.h"
#include <math.h>
 
#include "MdcGeomSvc/IMdcGeomSvc.h"
#include "MdcTables/MdcTables.h"
 
#include "TrkReco.h"
#include "TrkReco/TConformalFinder.h"
#include "TrkReco/TConformalFinder0.h"
#include "TrkReco/TCurlFinder.h"
#include "TrkReco/TFastFinder.h"
#include "TrkReco/TMDC.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TPerfectFinder.h"
#include "TrkReco/TTrack.h"
#include "TrkReco/TTrackBase.h"
#include "TrkReco/TTrackHEP.h"
#include "TrkReco/TTrackMC.h"
 
#include "EvTimeEvent/RecEsTime.h"
#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
// #include "MdcCalibFunSvc/MdcCalibFunSvc.h"
 
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IDataManagerSvc.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/StatusCode.h"
 
#include "EventModel/EventHeader.h"
#include "Identifier/Identifier.h"
#include "Identifier/MdcID.h"
#include "McTruth/McParticle.h"
#include "McTruth/MdcMcHit.h"
#include "MdcRawEvent/MdcDigi.h"
#include "RawDataProviderSvc/IRawDataProviderSvc.h"
#include "RawEvent/RawDataUtil.h"
#include "ReconEvent/ReconEvent.h"
 
#include "RawDataProviderSvc/MdcRawDataProvider.h" //jialk 20090624
 
#include <fstream>
#include <iostream>
#include <vector>
 
#include "BesTimerSvc/IBesTimerSvc.h"
// #include "BesTimerSvc/BesTimerSvc.h"
#include "BesTimerSvc/BesTimer.h"
 
using namespace std;
using namespace Event;
 
int   TrkRecoTest                = 0;
float TrkRecoHelixFitterChisqMax = 0;
DECLARE_COMPONENT( TrkReco )
TrkReco::TrkReco( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator )
    , _cdc( 0 )
    , _perfectFinder( 0 )
    , _rkfitter( "range fitter", false, 0, true )
    // , _rkfitter( 0 )
    , _confFinder( 0 )
    , _curlFinder( 0 )
    , _nEvents( 0 ) {
 
  std::cout << std::fixed << std::setprecision( 3 ); // yzhang debug for 720
  /// Initiate Parameters for TrkReco
  initPara();
 
  TrkRecoHelixFitterChisqMax = b_helixFitterChisqMax;
}
 
//************************************************************************Initialize
StatusCode TrkReco::initialize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in initialize()" << endmsg;
 
  _rkfitter.setBesFromGdml();
 
  StatusCode sc;
 
  /// Get RawDataProviderSvc
  IRawDataProviderSvc* irawDataProviderSvc;
  sc                  = service( "RawDataProviderSvc", irawDataProviderSvc );
  _rawDataProviderSvc = irawDataProviderSvc;
  if ( sc.isFailure() )
  {
    log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endmsg;
    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 = imdcCalibSvc;
  if ( sc.isFailure() ) { log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endmsg; }
 
  /// 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" << endmsg;
      return StatusCode::FAILURE;
    }
    m_timer[1] = m_timersvc->addItem( "Execution" );
    m_timer[1]->propName( "nExecution" );
  }
 
  return StatusCode::SUCCESS;
}
//************************************************************************MDCInit
TMDC* TrkReco::cdcInit( void ) {
  _cdcVersion = std::to_string( 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      = 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() {
  if ( !m_beginRun )
  {
    StatusCode sc = beginRun();
    if ( sc.isFailure() )
    {
      error() << "beginRun failed" << endmsg;
      return StatusCode::FAILURE;
    }
    m_beginRun = true;
  }
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in execute()" << endmsg;
 
  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 ( nullptr != aTrackCol )
      {
        dataManSvc->clearSubTree( "/Event/Recon/RecMdcTrackCol" );
        eventSvc()->unregisterObject( "/Event/Recon/RecMdcTrackCol" );
      }
      eventSvc()->findObject( "/Event/Recon/RecMdcHitCol", aRecHitCol );
      if ( nullptr != aRecHitCol )
      {
        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" << endmsg;
        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" << endmsg;
        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" << endmsg;
        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" << endmsg;
        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" << endmsg;
      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" << endmsg;
      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" << endmsg;
         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 << endmsg;
      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
        << endmsg;*/
 
      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" << endmsg;
      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" << endmsg;
        return StatusCode::SUCCESS;
      }
    }
 
    SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
    if ( !eventHeader )
    {
      log << MSG::FATAL << "Could not find Event Header" << endmsg;
      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" << endmsg;
      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" << endmsg;
      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" << endmsg;
      return StatusCode::FAILURE;
    }
    else
    {
      log << MSG::DEBUG << "RecMdcTrackCol retrieved of size " << mdcTracks->size() << endmsg;
      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 )
      {
        dataManSvc->clearSubTree( "/Event/Recon/RecMdcTrackCol" );
        eventSvc()->unregisterObject( "/Event/Recon/RecMdcTrackCol" );
      }
      eventSvc()->findObject( "/Event/Recon/RecMdcHitCol", aRecHitCol );
      if ( aRecHitCol )
      {
        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" << endmsg;
        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" << endmsg;
        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" << endmsg;
        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 ) {
  /// Initiate Parameters for TrkReco
  // _rkfitter = new TRungeFitter( "range fitter", false, 0, true );
 
  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()" << endmsg;
 
  /// 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 << endmsg;
        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() << endmsg;
      }
#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("");
  }
}
}*/