MdcxTrackFinder.cxx

Go to the documentation of this file.

//--------------------------------------------------------------------------
// File and Version Information:
//  $Id: MdcxTrackFinder.cxx,v 1.51 2022/02/21 04:34:53 maqm Exp $
//
// Description:
//  Class MdcxTrackFinder
//
// Environment:
//  Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Author List:
//  Steve Wagner            Original Author
//  Zhang Yao(zhangyao@ihep.ac.cn)
//
// Copyright Information:
//  Copyright (C) 1997      BEPCII
//
// History:
//  Migration for BESIII MDC
//
//------------------------------------------------------------------------
 
#include "GaudiKernel/IDataManagerSvc.h"
 
//-----------------------
// This Class's Header --
//-----------------------
#include "MdcxTrackFinder.h"
 
//-------------
// C Headers --
//-------------
#include <assert.h>
#include <math.h>
#include <time.h>
 
//---------------
// C++ Headers --
//---------------
#include <iostream>
 
//-------------------------------
// Collaborating Class Headers --
//-------------------------------
#include "EventModel/EventHeader.h"
#include "GaudiKernel/MsgStream.h"
#include "MdcGeom/MdcDetector.h"
#include "MdcRawEvent/MdcDigi.h"
#include "PatBField/BField.h"
 
#include "CLHEP/Alist/AIterator.h"
#include "EvTimeEvent/RecEsTime.h"
#include "Identifier/Identifier.h"
#include "Identifier/MdcID.h"
#include "McTruth/McParticle.h"
#include "McTruth/MdcMcHit.h"
#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
#include "MdcRecEvent/RecMdcHit.h"
#include "MdcRecEvent/RecMdcTrack.h"
#include "MdcRecoUtil/Pdt.h"
#include "MdcTrkRecon/MdcTrack.h"
#include "MdcxReco/MdcxAddHits.h"
#include "MdcxReco/MdcxFindSegs.h"
#include "MdcxReco/MdcxFindTracks.h"
#include "MdcxReco/MdcxFittedHel.h"
#include "MdcxReco/MdcxHel.h"
#include "MdcxReco/MdcxHit.h"
#include "MdcxReco/MdcxHits.h"
#include "MdcxReco/MdcxMergeDups.h"
#include "MdcxReco/MdcxSeg.h"
#include "MdcxReco/Mdcxprobab.h"
#include "RawDataProviderSvc/IRawDataProviderSvc.h"
#include "RawDataProviderSvc/MdcRawDataProvider.h"
#include "RawEvent/RawDataUtil.h"
#include "ReconEvent/ReconEvent.h"
#include "TrkBase/TrkErrCode.h"
#include "TrkBase/TrkExchangePar.h"
#include "TrkBase/TrkFit.h"
#include "TrkBase/TrkFitStatus.h"
#include "TrkBase/TrkHitList.h"
#include "TrkBase/TrkRecoTrk.h"
#include "TrkFitter/TrkHelixMaker.h"
#include "TrkFitter/TrkLineMaker.h"
// debug
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/NTuple.h"
#include "MdcGeom/Constants.h"
#include "MdcGeom/MdcTrkReconCut.h"
#include "MdcxReco/MdcxHistItem.h"
#include "MdcxReco/MdcxSegPatterns.h"
#include "TrigEvent/TrigData.h"
#include "TrkBase/TrkHotList.h"
 
#include "BesTimerSvc/BesTimer.h"
#include "BesTimerSvc/IBesTimerSvc.h"
 
using MdcPatRec::BField;
using std::cout;
using std::endl;
 
extern double MdcTrkReconCut_helix_fit[43];
//----------------
// Constructors --
//----------------
DECLARE_COMPONENT( MdcxTrackFinder )
MdcxTrackFinder::MdcxTrackFinder( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ), m_mdcCalibFunSvc( 0 ) {
  // input
  declareProperty( "pdtFile", m_pdtFile = "pdt.table" );
  // debug control
  declareProperty( "debug", m_debug = 0 );
  declareProperty( "hist", m_hist = 0 );
  declareProperty( "mcHist", m_mcHist = false );
  // cuts and control
  declareProperty( "cresol", m_cresol = 0.013 );
 
  declareProperty( "getDigiFlag", m_getDigiFlag = 0 );
  declareProperty( "maxMdcDigi", m_maxMdcDigi = 0 );
  declareProperty( "keepBadTdc", m_keepBadTdc = 0 );
  declareProperty( "dropHot", m_dropHot = 0 );
  declareProperty( "keepUnmatch", m_keepUnmatch = 0 );
  declareProperty( "salvageTrk", m_salvageTrk = false );
  declareProperty( "dropMultiHotInLayer", m_dropMultiHotInLayer = false );
  declareProperty( "dropTrkPt", m_dropTrkPt = -999. );
  declareProperty( "d0Cut", m_d0Cut = 999. );
  declareProperty( "z0Cut", m_z0Cut = 999. );
 
  declareProperty( "minMdcDigi", m_minMdcDigi = 0 );
  declareProperty( "countPropTime", m_countPropTime = true );
  declareProperty( "addHitCut", m_addHitCut = 5. );
  declareProperty( "dropHitsSigma", m_dropHitsSigma );
  declareProperty( "helixFitCut", m_helixFitCut );
  declareProperty( "minTrkProb", m_minTrkProb = 0.01 );
  declareProperty( "csmax4", m_csmax4 = 50. );
  declareProperty( "csmax3", m_csmax3 = 1. );
  declareProperty( "helixFitSigma", m_helixFitSigma = 5. );
  declareProperty( "maxRcsInAddSeg", m_maxRcsInAddSeg = 50. );
  declareProperty( "nSigAddHitTrk", m_nSigAddHitTrk = 5. );
  declareProperty( "maxProca", m_maxProca = 0.6 );
  declareProperty( "doSag", m_doSag = false );
  declareProperty( "lineFit", m_lineFit = false );
  // declareProperty("cosmicFit",     m_cosmicFit= false);
}
 
//--------------
// Destructor --
//--------------
MdcxTrackFinder::~MdcxTrackFinder() {
  if ( m_bfield ) delete m_bfield;
}
 
StatusCode MdcxTrackFinder::beginRun() {
  // Get Mdc Detector Geometry
  m_gm = MdcDetector::instance( m_doSag );
  if ( NULL == m_gm ) return StatusCode::FAILURE;
  MdcxHit::setMdcDetector( m_gm );
 
  return StatusCode::SUCCESS;
}
//--------------
// Operations --
//--------------
StatusCode MdcxTrackFinder::initialize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in initialize()" << endmsg;
 
  t_nTkTot = 0;
  for ( int i = 0; i < 20; i++ ) t_nTkNum[i] = 0;
 
  // m_flags.readPar(m_paramFile);
#ifdef MDCXTIMEDEBUG
  StatusCode tsc = service( "BesTimerSvc", m_timersvc );
  if ( tsc.isFailure() )
  {
    log << MSG::WARNING << name() << ": Unable to locate BesTimer Service" << endmsg;
    return StatusCode::FAILURE;
  }
  m_timer[0] = m_timersvc->addItem( "Execution" );
  m_timer[0]->propName( "Execution" );
  m_timer[1] = m_timersvc->addItem( "findSeg" );
  m_timer[1]->propName( "findSeg" );
  m_timer[2] = m_timersvc->addItem( "findTrack" );
  m_timer[2]->propName( "findTrack" );
  m_timer[3] = m_timersvc->addItem( "fitting" );
  m_timer[3]->propName( "fitting" );
#endif
 
  if ( m_helixFitCut.size() == 43 )
  {
    for ( int i = 0; i < 43; i++ )
    {
      // MdcTrkReconCut_helix_fit[i] = m_helixFitCut[i];
      TrkHelixFitter::nSigmaCut[i] = m_helixFitCut[i];
    }
  }
  MdcxParameters::debug          = m_debug;
  MdcxParameters::minTrkProb     = m_minTrkProb;
  MdcxParameters::csmax4         = m_csmax4;
  MdcxParameters::csmax3         = m_csmax3;
  MdcxParameters::helixFitSigma  = m_helixFitSigma;
  MdcxParameters::maxRcsInAddSeg = m_maxRcsInAddSeg;
  MdcxParameters::nSigAddHitTrk  = m_nSigAddHitTrk;
  MdcxParameters::maxProca       = m_maxProca;
  TrkHelixFitter::m_debug        = ( m_debug > 7 );
  Pdt::readMCppTable( m_pdtFile );
  MdcxFittedHel::debug = m_debug;
 
  // Get MdcCalibFunSvc
  StatusCode sc = service( "MdcCalibFunSvc", m_mdcCalibFunSvc );
  if ( sc.isFailure() )
  {
    log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endmsg;
    return StatusCode::FAILURE;
  }
  MdcxHit::setMdcCalibFunSvc( m_mdcCalibFunSvc );
  MdcxHit::setCountPropTime( m_countPropTime );
 
  // Get RawDataProviderSvc
  sc = service( "RawDataProviderSvc", m_rawDataProviderSvc );
  if ( sc.isFailure() )
  {
    log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endmsg;
    return StatusCode::FAILURE;
  }
 
  // Initailize magnetic filed
  sc = service( "MagneticFieldSvc", m_pIMF );
  if ( sc != StatusCode::SUCCESS )
  { log << MSG::ERROR << "Unable to open Magnetic field service" << endmsg; }
  m_bfield = new BField( m_pIMF );
  log << MSG::INFO << "field z = " << m_bfield->bFieldNominal() << endmsg;
  m_context = new TrkContextEv( m_bfield );
 
  if ( m_hist ) { bookNTuple(); }
  if ( m_dropHitsSigma.size() == 43 )
  {
    for ( int ii = 0; ii < 43; ii++ )
    { MdcxParameters::dropHitsSigma[ii] = m_dropHitsSigma[ii]; }
  }
 
  return StatusCode::SUCCESS;
}
 
StatusCode MdcxTrackFinder::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;
 
  b_saveEvent = false;
  setFilterPassed( b_saveEvent );
 
#ifdef MDCXTIMEDEBUG
  m_timer[0]->start();
  m_timer[1]->start();
#endif
 
  nTk      = 0;
  t_nTdsTk = 0;
  t_nDigi  = 0;
  t_nSeg   = 0; // yzhang  for fill
  //------------------------------------
  // Get event No.
  //------------------------------------
  SmartDataPtr<Event::EventHeader> evtHead( eventSvc(), "/Event/EventHeader" );
  if ( !evtHead )
  {
    log << MSG::FATAL << "Could not retrieve event header" << endmsg;
    return StatusCode::FAILURE;
  }
  m_eventNo = evtHead->eventNumber();
  if ( m_debug > 0 )
    std::cout << "x evt: " << evtHead->runNumber() << " " << m_eventNo << std::endl;
  long t_evtNo = m_eventNo;
  g_eventNo    = m_eventNo;
  // if (t_evtNo % 1000 == 0) std::cout << "x evt: " << t_evtNo << std::endl;
  IDataManagerSvc* dataManSvc;
  DataObject*      aTrackCol;
  DataObject*      aHitCol;
  if ( !m_salvageTrk )
  {
    SmartIF<IDataManagerSvc> dataManSvc( eventSvc() );
    eventSvc()->findObject( "/Event/Recon/RecMdcTrackCol", aTrackCol );
    if ( aTrackCol )
    {
      dataManSvc->clearSubTree( "/Event/Recon/RecMdcTrackCol" );
      eventSvc()->unregisterObject( "/Event/Recon/RecMdcTrackCol" );
    }
    eventSvc()->findObject( "/Event/Recon/RecMdcHitCol", aHitCol );
    if ( aHitCol )
    {
      dataManSvc->clearSubTree( "/Event/Recon/RecMdcHitCol" );
      eventSvc()->unregisterObject( "/Event/Recon/RecMdcHitCol" );
    }
  }
 
  //------------------------------------
  // Initialize track collection in TDS
  //------------------------------------
  DataObject* aReconEvent;
  eventSvc()->findObject( "/Event/Recon", aReconEvent );
  if ( aReconEvent == NULL )
  {
    aReconEvent = new ReconEvent();
    sc          = eventSvc()->registerObject( "/Event/Recon", aReconEvent );
    if ( sc != StatusCode::SUCCESS )
    {
      log << MSG::FATAL << "Could not register ReconEvent" << endmsg;
      return StatusCode::FAILURE;
    }
  }
  RecMdcTrackCol* trackList;
  eventSvc()->findObject( "/Event/Recon/RecMdcTrackCol", aTrackCol );
  if ( aTrackCol ) { trackList = dynamic_cast<RecMdcTrackCol*>( aTrackCol ); }
  else
  {
    trackList = new RecMdcTrackCol;
    sc        = eventSvc()->registerObject( EventModel::Recon::RecMdcTrackCol, trackList );
    if ( !sc.isSuccess() )
    {
      log << MSG::FATAL << " Could not register RecMdcTrack collection" << endmsg;
      return StatusCode::FAILURE;
    }
  }
  RecMdcHitCol* hitList;
  eventSvc()->findObject( "/Event/Recon/RecMdcHitCol", aHitCol );
  if ( aHitCol ) { hitList = dynamic_cast<RecMdcHitCol*>( aHitCol ); }
  else
  {
    hitList = new RecMdcHitCol;
    sc      = eventSvc()->registerObject( EventModel::Recon::RecMdcHitCol, hitList );
    if ( !sc.isSuccess() )
    {
      log << MSG::FATAL << " Could not register RecMdcHit collection" << endmsg;
      return StatusCode::FAILURE;
    }
  }
 
  //------------------------------------
  // Initialize hit collection in TDS
  //------------------------------------
  DataObject* pnode = 0;
  sc                = eventSvc()->retrieveObject( "/Event/Hit", pnode );
  if ( !sc.isSuccess() )
  {
    pnode = new DataObject;
    sc    = eventSvc()->registerObject( "/Event/Hit", pnode );
    if ( !sc.isSuccess() )
    {
      log << MSG::FATAL << " Could not register /Event/Hit branch " << endmsg;
      return StatusCode::FAILURE;
    }
  }
  SmartDataPtr<MdcHitCol> m_hitCol( eventSvc(), "/Event/Hit/MdcHitCol" );
  if ( !m_hitCol )
  {
    m_hitCol = new MdcHitCol;
    sc       = eventSvc()->registerObject( "/Event/Hit/MdcHitCol", m_hitCol );
    if ( !sc.isSuccess() )
    {
      log << MSG::FATAL << " Could not register hit collection" << endmsg;
      return StatusCode::FAILURE;
    }
  }
 
  //------------------------------------
  // Get bunch time t0 (ns) and timing
  //------------------------------------
  m_bunchT0 = -999.;
  SmartDataPtr<RecEsTimeCol> aevtimeCol( eventSvc(), "/Event/Recon/RecEsTimeCol" );
  if ( !aevtimeCol || aevtimeCol->size() == 0 )
  {
    log << MSG::WARNING << "evt " << m_eventNo << " Could not find RecEsTimeCol" << endmsg;
    return StatusCode::SUCCESS;
  }
 
  RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
  for ( ; iter_evt != aevtimeCol->end(); iter_evt++ )
  {
    m_bunchT0 = ( *iter_evt )->getTest();
    m_t0Stat  = ( *iter_evt )->getStat();
    if ( m_debug > 1 )
      std::cout << name() << " " << t_evtNo << " t0  " << m_bunchT0 << " t0Stat  " << m_t0Stat
                << std::endl;
    if ( ( m_t0Stat == 0 ) || ( m_bunchT0 < 0. ) || ( m_bunchT0 > 9999.0 ) )
    {
      log << MSG::WARNING << "Skip evt:" << m_eventNo << " by  t0 = " << m_bunchT0 << endmsg;
      // return StatusCode::SUCCESS;
    }
  }
  if ( m_debug > 1 )
    std::cout << name() << " evt " << t_evtNo << " t0  " << m_bunchT0 << " t0Stat  "
              << m_t0Stat << std::endl;
  int                    trigtiming = -10;
  SmartDataPtr<TrigData> trigData( eventSvc(), "/Event/Trig/TrigData" );
  if ( trigData )
  {
    log << MSG::INFO << "Trigger conditions 0--43:" << endmsg;
    for ( int i = 0; i < 48; i++ )
    {
      log << MSG::INFO << trigData->getTrigCondName( i ) << " ---- "
          << trigData->getTrigCondition( i ) << endmsg;
    }
    for ( int i = 0; i < 16; i++ )
      log << MSG::INFO << "Trigger channel " << i << ": " << trigData->getTrigChannel( i )
          << endmsg;
    m_timing = trigData->getTimingType();
    // cout<<"-----------------trigger timing type-----------------------: "<<trigtiming<<endl;
    log << MSG::INFO << "Tigger Timing type: " << trigtiming << endmsg;
  }
 
  //------------------------------------
  //  Initialize MdcxHits
  //------------------------------------
  m_mdcxHits.reset();
  uint32_t getDigiFlag = 0;
  getDigiFlag += m_maxMdcDigi;
  if ( m_dropHot || m_salvageTrk ) getDigiFlag |= MdcRawDataProvider::b_dropHot;
  if ( m_keepBadTdc ) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
  if ( m_keepUnmatch ) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
  mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec( getDigiFlag );
  t_nDigi    = mdcDigiVec.size();
 
  // fill Mc truth
  // if(m_hist) fillMcTruth();
 
  // skip event by hit numbe
  if ( t_nDigi < m_minMdcDigi )
  {
    if ( 0 == t_nDigi ) { log << MSG::WARNING << " No hits in MdcDigiVec" << endmsg; }
    log << MSG::WARNING << " Skip this event for MdcDigiVec.size() < " << m_minMdcDigi
        << endmsg;
    return StatusCode::SUCCESS;
  }
  m_mdcxHits.create( mdcDigiVec, m_bunchT0, m_cresol );
  const HepAList<MdcxHit>& dchitlist = m_mdcxHits.GetMdcxHitList();
 
  if ( m_debug > 2 ) m_mdcxHits.print( std::cout, 6796 );
 
  //--------------------------------------------
  //  Make segments (MdcxSeg's) out of MdcxHit's
  //--------------------------------------------
  MdcxFindSegs             dcsegs( dchitlist, m_debug );
  const HepAList<MdcxSeg>& seglist = dcsegs.GetMdcxSeglist();
  if ( m_debug > 1 ) { dumpMdcxSegs( seglist ); }
  t_nSeg = seglist.length();
  // if(m_hist){ fillMdcxSegs(seglist);}
 
#ifdef MDCXTIMEDEBUG
  m_timer[1]->stop();
  m_timer[2]->start();
#endif
  //--------------------------------------------
  //  Make tracks (MdcxFittedHel's) out of MdcxSeg's
  //--------------------------------------------
  MdcxFindTracks           dctrks( seglist, m_debug );
  HepAList<MdcxFittedHel>& firsttrkl = (HepAList<MdcxFittedHel>&)dctrks.GetMdcxTrklist();
  if ( m_debug > 1 ) dumpTrackList( firsttrkl );
 
#ifdef MDCXTIMEDEBUG
  m_timer[2]->stop();
  m_timer[3]->start();
#endif
  // if(m_hist){ fillTrkl(firsttrkl);}
 
  // if(m_debug>1){
  //   std::cout << "dchitlist after find tracks before MergeDups, nhits=" <<
  //   dchitlist.length() << std::endl; for (int ii = 0; ii < dchitlist.length(); ii++) {
  //     dchitlist[ii]->print(std::cout, ii);
  //   }
  //   std::cout<<std::endl;
  // }
  MdcxMergeDups            dcmergeem( firsttrkl, m_debug );
  HepAList<MdcxFittedHel>& trkl = (HepAList<MdcxFittedHel>&)dcmergeem.GetMergedTrklist();
 
  // if (m_debug > 1 ){
  //   cout << "MdcxTrackFinder: after MergeDups, have "
  //   << trkl.length() << " track(s). nhits=" << dchitlist.length() << endl;
  //   for (int ii = 0; ii < dchitlist.length(); ii++) {
  //     dchitlist[ii]->print(std::cout, ii);
  //   }
  //   std::cout<<std::endl;
  // }
 
  //---------------------------------------------------------
  //  Put my tracks into official fitter and store to TDS
  //----------------------------------------------------------
 
  sc = FitMdcxTrack( trkl, dchitlist, m_hitCol, trackList, hitList );
  if ( !sc.isSuccess() ) { return StatusCode::SUCCESS; }
  t_nTdsTk = trackList->size();
 
  t_nTkTot += trackList->size();
  if ( t_nTdsTk < 20 ) t_nTkNum[t_nTdsTk]++;
 
#ifdef MDCXTIMEDEBUG
  m_timer[0]->stop();
  m_timer[3]->stop();
#endif
 
  if ( m_hist ) fillEvent();
  if ( m_debug > 0 )
  {
    DataObject* pNode;
    eventSvc()->retrieveObject( "/Event/Recon/RecMdcTrackCol", pNode );
    RecMdcTrackCol* tmpTrackCol = dynamic_cast<RecMdcTrackCol*>( pNode );
    eventSvc()->retrieveObject( "/Event/Recon/RecMdcHitCol", pNode );
    int nTdsTk = 0;
    if ( tmpTrackCol ) nTdsTk = tmpTrackCol->size();
 
    // if (t_evtNo % 1000 == 0) {
    std::cout << "MdcxTrackFinder: evtNo " << m_eventNo << " t0=" << m_bunchT0 << " Found "
              << trkl.length() << " keep " << t_nTdsTk << " finialy keep " << nTdsTk;
 
    int ndelete = 0;
    trkl.length() - trackList->size();
    if ( ndelete > 0 ) std::cout << " delete " << ndelete;
    std::cout << " track(s)" << endl;
    //}
 
    if ( m_debug > 1 ) dumpTdsTrack( tmpTrackCol );
    if ( m_debug > 1 ) dumpTrack( tmpTrackCol );
    // dumpTdsHits(tmpHitCol);
  }
  if ( ( trackList->size() != 4 ) ) b_saveEvent = true;
  setFilterPassed( b_saveEvent );
 
  return StatusCode::SUCCESS;
}
 
StatusCode MdcxTrackFinder::finalize() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "in finalize()" << endmsg;
 
  std::cout << " --after " << name() << " keep " << t_nTkTot << " tracks  " << std::endl;
  for ( int i = 0; i < 20; i++ )
  {
    if ( t_nTkNum[i] > 0 ) std::cout << " nTk=" << i << " " << t_nTkNum[i] << std::endl;
  }
 
  // tot evtNo, trkNum
  return StatusCode::SUCCESS;
}
 
void MdcxTrackFinder::MdcxHitsToHots( MdcxHel& mdcxHelix, const HepAList<MdcxHit>& mdcxHits,
                                      TrkHitList* m_trkHitList, MdcHitCol* mdcHitCol ) {
 
  if ( 0 == mdcxHits.length() ) return;
 
  int ihits = 0;
  while ( mdcxHits[ihits] )
  {
    int           ambig  = 0; //=mdcxHelix.Doca_Samb();//yzhang delete 2011-10-13
    const MdcHit* newhit = mdcxHits[ihits]->getMdcHit();
    if ( 0 == newhit )
    {
      const MdcDigi* theDigi = mdcxHits[ihits]->getDigi();
      int            layer   = MdcID::layer( mdcxHits[ihits]->getDigi()->identify() );
      int            wire    = MdcID::wire( mdcxHits[ihits]->getDigi()->identify() );
      m_digiMap[layer][wire] = mdcxHits[ihits]->getDigi();
      MdcHit* thehit         = new MdcHit( theDigi, m_gm );
      thehit->setCalibSvc( m_mdcCalibFunSvc );
      thehit->setCountPropTime( m_countPropTime );
      // thehit->setCosmicFit(m_cosmicFit);
 
      mdcHitCol->push_back( thehit );
      newhit = thehit;
    }
    MdcRecoHitOnTrack temp( *newhit, ambig, m_bunchT0 ); // m_bunchT0 nano second here
    MdcHitOnTrack*    newhot = &temp;
    // double fltLen = mdcxHelix.Doca_FLen();
    // std::cout<<" fltLen-- "<<fltLen<<"   "<<std::endl;
    // newhot->setFltLen(fltLen);
    // newhot->print(std::cout);
    // std::cout<< __FILE__ << " ambig  " << ambig << " append  "<<std::endl;
 
    // Store MdcxHits to TrkHitList
    m_trkHitList->appendHot( newhot ); // yzhang TEMP FIXME
    newhot->setUsability( true );      // yzhang add 2011-10-12 TEMP
    newhot->setActivity( true );       // yzhang add 2011-10-12 TEMP
 
    ihits++;
  }
}
 
StatusCode MdcxTrackFinder::FitMdcxTrack( HepAList<MdcxFittedHel>& trkl,
                                          const HepAList<MdcxHit>& dchitlist,
                                          MdcHitCol* m_hitCol, RecMdcTrackCol* trackList,
                                          RecMdcHitCol* hitList ) {
  StatusCode sc;
  MsgStream  log( msgSvc(), name() );
 
  //--Add Hit to MdcxFittedHel
  // if(m_debug>1){
  //  std::cout << "dchitlist before addHits nhits=" << dchitlist.length() << std::endl;
  //  for (int ii = 0; ii < dchitlist.length(); ii++) {
  //    dchitlist[ii]->print(std::cout, ii);
  //  }
  //  std::cout<<std::endl;
  //}
  MdcxAddHits dcaddem( trkl, dchitlist, m_addHitCut );
  // if (m_debug > 1){
  //   cout << "MdcxTrackFinder: after AddHits, have "
  //     << trkl.length() << " track(s). nhits=" << dchitlist.length() << endl;
  //   for (int ii = 0; ii < dchitlist.length(); ii++) {
  //     dchitlist[ii]->print(std::cout, ii);
  //   }
  //   std::cout<<std::endl;
  // }
 
  TrkLineMaker  linefactory;
  TrkHelixMaker helixfactory;
 
  int tkLen = trkl.length(); // FIXME
  for ( int kk = 0; kk < tkLen; kk++ )
  {
    const HepAList<MdcxHit>& xhits = trkl[kk]->XHitList();
    if ( m_debug > 2 )
    {
      std::cout << __FILE__ << " FitMdcxTrack " << kk << std::endl;
      for ( int i = 0; i < xhits.length(); i++ ) { xhits[i]->print( std::cout ); }
      std::cout << std::endl;
    }
 
    if ( m_debug > 2 ) std::cout << " before add hits nhits=" << xhits.length() << std::endl;
    HepAList<MdcxHit> xass = dcaddem.GetAssociates( kk );
    if ( m_debug > 2 ) { std::cout << " after,add " << xass.length() << " hits" << std::endl; }
 
    MdcxFittedHel  mdcxHelix = *trkl[kk];
    double         thechisq  = mdcxHelix.Chisq();
    TrkExchangePar tt( -mdcxHelix.D0(), mdcxHelix.Phi0(), mdcxHelix.Omega(), -mdcxHelix.Z0(),
                       -mdcxHelix.Tanl() );
    TrkRecoTrk*    aTrack;
    int            nparm;
 
    if ( m_lineFit )
    {
      /// m_bunchT0 in "ns" here, but "second" in TrkLineMaker&TrkHelixMaker
      aTrack = linefactory.makeTrack( tt, thechisq, *m_context, m_bunchT0 * 1.e-9 );
      nparm  = 4;
      linefactory.setFlipAndDrop( *aTrack, true, true );
    }
    else
    {
      aTrack = helixfactory.makeTrack( tt, thechisq, *m_context, m_bunchT0 * 1.e-9 );
      nparm  = 5;
      helixfactory.setFlipAndDrop( *aTrack, true, true );
    }
 
    TrkHitList* m_trkHitList = aTrack->hits();
    if ( 0 == m_trkHitList )
    {
      delete aTrack;
      aTrack = NULL;
      continue;
    }
 
    MdcxHitsToHots( mdcxHelix, xhits, m_trkHitList, m_hitCol );
    // std::cout<<"xhits---------------------"<<std::endl;//debug
    // m_trkHitList->hotList().printAll(cout);//debug
    MdcxHitsToHots( mdcxHelix, xass, m_trkHitList, m_hitCol );
    // std::cout<<"xass----------------------"<<std::endl;//debug
    // m_trkHitList->hotList().printAll(cout);//debug
    // std::cout<<"size "<<m_trkHitList->hotList().nHit()<<std::endl;
    // int beforDrop = m_trkHitList->hotList().nHit();
    if ( m_dropMultiHotInLayer ) dropMultiHotInLayer( m_trkHitList ); // yzhang debug FIXME
    // int afterDrop = m_trkHitList->hotList().nHit();
    // std::cout<<"drop "<<beforDrop-afterDrop<<" keep:"<<afterDrop<<::endl;
    if ( m_debug > 1 ) { std::cout << "== put to official fitter  " << std::endl; }
    TrkErrCode err = m_trkHitList->fit();
 
    if ( m_debug > 1 )
    {
      std::cout << "== after official fitter  " << std::endl;
      aTrack->printAll( std::cout );
    }
 
    const TrkFit* theFit = aTrack->fitResult();
    float         rcs    = 10000.0;
 
    if ( theFit )
    {
      int ndof = theFit->nActive() - nparm;
      if ( ndof > 0 ) rcs = theFit->chisq() / float( ndof );
      if ( m_debug > 1 )
      {
        if ( 4 == nparm ) cout << " TrkLineMaker";
        else cout << " TrkHelixMaker";
        cout << " trkNo. " << kk << " success " << err.success() << " rcs " << rcs << " chi2 "
             << theFit->chisq() << " nactive  " << theFit->nActive() << endl;
      }
    }
 
    if ( ( 1 == err.success() ) && ( rcs < 20.0 ) )
    {
      if ( m_debug > 1 ) std::cout << "== fit success " << std::endl; // yzhang debug
      if ( 4 == nparm ) { linefactory.setFlipAndDrop( *aTrack, false, false ); }
      else { helixfactory.setFlipAndDrop( *aTrack, false, false ); }
      //-------------Stick the found tracks into the list in RecMdcTrackCol--------
      if ( m_debug > 1 ) { cout << "MdcxTrackFinder: accept a track " << endl; }
      // update history
      aTrack->status()->addHistory( err, name().c_str() ); // yzhang FIXME
      store( aTrack, trackList, hitList );                 // aTrack have been deleted
    }
    else if ( ( 2 == err.success() ) && ( rcs < 150.0 ) )
    { //////// zoujh
      if ( m_debug > 1 )
        std::cout << "== fit success = 2, refit now" << std::endl; // yzhang debug
      int nrefit = 0;
      while ( nrefit++ < 5 )
      {
        if ( m_debug > 1 ) std::cout << "refit time " << nrefit << std::endl;
        err = m_trkHitList->fit();
        if ( err.success() == 1 ) break;
      }
      if ( err.success() == 1 )
      {
        if ( 4 == nparm ) { linefactory.setFlipAndDrop( *aTrack, false, false ); }
        else { helixfactory.setFlipAndDrop( *aTrack, false, false ); }
        //-------------Stick the found tracks into the list in RecMdcTrackCol--------
        // if (m_debug>1) { cout << "MdcxTrackFinder: accept a track and store to TDS" << endl;
        // } update history
        aTrack->status()->addHistory( err, name().c_str() ); // yzhang FIXME
        store( aTrack, trackList, hitList );                 // aTrack have been deleted
      } //////////////////////////zoujh
    }
    else
    {
      if ( m_debug > 1 )
      {
        std::cout << "== fit faild " << std::endl; // yzhang debug
        err.print( cout );
        cout << endl;
      }
      delete aTrack;
      aTrack = NULL;
      //---------------------------------------
      //  Fit no good; try a better input helix
      //---------------------------------------
      if ( m_debug > 1 )
        std::cout << "== Fit no good; try a better input helix" << std::endl; // yzhang debug
      mdcxHelix.Grow( *trkl[kk], xass );
      mdcxHelix.VaryRes();
      mdcxHelix.SetChiDofBail( 1500 ); // yzhang add 2009-11-03
      int fail = mdcxHelix.ReFit();
      if ( m_debug > 1 ) std::cout << __FILE__ << " refit fail:" << fail << std::endl;
      if ( !mdcxHelix.Fail() )
      {
        const HepAList<MdcxHit>& bxhits = mdcxHelix.XHitList();
        thechisq                        = mdcxHelix.Chisq();
        TrkExchangePar tb( mdcxHelix.D0(), mdcxHelix.Phi0(), mdcxHelix.Omega(), mdcxHelix.Z0(),
                           mdcxHelix.Tanl() );
        TrkRecoTrk*    bTrack;
        if ( 4 == nparm )
        {
          bTrack = linefactory.makeTrack( tb, thechisq, *m_context, m_bunchT0 * 1.e-9 );
          linefactory.setFlipAndDrop( *bTrack, false, false );
        }
        else
        {
          bTrack = helixfactory.makeTrack( tb, thechisq, *m_context, m_bunchT0 * 1.e-9 );
          helixfactory.setFlipAndDrop( *bTrack, false, false );
        }
        TrkHitList* bhits = bTrack->hits();
        if ( 0 == bhits )
        {
          delete bTrack;
          bTrack = NULL;
          continue;
        }
 
        MdcxHitsToHots( mdcxHelix, bxhits, bhits, m_hitCol );
        TrkErrCode    berr = bhits->fit();
        const TrkFit* bFit = bTrack->fitResult();
        rcs                = 10000.0;
        if ( bFit )
        {
          int ndof = bFit->nActive() - nparm;
          if ( ndof > 0 ) rcs = bFit->chisq() / float( ndof );
          if ( m_debug > 1 )
          {
            if ( 4 == nparm ) cout << " TrkLineMaker";
            else cout << " TrkHelixMaker";
            cout << " success trkNo. " << kk << " status " << berr.success() << " rcs " << rcs
                 << " chi2 " << bFit->chisq() << " nactive  " << bFit->nActive() << endl;
          }
        }
        if ( ( 1 == berr.success() ) && ( rcs < 50.0 ) )
        {
          // update history
          bTrack->status()->addHistory( berr, name().c_str() ); // yzhang FIXME
          if ( m_debug > 1 )
          {
            cout << "MdcxTrackFinder: accept b track and store to TDS" << endl;
            bTrack->printAll( cout );
          }
          store( bTrack, trackList, hitList ); // bTrack have been deleted
        }
        else
        {
          if ( m_debug > 1 )
          {
            cout << " fit failed " << endl;
            berr.print( cout );
            cout << endl;
          }
          if ( bTrack )
          {
            delete bTrack;
            bTrack = NULL;
          }
        }
      }
      else
      {
        // cout<< " grow and refit failed "<<endl;
      }
    }
  }
  if ( m_debug > 1 ) dumpTdsTrack( trackList );
  return StatusCode::SUCCESS;
 
} // end of FitMdcxTrack
 
void MdcxTrackFinder::store( TrkRecoTrk* aTrack, RecMdcTrackCol* trackList,
                             RecMdcHitCol* hitList ) {
  assert( aTrack );
  nTk++;
  int            trackId = trackList->size();
  TrkExchangePar helix   = aTrack->fitResult()->helix( 0. );
  if ( m_dropTrkPt > 0. && ( aTrack->fitResult()->pt() < m_dropTrkPt ) )
  {
    if ( m_debug > 1 )
      std::cout << "MdcxTrackFinder delete track by pt " << aTrack->fitResult()->pt()
                << "<ptCut " << m_dropTrkPt << std::endl;
    return;
  }
 
  if ( ( ( fabs( helix.d0() ) > m_d0Cut ) || ( fabs( helix.z0() ) > m_z0Cut ) ) )
  {
    if ( m_debug > 1 )
      std::cout << name() << " delete track by d0 " << helix.d0() << ">d0Cut " << m_d0Cut
                << " or z0 " << helix.z0() << " >z0Cut " << m_z0Cut << std::endl;
    return;
  }
 
  if ( m_hist ) fillTrack( aTrack );
  MdcTrack mdcTrack( aTrack ); // aTrack have been deleted in ~MdcTrack()
  // tkStat:  0,PatRec 1,MdcxReco 2,Tsf 3,CurlFinder -1,Combined cosmic
  int tkStat     = 1;
  int nHitbefore = hitList->size();
 
  mdcTrack.storeTrack( trackId, trackList, hitList, tkStat );
  int nHitAfter = hitList->size();
  if ( nHitAfter - nHitbefore < 10 ) setFilterPassed( true );
}
 
void MdcxTrackFinder::dumpTrack( RecMdcTrackCol* trackList ) {
  RecMdcTrackCol::iterator i_tk = trackList->begin();
  for ( ; i_tk != trackList->end(); i_tk++ ) { printTrack( *( i_tk ) ); }
} // dumpTrack
 
void MdcxTrackFinder::printTrack( RecMdcTrack* tk ) {
  // yzhang debug
  std::cout << " MdcTrack Id:" << tk->trackId() << " q:" << tk->charge() << std::endl;
  std::cout << "dr   Fi0   Cpa   Dz   Tanl   Chi2   Ndf   nSt   FiTerm  poca" << std::endl;
  std::cout << "(" << setw( 5 ) << tk->helix( 0 ) << "," << setw( 5 ) << tk->helix( 1 ) << ","
            << setw( 5 ) << tk->helix( 2 ) << "," << setw( 5 ) << tk->helix( 3 ) << ","
            << setw( 5 ) << tk->helix( 4 ) << ")" << setw( 5 ) << tk->chi2() << setw( 4 )
            << tk->ndof() << setw( 4 ) << tk->getNhits() << setw( 4 ) << tk->nster()
            << setw( 5 ) << tk->getFiTerm() << tk->poca() << std::endl;
  std::cout << " ErrMat  " << tk->err() << std::endl;
 
  std::cout << "hitId  tkId  (l,w) fltLen  lr  dt  ddl  tdc " //<<chi2Add
            << "doca  entr  z  tprop stat  " << std::endl;
 
  HitRefVec           hl = tk->getVecHits();
  HitRefVec::iterator it = hl.begin();
  for ( ; it != hl.end(); ++it )
  {
    RecMdcHit*      h         = *it;
    int             layer     = MdcID::layer( h->getMdcId() );
    double          _vprop    = ( layer < 8 ) ? Constants::vpropInner : Constants::vpropOuter;
    const MdcLayer* _layerPtr = m_gm->Layer( layer );
    double          _zlen     = _layerPtr->zLength();
    double          z         = h->getZhit();
    double          tprop;
    if ( 0 == layer % 2 )
    {
      tprop = ( 0.5 * _zlen + z ) / _vprop; // odd
    }
    else
    {
      tprop = ( 0.5 * _zlen - z ) / _vprop; // even
    }
    // build the sense wires
    // const MdcSWire* wire =
    // m_gm->Wire(MdcID::layer(h->getMdcId()),MdcID::wire(h->getMdcId())); double z =
    // wire->zForw(); while(z<wire->zRear()){ HepPoint3D pos; Hep3Vector dir;
    // if(!(wire->getTraj()==NULL)){
    // wire->getTraj()->getInfo(z,pos,dir);
    //}
    // std::cout<<"("<< wire->layer()->layNum()<<","<<wire->cell()<<" "<<wire->Id()<<")";
    // std::cout<<" z, sag:"<<z
    //<<", "<<wire->getTraj()->deltaY(z-wire->zForw())<<std::endl;
    ////<<" pos:"<<pos <<" dir:"<<dir<<std::endl;
    // z+=1.;
    // }
    std::cout << setw( 3 ) << h->getId() << setw( 2 ) << h->getTrkId() << "("
              << MdcID::layer( h->getMdcId() ) << "," << MdcID::wire( h->getMdcId() ) << ")"
              << setw( 10 ) << h->getFltLen() << setw( 2 ) << h->getFlagLR() << setw( 10 )
              << h->getDriftT() << setw( 12 ) << h->getDriftDistLeft()
              <<                    // setw(8) << h->getErrDriftDistRight() <<
        setw( 8 ) << h->getTdc() << // setw(10) << h->getChisqAdd() <<
        setw( 10 ) << h->getDoca() << setw( 10 ) << h->getEntra() << setw( 10 ) << h->getZhit()
              << setw( 10 ) << tprop << setw( 2 ) << h->getStat() << std::endl;
  }
} // print track
 
void MdcxTrackFinder::bookNTuple() {
  MsgStream  log( msgSvc(), name() );
  StatusCode sc;
  if ( !m_hist ) return;
 
  g_poison      = histoSvc()->book( "poison", "poison", 43, 0, 42, 288, 0, 288 );
  g_csmax4      = histoSvc()->book( "csmax4", "csmax4", 100, 0, 100 );
  g_csmax3      = histoSvc()->book( "csmax3", "csmax3", 50000, 0, 500000 );
  g_omegag      = histoSvc()->book( "omegag", "omegag", 1000, 0., 1. );
  g_dPhiAU      = histoSvc()->book( "dPhiAU", "dPhiAU", 1000, 0., 3.5 );
  g_dPhiAU_0    = histoSvc()->book( "dPhiAU_0", "dPhiAU_0", 1000, 0., 3.5 );
  g_dPhiAU_1    = histoSvc()->book( "dPhiAU_1", "dPhiAU_1", 1000, 0., 3.5 );
  g_dPhiAU_5    = histoSvc()->book( "dPhiAU_5", "dPhiAU_5", 1000, 0., 3.5 );
  g_dPhiAU_7    = histoSvc()->book( "dPhiAU_7", "dPhiAU_7", 1000, 0., 3.5 );
  g_dPhiAV      = histoSvc()->book( "dPhiAV", "dPhiAV", 1000, 0., 3.5 );
  g_addSegPhi   = histoSvc()->book( "addSegPhi", "addSegPhi", 1000, 0., 3.5 );
  g_dPhiAV_0    = histoSvc()->book( "dPhiAV_0", "dPhiAV_0", 1000, 0., 3.5 );
  g_dPhiAV_1    = histoSvc()->book( "dPhiAV_1", "dPhiAV_1", 1000, 0., 3.5 );
  g_dPhiAV_6    = histoSvc()->book( "dPhiAV_6", "dPhiAV_6", 1000, 0., 3.5 );
  g_dPhiAV_8    = histoSvc()->book( "dPhiAV_8", "dPhiAV_8", 1000, 0., 3.5 );
  g_trkllmk     = histoSvc()->book( "trkllmk", "trkllmk", 43, 0., 42 );
  g_trklgood    = histoSvc()->book( "trklgood", "trklgood", 43, 0., 42 );
  g_trklcircle  = histoSvc()->book( "trklcircle", "trklcircle", 43, 0., 42 );
  g_trklhelix   = histoSvc()->book( "trklhelix", "trklhelix", 43, 0., 42 );
  g_trkldrop1   = histoSvc()->book( "trkldrop1", "trkldrop1", 43, 0., 42 );
  g_trkldrop2   = histoSvc()->book( "trkldrop2", "trkldrop2", 43, 0., 42 );
  g_trklappend1 = histoSvc()->book( "trklappend1", "trklappend1", 43, 0., 42 );
  g_trklappend2 = histoSvc()->book( "trklappend2", "trklappend2", 43, 0., 42 );
  g_trklappend3 = histoSvc()->book( "trklappend3", "trklappend3", 43, 0., 42 );
  // g_fitOmega= histoSvc()->book( "fitOmega", "fitOmega",200,0.,100.);
  g_trklfirstProb = histoSvc()->book( "trklfirstProb", "trklfirstProb", 200, 0., 2. );
  g_trkltemp      = histoSvc()->book( "trkltemp", "trkltemp", 200, -3.5, 3.5 );
 
  g_trklproca = histoSvc()->book( "trklproca", "trklproca", 100, 0., 5. );
  g_trkle     = histoSvc()->book( "trkle", "trkle", 200, 0., 0.025 );
  g_trkld     = histoSvc()->book( "trkld", "trkld", 200, -1.2, 1.2 );
  g_trklel    = histoSvc()->book( "trklel", "trklel", 200, 0., 0.025, 43, 0, 43 );
  g_trkldl    = histoSvc()->book( "trkldl", "trkldl", 200, -1.2, 1.2, 43, 0, 43 );
  g_trkldoca  = histoSvc()->book( "trkldoca", "trkldoca", 200, -1.2, 1.2 );
  g_trkllayer = histoSvc()->book( "trkllayer", "trkllayer", 43, 0, 43 );
  g_dropHitsSigma =
      histoSvc()->book( "dropHitsSigma", "dropHitsSigma", 43, 0, 43, 100, 0, 11 );
  g_addHitCut   = histoSvc()->book( "addHitCut", "addHitCut", 100, 0, 10 );
  g_addHitCut2d = histoSvc()->book( "addHitCut2d", "addHitCut2d", 43, 0, 43, 100, 0, 10 );
  // g_addSegPhiDiff = histoSvc()->book( "addSegPhiDiff", "addSegPhiDiff",100,-6.28,6.28);
 
  NTuplePtr nt1( ntupleSvc(), "MdcxReco/rec" );
  if ( nt1 ) { m_xtuple1 = nt1; }
  else
  {
    m_xtuple1 = ntupleSvc()->book( "MdcxReco/rec", CLID_ColumnWiseTuple,
                                   "MdcxReco reconsturction results" );
    if ( m_xtuple1 )
    {
      sc = m_xtuple1->addItem( "t0", m_xt0 );
      sc = m_xtuple1->addItem( "timing", m_xtiming );
      sc = m_xtuple1->addItem( "t0Stat", m_xt0Stat );
      sc = m_xtuple1->addItem( "t0Truth", m_xt0Truth );
 
      sc = m_xtuple1->addItem( "nSlay", m_xnSlay );
      sc = m_xtuple1->addItem( "p", m_xp );
      sc = m_xtuple1->addItem( "pt", m_xpt );
      sc = m_xtuple1->addItem( "pz", m_xpz );
      sc = m_xtuple1->addItem( "d0", m_xd0 );
      sc = m_xtuple1->addItem( "phi0", m_xphi0 );
      sc = m_xtuple1->addItem( "cpa", m_xcpa );
      sc = m_xtuple1->addItem( "z0", m_xz0 );
      sc = m_xtuple1->addItem( "tanl", m_xtanl );
      sc = m_xtuple1->addItem( "q", m_xq );
      sc = m_xtuple1->addItem( "pocax", m_xpocax );
      sc = m_xtuple1->addItem( "pocay", m_xpocay );
      sc = m_xtuple1->addItem( "pocaz", m_xpocaz );
 
      sc = m_xtuple1->addItem( "evtNo", m_xevtNo );
      sc = m_xtuple1->addItem( "nSt", m_xnSt );
      sc = m_xtuple1->addItem( "nAct", m_xnAct );
      sc = m_xtuple1->addItem( "nDof", m_xnDof );
      sc = m_xtuple1->addItem( "chi2", m_xchi2 );
      sc = m_xtuple1->addItem( "tkId", m_xtkId );
      sc = m_xtuple1->addItem( "nHit", m_xnHit, 0, 7000 ); // # of hit/rec track
      sc = m_xtuple1->addIndexedItem( "resid", m_xnHit, m_xresid );
      sc = m_xtuple1->addIndexedItem( "sigma", m_xnHit, m_xsigma );
      sc = m_xtuple1->addIndexedItem( "driftD", m_xnHit, m_xdriftD );
      sc = m_xtuple1->addIndexedItem( "driftT", m_xnHit, m_xdriftT );
      sc = m_xtuple1->addIndexedItem( "doca", m_xnHit, m_xdoca );
      sc = m_xtuple1->addIndexedItem( "entra", m_xnHit, m_xentra );
      sc = m_xtuple1->addIndexedItem( "ambig", m_xnHit, m_xambig );
      sc = m_xtuple1->addIndexedItem( "fltLen", m_xnHit, m_xfltLen );
      sc = m_xtuple1->addIndexedItem( "tof", m_xnHit, m_xtof );
      sc = m_xtuple1->addIndexedItem( "act", m_xnHit, m_xact );
      sc = m_xtuple1->addIndexedItem( "tdc", m_xnHit, m_xtdc );
      sc = m_xtuple1->addIndexedItem( "adc", m_xnHit, m_xadc );
      sc = m_xtuple1->addIndexedItem( "layer", m_xnHit, m_xlayer );
      sc = m_xtuple1->addIndexedItem( "wire", m_xnHit, m_xwire );
      sc = m_xtuple1->addIndexedItem( "x", m_xnHit, m_xx );
      sc = m_xtuple1->addIndexedItem( "y", m_xnHit, m_xy );
      sc = m_xtuple1->addIndexedItem( "z", m_xnHit, m_xz );
    }
    else
    { // did not manage to book the N tuple....
      log << MSG::ERROR << "    Cannot book N-tuple: MdcxReco/rec" << endmsg;
      // return;
    }
  } // end book of nt1
  NTuplePtr nt4( ntupleSvc(), "MdcxReco/evt" );
  if ( nt4 ) { m_xtupleEvt = nt4; }
  else
  {
    m_xtupleEvt =
        ntupleSvc()->book( "MdcxReco/evt", CLID_ColumnWiseTuple, "MdcxReco event data" );
    if ( m_xtupleEvt )
    {
      sc = m_xtupleEvt->addItem( "evtNo", m_xt4EvtNo );
      sc = m_xtupleEvt->addItem( "nRecTk", m_xt4nRecTk );
      sc = m_xtupleEvt->addItem( "nTdsTk", m_xt4nTdsTk );
      sc = m_xtupleEvt->addItem( "t0", m_xt4t0 );
      sc = m_xtupleEvt->addItem( "t0Stat", m_xt4t0Stat );
      sc = m_xtupleEvt->addItem( "t0Truth", m_xt4t0Truth );
      sc = m_xtupleEvt->addItem( "time", m_xt4time );
      sc = m_xtupleEvt->addItem( "timeSeg", m_xt4timeSeg );
      sc = m_xtupleEvt->addItem( "timeTrack", m_xt4timeTrack );
      sc = m_xtupleEvt->addItem( "timeFit", m_xt4timeFit );
      sc = m_xtupleEvt->addItem( "nSeg", m_xt4nSeg );
      sc = m_xtupleEvt->addItem( "nDigi", m_xt4nDigi, 0, 7000 ); // # of hit/mc track
      sc = m_xtupleEvt->addIndexedItem( "layer", m_xt4nDigi, m_xt4Layer );
      sc = m_xtupleEvt->addIndexedItem( "rt", m_xt4nDigi, m_xt4Time );
      sc = m_xtupleEvt->addIndexedItem( "rc", m_xt4nDigi, m_xt4Charge );
      // sc = m_xtupleEvt->addIndexedItem ("rawHit", m_xt4nDigi, m_xt4rawHit);
      // sc = m_xtupleEvt->addIndexedItem ("recHit", m_xt4nDigi, m_xt4recHit);
    }
    else
    { // did not manage to book the N tuple....
      log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/evt" << endmsg;
      // return;
    }
  } // end of book nt4
 
  // book tuple of segment
  NTuplePtr ntSeg( ntupleSvc(), "MdcxReco/seg" );
  if ( ntSeg ) { m_xtupleSeg = ntSeg; }
  else
  {
    m_xtupleSeg = ntupleSvc()->book( "MdcxReco/seg", CLID_ColumnWiseTuple,
                                     "MdcxTrackFinder segment data" );
    if ( m_xtupleSeg )
    {
      sc = m_xtupleSeg->addItem( "sl", m_xtsSl );
      sc = m_xtupleSeg->addItem( "d0", m_xtsD0 );
      sc = m_xtupleSeg->addItem( "omega", m_xtsOmega );
      sc = m_xtupleSeg->addItem( "phi0", m_xtsPhi0 );
      sc = m_xtupleSeg->addItem( "d0sl", m_xtsD0_sl_approx );
      sc = m_xtupleSeg->addItem( "phi0sl", m_xtsPhi0_sl_approx );
      sc = m_xtupleSeg->addItem( "xbbrrf", m_xtsXline_bbrrf );
      sc = m_xtupleSeg->addItem( "ybbrrf", m_xtsYline_bbrrf );
      sc = m_xtupleSeg->addItem( "xslope", m_xtsXline_slope );
      sc = m_xtupleSeg->addItem( "yslope", m_xtsYline_slope );
      sc = m_xtupleSeg->addItem( "pat", m_xtsPat );
      sc = m_xtupleSeg->addItem( "chisq", m_xtsChisq );
      sc = m_xtupleSeg->addItem( "nDigi", m_xtsNDigi, 0, 20 );
      sc = m_xtupleSeg->addIndexedItem( "layer", m_xtsNDigi, m_xtsLayer );
      sc = m_xtupleSeg->addIndexedItem( "wire", m_xtsNDigi, m_xtsWire );
      sc = m_xtupleSeg->addIndexedItem( "inSeg", m_xtsNDigi, m_xtsInSeg );
    }
    else
    {
      log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/evt" << endmsg;
      // return;
    }
  }
  NTuplePtr nt5( ntupleSvc(), "MdcxReco/trkl" );
  if ( nt5 ) { m_xtupleTrkl = nt5; }
  else
  {
    m_xtupleTrkl =
        ntupleSvc()->book( "MdcxReco/trkl", CLID_RowWiseTuple, "MdcxReco track info" );
    if ( m_xtupleTrkl )
    {
      sc = m_xtupleTrkl->addItem( "layer", m_xt5Layer );
      sc = m_xtupleTrkl->addItem( "wire", m_xt5Wire );
    }
    else
    {
      log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/trkl" << endmsg;
      // return;
    }
  }
 
  NTuplePtr ntCsmcSew( ntupleSvc(), "MdcxReco/csmc" );
  if ( ntCsmcSew ) { m_xtupleCsmcSew = ntCsmcSew; }
  else
  {
    m_xtupleCsmcSew = ntupleSvc()->book( "MdcxReco/csmc", CLID_ColumnWiseTuple,
                                         "MdcxReco reconsturction results" );
    if ( m_xtupleCsmcSew )
    {
      sc = m_xtupleCsmcSew->addItem( "dD0", m_csmcD0 );
      sc = m_xtupleCsmcSew->addItem( "dPhi0", m_csmcPhi0 );
      sc = m_xtupleCsmcSew->addItem( "dZ0", m_csmcZ0 );
      sc = m_xtupleCsmcSew->addItem( "dOmega", m_csmcOmega );
      sc = m_xtupleCsmcSew->addItem( "dPt", m_csmcPt );
      sc = m_xtupleCsmcSew->addItem( "dTanl", m_csmcTanl );
    }
    else
    {
      log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/csmc" << endmsg;
      // return;
    }
  }
  NTuplePtr ntSegAdd1( ntupleSvc(), "MdcxReco/addSeg1" );
  if ( ntSegAdd1 ) { m_xtupleAddSeg1 = ntSegAdd1; }
  else
  {
    m_xtupleAddSeg1 =
        ntupleSvc()->book( "MdcxReco/addSeg1", CLID_ColumnWiseTuple, "MdcxReco event data" );
    if ( m_xtupleAddSeg1 )
    {
      sc = m_xtupleAddSeg1->addItem( "same", m_addSegSame );
      sc = m_xtupleAddSeg1->addItem( "seedSl", m_addSegSeedSl );
      sc = m_xtupleAddSeg1->addItem( "seedPhi", m_addSegSeedPhi );
      sc = m_xtupleAddSeg1->addItem( "seedPhiLay", m_addSegSeedPhiLay );
      sc = m_xtupleAddSeg1->addItem( "seedLen", m_addSegSeedLen );
      sc = m_xtupleAddSeg1->addItem( "seedD0", m_addSegSeedD0 );
      sc = m_xtupleAddSeg1->addItem( "seedPhi0", m_addSegSeedPhi0 );
      sc = m_xtupleAddSeg1->addItem( "addSl", m_addSegAddSl );
      sc = m_xtupleAddSeg1->addItem( "addPhi", m_addSegAddPhi );
      sc = m_xtupleAddSeg1->addItem( "addPhiLay", m_addSegAddPhiLay );
      sc = m_xtupleAddSeg1->addItem( "addLen", m_addSegAddLen );
      sc = m_xtupleAddSeg1->addItem( "addD0", m_addSegAddD0 );
      sc = m_xtupleAddSeg1->addItem( "addPhi0", m_addSegAddPhi0 );
    }
    else
    {
      log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/addSeg1" << endmsg;
      // return;
    }
  }
 
  NTuplePtr ntSegComb( ntupleSvc(), "MdcxReco/segComb" );
  if ( ntSegComb ) { m_xtupleSegComb = ntSegComb; }
  else
  {
    m_xtupleSegComb =
        ntupleSvc()->book( "MdcxReco/segComb", CLID_ColumnWiseTuple, "MdcxReco event data" );
    if ( m_xtupleSegComb )
    {
      sc = m_xtupleSegComb->addItem( "evtNo", m_segCombEvtNo );
      sc = m_xtupleSegComb->addItem( "omega", m_segCombOmega );
      sc = m_xtupleSegComb->addItem( "sameAU", m_segCombSameAU );
      sc = m_xtupleSegComb->addItem( "sameUV", m_segCombSameUV );
      sc = m_xtupleSegComb->addItem( "dlenAU", m_segCombDLenAU );
      sc = m_xtupleSegComb->addItem( "dlenUV", m_segCombDLenUV );
      sc = m_xtupleSegComb->addItem( "slA", m_segCombSlA );
      sc = m_xtupleSegComb->addItem( "slU", m_segCombSlU );
      sc = m_xtupleSegComb->addItem( "slV", m_segCombSlV );
      sc = m_xtupleSegComb->addItem( "phiA", m_segCombPhiA );
      sc = m_xtupleSegComb->addItem( "phiU", m_segCombPhiU );
      sc = m_xtupleSegComb->addItem( "phiV", m_segCombPhiV );
    }
    else
    {
      log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/segComb" << endmsg;
      // return;
    }
  }
 
  NTuplePtr ntDropHits( ntupleSvc(), "MdcxReco/dropHits" );
  if ( ntDropHits ) { m_xtupleDropHits = ntDropHits; }
  else
  {
    m_xtupleDropHits =
        ntupleSvc()->book( "MdcxReco/dropHits", CLID_ColumnWiseTuple, "MdcxReco event data" );
    if ( m_xtupleDropHits )
    {
      sc = m_xtupleDropHits->addItem( "evtNo", m_segDropHitsEvtNo );
      sc = m_xtupleDropHits->addItem( "layer", m_segDropHitsLayer );
      sc = m_xtupleDropHits->addItem( "wire", m_segDropHitsWire );
      sc = m_xtupleDropHits->addItem( "pull", m_segDropHitsPull );
      sc = m_xtupleDropHits->addItem( "doca", m_segDropHitsDoca );
      sc = m_xtupleDropHits->addItem( "sigma", m_segDropHitsSigma );
      sc = m_xtupleDropHits->addItem( "drift", m_segDropHitsDrift );
      sc = m_xtupleDropHits->addItem( "mcTkId", m_segDropHitsMcTkId );
    }
    else
    {
      log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/dropHits" << endmsg;
      // return;
    }
  }
  NTuplePtr ntAddSeg2( ntupleSvc(), "MdcxReco/addSeg2" );
  if ( ntAddSeg2 ) { m_xtupleAddSeg2 = ntAddSeg2; }
  else
  {
    m_xtupleAddSeg2 =
        ntupleSvc()->book( "MdcxReco/addSeg2", CLID_ColumnWiseTuple, "MdcxReco event data" );
    if ( m_xtupleAddSeg2 )
    {
      sc = m_xtupleAddSeg2->addItem( "evtNo", m_addSegEvtNo );
      sc = m_xtupleAddSeg2->addItem( "slayer", m_addSegSlayer );
      sc = m_xtupleAddSeg2->addItem( "poca", m_addSegPoca );
      sc = m_xtupleAddSeg2->addItem( "len", m_addSegLen );
    }
    else
    {
      log << MSG::ERROR << "Cannot book N-tuple: MdcxReco/addSeg2" << endmsg;
      // return;
    }
  }
  //--------------end of book ntuple------------------
}
 
void MdcxTrackFinder::fillEvent() {
  //-----------get raw digi-----------------------
  if ( m_xtupleEvt == NULL ) return;
#ifdef MDCXTIMEDEBUG
  m_xt4time      = m_timer[0]->elapsed();
  m_xt4timeSeg   = m_timer[1]->elapsed();
  m_xt4timeTrack = m_timer[2]->elapsed();
  m_xt4timeFit   = m_timer[3]->elapsed();
#endif
  m_xt4EvtNo                 = m_eventNo;
  m_xt4t0                    = m_bunchT0;
  m_xt4t0Stat                = m_t0Stat;
  m_xt4t0Truth               = m_t0Truth;
  m_xt4nRecTk                = nTk;
  m_xt4nTdsTk                = t_nTdsTk;
  m_xt4nDigi                 = t_nDigi;
  m_xt4nSeg                  = t_nSeg;
  int                  iDigi = 0;
  MdcDigiCol::iterator iter  = mdcDigiVec.begin();
  for ( ; iDigi < t_nDigi; iter++ )
  {
    int l             = MdcID::layer( ( *iter )->identify() );
    m_xt4Layer[iDigi] = l;
    // int w = MdcID::wire((*iter)->identify());
    m_xt4Time[iDigi]   = RawDataUtil::MdcTime( ( *iter )->getTimeChannel() );
    m_xt4Charge[iDigi] = RawDataUtil::MdcCharge( ( *iter )->getChargeChannel() );
    // m_xt4rawHit[l]++;
    iDigi++;
  } // end for iter
  m_xtupleEvt->write();
}
 
void MdcxTrackFinder::fillTrack( TrkRecoTrk* atrack ) {
 
  if ( !m_xtuple1 ) return;
  //-----------get MC truth data-------------------
  m_xevtNo               = m_eventNo;
  int recHitMap[43][288] = { 0 };
  // int nTk = (*m_tracks).nTrack();
  // for (int itrack = 0; itrack < nTk; itrack++) {
  if ( atrack == NULL ) return;
 
  const TrkFit* fitResult = atrack->fitResult();
  if ( fitResult == 0 ) return; // check the fit worked
 
  // for fill ntuples
  int nSt     = 0; // int nSeg=0;
  int seg[11] = { 0 };
  int segme;
  //-----------hit list-------------
  HitRefVec         hitRefVec;
  const TrkHitList* hitlist = atrack->hits();
 
  TrkHitList::hot_iterator hot            = hitlist->begin();
  int                      layerCount[43] = { 0 };
  int                      i              = 0;
  for ( ; hot != hitlist->end(); hot++ )
  {
 
    const MdcRecoHitOnTrack* recoHot;
    recoHot     = dynamic_cast<const MdcRecoHitOnTrack*>( &( *hot ) );
    int layer   = recoHot->mdcHit()->layernumber();
    int wire    = recoHot->mdcHit()->wirenumber();
    m_xlayer[i] = layer;
    // m_xt4recHit[layer]++;//fill rec hit for hit eff
    m_xwire[i] = wire;
    layerCount[layer]++;
    recHitMap[layer][wire]++;
    m_xambig[i] = recoHot->wireAmbig(); // wire ambig
    // fltLen
    double fltLen = recoHot->fltLen();
    m_xfltLen[i]  = fltLen;
    double tof    = recoHot->getParentRep()->arrivalTime( fltLen );
    // position
    HepPoint3D pos;
    Hep3Vector dir;
    recoHot->getParentRep()->traj().getInfo( fltLen, pos, dir );
    m_xx[i]      = pos.x();
    m_xy[i]      = pos.y();
    m_xz[i]      = pos.z();
    m_xdriftT[i] = recoHot->mdcHit()->driftTime( tof, pos.z() );
    m_xtof[i]    = tof;
    m_xdriftD[i] = recoHot->drift();
    m_xsigma[i]  = recoHot->hitRms();
    m_xtdc[i]    = recoHot->rawTime();
    m_xadc[i]    = recoHot->mdcHit()->charge();
    m_xdoca[i]   = recoHot->dcaToWire(); // sign w.r.t. dirft() FIXME
    m_xentra[i]  = recoHot->entranceAngle();
    // m_xentraHit[i]    = recoHot->entranceAngleHit();
    m_xact[i] = recoHot->isActive();
    // resid
    double res = 999., rese = 999.;
    if ( recoHot->resid( res, rese, false ) ) {}
    else {}
    m_xresid[i] = res;
    // for n seg
    segme = 0;
    if ( layer > 0 ) { segme = ( layer - 1 ) / 4; }
    seg[segme]++;
    if ( recoHot->layer()->view() ) { ++nSt; }
    i++;
  } // end fill hit
 
  int nSlay = 0;
  for ( int i = 0; i < 11; i++ )
  {
    if ( seg[i] > 0 ) nSlay++;
  }
  m_xnSlay = nSlay;
 
  //------------fill track result-------------
  // m_xtkId = itrack;
  // track parameters at closest approach to beamline
  double         fltLenPoca = 0.0;
  TrkExchangePar helix      = fitResult->helix( fltLenPoca );
  m_xphi0                   = helix.phi0();
  m_xtanl                   = helix.tanDip();
  m_xz0                     = helix.z0();
  m_xd0                     = helix.d0();
  if ( fabs( m_xz0 ) > 20 || fabs( m_xd0 ) > 2 )
  {
    // b_saveEvent = true;
    if ( m_debug > 1 )
      std::cout << "evt:" << m_eventNo << " BigVtx:"
                << " d0 " << helix.d0() << " z0 " << helix.z0() << std::endl;
  }
  m_xpt = fitResult->pt();
  if ( m_xpt > 0.00001 ) { m_xcpa = fitResult->charge() / fitResult->pt(); }
  // momenta and position
  CLHEP::Hep3Vector p1 = fitResult->momentum( fltLenPoca );
  double            px, py, pz, pxy;
  pxy             = fitResult->pt();
  px              = p1.x();
  py              = p1.y();
  pz              = p1.z();
  m_xp            = p1.mag();
  m_xpz           = pz;
  HepPoint3D poca = fitResult->position( fltLenPoca );
  m_xpocax        = poca.x();
  m_xpocay        = poca.y();
  m_xpocaz        = poca.z();
 
  m_xq    = fitResult->charge();
  m_xnAct = fitResult->nActive();
  m_xnHit = hitlist->nHit();
  m_xnDof = fitResult->nDof();
  m_xnSt  = nSt;
  m_xchi2 = fitResult->chisq();
  // for (int l=0;l<43;l++) m_xlayerCount[l] = layerCount[l];
  m_xt0      = m_bunchT0;
  m_xtiming  = m_timing;
  m_xt0Stat  = m_t0Stat;
  m_xt0Truth = m_t0Truth;
  m_xtuple1->write();
  //}//end of loop rec trk list
 
} // end of
 
void MdcxTrackFinder::dumpMdcxSegs( const HepAList<MdcxSeg>& segList ) const {
  cout << name() << " found " << segList.length() << " segs :" << endl;
  for ( int i = 0; i < segList.length(); i++ )
  {
    std::cout << i << " ";
    segList[i]->printSeg();
  }
}
 
void MdcxTrackFinder::fillMdcxSegs( const HepAList<MdcxSeg>& segList ) const {
  if ( !m_xtupleSeg ) return;
 
  int cellMax[43] = { 40,  44,  48,  56,  64,  72,  80,  80,  76,  76,  88,
                      88,  100, 100, 112, 112, 128, 128, 140, 140, 160, 160,
                      160, 160, 176, 176, 176, 176, 208, 208, 208, 208, 240,
                      240, 240, 240, 256, 256, 256, 256, 288, 288, 288 };
  // Fill hits of every layer after segment finding
  int hitInSegList[43][288];
  for ( int ii = 0; ii < 43; ii++ )
  {
    for ( int jj = 0; jj < cellMax[ii]; jj++ ) { hitInSegList[ii][jj] = 0; }
  }
  for ( int i = 0; i < segList.length(); i++ )
  {
    MdcxSeg* aSeg = segList[i];
    for ( int ihit = 0; ihit < aSeg->Nhits(); ihit++ )
    {
      const MdcxHit* hit   = aSeg->XHitList()[ihit];
      int            layer = hit->Layer();
      int            wire  = hit->WireNo();
      hitInSegList[layer][wire]++;
    }
  }
  for ( int i = 0; i < segList.length(); i++ )
  {
    MdcxSeg* aSeg = segList[i];
    if ( aSeg == NULL ) continue;
    m_xtsSl             = aSeg->SuperLayer();
    m_xtsD0             = aSeg->D0();
    m_xtsOmega          = aSeg->Omega();
    m_xtsPhi0           = aSeg->Phi0();
    m_xtsD0_sl_approx   = aSeg->D0_sl_approx();
    m_xtsPhi0_sl_approx = aSeg->Phi0_sl_approx();
    m_xtsXline_bbrrf    = aSeg->Xline_bbrrf();
    m_xtsYline_bbrrf    = aSeg->Yline_bbrrf();
    m_xtsXline_slope    = aSeg->Xline_slope();
    m_xtsYline_slope    = aSeg->Yline_slope();
    int patIndex        = -1;
    for ( int ii = 0; ii < 14; ii++ )
    {
      if ( aSeg->Pat() == (int)MdcxSegPatterns::patt4[ii] )
      {
        patIndex = ii;
        break;
      }
    }
    for ( int ii = 0; ii < 20; ii++ )
    {
      if ( aSeg->Pat() == (int)MdcxSegPatterns::patt3[ii] )
      {
        patIndex = ii + 15;
        break;
      }
    }
    m_xtsPat   = patIndex;
    m_xtsChisq = aSeg->Chisq();
    m_xtsNDigi = aSeg->Nhits();
    for ( int ihit = 0; ihit < aSeg->Nhits(); ihit++ )
    {
      const MdcxHit* hit   = aSeg->XHitList()[ihit];
      int            layer = hit->Layer();
      int            wire  = hit->WireNo();
      m_xtsLayer[ihit]     = layer;
      m_xtsWire[ihit]      = wire;
      m_xtsInSeg[ihit]     = hitInSegList[layer][wire];
      // std::cout<< "hitInSegList   "<<hitInSegList[layer][wire] << std::endl;//yzhang debug
    }
    m_xtupleSeg->write();
  }
 
} // end of fillMdcxSegs
 
void MdcxTrackFinder::dumpTrackList( const HepAList<MdcxFittedHel>& trackList ) const {
  std::cout << "dump track list   " << std::endl;
  for ( int i = 0; i < trackList.length(); i++ )
  {
    std::cout << "track  " << i << std::endl;
    for ( int ihit = 0; ihit < trackList[i]->Nhits(); ihit++ )
    {
      const MdcxHit* hit   = trackList[i]->XHitList()[ihit];
      int            layer = hit->Layer();
      int            wire  = hit->WireNo();
      std::cout << " (" << layer << "," << wire << ") ";
    }
    std::cout << std::endl;
  }
} // end of dumpTrackList
 
void MdcxTrackFinder::fillTrkl( const HepAList<MdcxFittedHel>& firsttrkl ) const {
  if ( !m_xtupleTrkl ) return;
  int nDigi = 0;
  int iDigi = 0;
  for ( int i = 0; i < firsttrkl.length(); i++ ) { nDigi += firsttrkl[i]->Nhits(); }
  for ( int i = 0; i < firsttrkl.length(); i++ )
  {
    for ( int ihit = 0; ihit < firsttrkl[i]->Nhits(); ihit++ )
    {
      const MdcxHit* hit   = firsttrkl[i]->XHitList()[ihit];
      int            layer = hit->Layer();
      int            wire  = hit->WireNo();
      m_xt5Layer           = layer;
      m_xt5Wire            = wire;
      m_xtupleTrkl->write();
      iDigi++;
    }
  }
} // end of fillTrkl
 
void MdcxTrackFinder::fillMcTruth() {
  MsgStream  log( msgSvc(), name() );
  StatusCode sc;
  // Initialize
  for ( int ii = 0; ii < 43; ii++ )
  {
    for ( int jj = 0; jj < 288; jj++ ) { haveDigi[ii][jj] = -2; }
  }
  int nDigi = mdcDigiVec.size();
  for ( int iDigi = 0; iDigi < nDigi; iDigi++ )
  {
    int l = MdcID::layer( ( mdcDigiVec[iDigi] )->identify() );
    int w = MdcID::wire( ( mdcDigiVec[iDigi] )->identify() );
    // haveDigi[l][w]=(mdcDigiVec[iDigi])->getTrackIndex();
    haveDigi[l][w] = 1;
  }
 
  if ( m_mcHist )
  {
    //------------------get event start time truth-----------
    m_t0Truth = -10.;
    SmartDataPtr<Event::McParticleCol> mcParticleCol( eventSvc(), "/Event/MC/McParticleCol" );
    if ( !mcParticleCol ) { log << MSG::INFO << "Could not retrieve McParticelCol" << endmsg; }
    else
    {
      Event::McParticleCol::iterator iter_mc = mcParticleCol->begin();
      for ( ; iter_mc != mcParticleCol->end(); iter_mc++ )
      {
        if ( ( *iter_mc )->primaryParticle() )
        {
          m_t0Truth = ( *iter_mc )->initialPosition().t();
          // px = (*iter_mc)->initialFourMomentum().x()/1000.;//GeV
          // py = (*iter_mc)->initialFourMomentum().y()/1000.;//GeV
          // pz = (*iter_mc)->initialFourMomentum().z()/1000.;//GeV
        }
      }
    }
  }
  //------------------Retrieve MC truth MdcMcHit------------
  /*
     SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
     if (!mcMdcMcHitCol) {
     log << MSG::WARNING << "Could not find MdcMcHit" << endmsg;
     }else{
     Event::MdcMcHitCol::iterator iter_mchit = mcMdcMcHitCol->begin();
     for (;iter_mchit != mcMdcMcHitCol->end(); iter_mchit++ ) {
     const Identifier id= (*iter_mchit)->identify();
     int layer = MdcID::layer(id);
     int wire = MdcID::wire(id);
     mcDrift[layer][wire] = (*iter_mchit)->getDriftDistance();  //drift in MC.
     mcLR[layer][wire] = (*iter_mchit)->getPositionFlag();
     mcX[layer][wire] = (*iter_mchit)->getPositionX();
     mcY[layer][wire] = (*iter_mchit)->getPositionY();
     mcZ[layer][wire] = (*iter_mchit)->getPositionZ();
     if (mcLR[layer][wire] == 0) mcLR[layer][wire] = -1;
     }
     }
     */
}
 
void MdcxTrackFinder::dropMultiHotInLayer( TrkHitList* trkHitList ) {
  double tdr[43];
  double tdr_wire[43];
  for ( int i = 0; i < 43; i++ ) { tdr[i] = 9999.; }
 
  // make flag
  TrkHotList::hot_iterator hotIter = trkHitList->hotList().begin();
  while ( hotIter != trkHitList->hotList().end() )
  {
    MdcHitOnTrack* hot = const_cast<MdcHitOnTrack*>( &( *hotIter->mdcHitOnTrack() ) );
 
    // driftTime(tof,z)
    double dt    = hot->mdcHit()->driftTime( 0., 0. );
    int    layer = hot->mdcHit()->layernumber();
    int    wire  = hot->mdcHit()->wirenumber();
 
    // std::cout<<__FILE__<<" "<<dt<<  std::endl;
    if ( dt < tdr[layer] )
    {
      tdr[layer]      = dt;
      tdr_wire[layer] = wire;
    }
    hotIter++;
  }
 
  // std::cout<<" tdr wire ";
  // for(int i=0;i<43;i++){
  // std::cout<<i<<","<<tdr_wire[i]<<" tdr="<<tdr[i]<<"\n";
  // }
  //  inactive multi hit
  hotIter = trkHitList->hotList().begin();
  while ( hotIter != trkHitList->hotList().end() )
  {
    int layer = hotIter->mdcHitOnTrack()->mdcHit()->layernumber();
    int wire  = hotIter->mdcHitOnTrack()->mdcHit()->wirenumber();
    // double dt = hotIter->mdcHitOnTrack()->mdcHit()->driftTime(0.,0.)-m_bunchT0;
 
    if ( ( tdr[layer] < 9998. ) && ( tdr_wire[layer] != wire ) )
    {
      MdcHitOnTrack* hot = const_cast<MdcHitOnTrack*>( &( *hotIter->mdcHitOnTrack() ) );
      hot->setActivity( false );
      trkHitList->removeHit( hotIter->mdcHitOnTrack()->mdcHit() );
      // std::cout<<__FILE__<<" inactive "<< layer<<" "<<wire<<" dt "<<dt <<  std::endl;
    }
    else { hotIter++; }
  }
}
void MdcxTrackFinder::dumpTdsTrack( RecMdcTrackCol* trackList ) {
  std::cout << "tksize = " << trackList->size() << std::endl; // yzhang debug
  RecMdcTrackCol::iterator it = trackList->begin();
  for ( ; it != trackList->end(); it++ )
  {
    RecMdcTrack* tk = *it;
    std::cout << "//====RecMdcTrack " << tk->trackId() << "====:" << std::endl;
    cout << " d0 " << tk->helix( 0 ) << " phi0 " << tk->helix( 1 ) << " cpa " << tk->helix( 2 )
         << " z0 " << tk->helix( 3 ) << " tanl " << tk->helix( 4 ) << endl;
    std::cout << " q " << tk->charge() << " theta " << tk->theta() << " phi " << tk->phi()
              << " x0 " << tk->x() << " y0 " << tk->y() << " z0 " << tk->z() << " r0 "
              << tk->r() << endl;
    std::cout << " p " << tk->p() << " pt " << tk->pxy() << " px " << tk->px() << " py "
              << tk->py() << " pz " << tk->pz() << endl;
    std::cout << " tkStat " << tk->stat() << " chi2 " << tk->chi2() << " ndof " << tk->ndof()
              << " nhit " << tk->getNhits() << " nst " << tk->nster() << endl;
 
    int nhits = tk->getVecHits().size();
    std::cout << nhits << " Hits: " << std::endl;
    for ( int ii = 0; ii < nhits; ii++ )
    {
      Identifier id( tk->getVecHits()[ii]->getMdcId() );
      int        layer = MdcID::layer( id );
      int        wire  = MdcID::wire( id );
      cout << "(" << layer << "," << wire << "," << tk->getVecHits()[ii]->getStat()
           << ",lr:" << tk->getVecHits()[ii]->getFlagLR() << ") ";
    } // end of hit list
    std::cout << "  " << std::endl;
  } // end of tk list
  std::cout << "  " << std::endl;
 
} // end of dumpTdsTrack
 
void MdcxTrackFinder::dumpTdsHits( RecMdcHitCol* hitList ) {
  std::cout << __FILE__ << " " << __LINE__ << " All hits in TDS, nhit=" << hitList->size()
            << std::endl;
  RecMdcHitCol::iterator it = hitList->begin();
  for ( ; it != hitList->end(); it++ )
  {
    RecMdcHit* h = ( *it );
    Identifier id( h->getMdcId() );
    int        layer = MdcID::layer( id );
    int        wire  = MdcID::wire( id );
    cout << "(" << layer << "," << wire << ") lr:" << h->getFlagLR()
         << " stat:" << h->getStat() << " tk:" << h->getTrkId() << " doca:" << setw( 10 )
         << h->getDoca() << std::endl;
  } // end of hit list
}