Hough.cxx

Go to the documentation of this file.

#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IDataManagerSvc.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/IPartPropSvc.h"
 
#include "HoughTransAlg/Hough.h"
#include "EventModel/EventHeader.h"
#include "EvTimeEvent/RecEsTime.h"
#include "McTruth/McParticle.h"
#include "HepPDT/ParticleDataTable.hh"
#include "MdcTrkRecon/MdcTrackParams.h"
#include "MdcRecoUtil/Pdt.h"
#include "Identifier/Identifier.h"
#include "Identifier/CgemID.h"
#include "Identifier/MdcID.h"
#include "MdcGeom/MdcDetector.h"
#include "MdcData/MdcRecoHitOnTrack.h"
#include "CgemData/CgemHitOnTrack.h"
#include "TrkBase/TrkRep.h"
#include "TrkBase/TrkDifTraj.h"
#include "TrkBase/TrkExchangePar.h"
 
#include <math.h>
#include <iostream>
#include <fstream>
 
class TrkRep;
using namespace std;
 
HoughFinder::HoughFinder(const std::string& name, ISvcLocator* pSvcLocator):Algorithm(name, pSvcLocator)
{
    declareProperty("debug",        m_debug = 0);
    declareProperty("cgem",         m_cgem = 1);
    declareProperty("mcTruth",      m_mcTruth = 0);
    declareProperty("fillNTuple",   m_fillNTuple= 0);
    declareProperty("trackCharge",  m_trackCharge = 0);
    declareProperty("findPeakMethod",m_findPeakMethod=0);
 
    //declareProperty("nStep_xy",     m_nStep_xy = 1);
    //declareProperty("xBin_xy",      m_xBin_xy);
    //declareProperty("yBin_xy",      m_yBin_xy);
    //declareProperty("xExtend_xy",   m_xExtend_xy);
    //declareProperty("yExtend_xy",   m_yExtend_xy);
    //declareProperty("nVote_xy",     m_nVote_xy);
    //declareProperty("nRMS_xy",      m_nRMS_xy);
    //declareProperty("xInclude_xy",  m_xInclude_xy = 2);
    //declareProperty("yInclude_xy",  m_yInclude_xy = 2);
 
    //declareProperty("nStep_sz",     m_nStep_sz = 1);
    //declareProperty("xBin_sz",      m_xBin_sz);
    //declareProperty("yBin_sz",      m_yBin_sz);
    //declareProperty("xExtend_sz",   m_xExtend_sz);
    //declareProperty("yExtend_sz",   m_yExtend_sz);
    //declareProperty("nVote_sz",     m_nVote_sz);
    //declareProperty("nRMS_sz",      m_nRMS_sz);
    //declareProperty("xInclude_sz",  m_xInclude_sz = 2);
    //declareProperty("yInclude_sz",  m_yInclude_sz = 2);
 
    //read raw data setup
    declareProperty("keepBadTdc",   m_keepBadTdc = 0);
    declareProperty("dropHot",      m_dropHot = 0);
    declareProperty("keepUnmatch",  m_keepUnmatch = 0);
    declareProperty("driftTimeUpLimit",  m_driftTimeUpLimit = 1500);
    //declareProperty("",m_ = 0);
    declareProperty("pdtFile",      m_pdtFile = "pdt.table");
    
    // tests by llwang
    declareProperty("nBinTheta",    m_nBinTheta=100 );
    declareProperty("nBinRho",      m_nBinRho=50 );
    declareProperty("rhoRange",     m_rhoRange = 0.1);
    declareProperty("ExtPeak_theta",m_ExtPeak_theta= 4);
    declareProperty("ExtPeak_rho",  m_ExtPeak_rho= 4);
    declareProperty("shareHitRate", m_shareHitRate = 0.7);
    
    declareProperty("nBinTanl",     m_nBinTanl = 100);
    declareProperty("nBinDz",       m_nBinDz = 100);
    declareProperty("ExtPeak_tanl", m_ExtPeak_tanl = 1);
    declareProperty("ExtPeak_dz",   m_ExtPeak_dz = 1);
    declareProperty("filter",       m_filter= 0);
    declareProperty("eventFile",    m_evtFile= "EventList.txt");
    declareProperty("fitFlag",      m_fitFlag= 3);
    declareProperty("storeFlag",    storeFlag=1);
    declareProperty("checkHits",    m_checkHits=1);
    declareProperty("Layer",        Layer=20);
    declareProperty("clearTrack",   m_clearTrack=1);
    declareProperty("useCgemInGlobalFit",m_useCgemInGlobalFit=3);
    declareProperty("printFlag",    printFlag=0);
    declareProperty("removeNOuterHits", m_removeNOuterHits = 0);
 
    // --- cuts for circle
    declareProperty("circle_chiCut", m_circle_chiCut=25.0);
 
    // --- cuts to combine tracks
    declareProperty("drCut_combine",   m_drCut=2.0);
    declareProperty("phi0Cut_combine", m_phi0Cut=0.2);
    declareProperty("kappaCut_combine",m_kappaCut=1.5);
    declareProperty("dzCut_combine",   m_dzCut=5.0);
    declareProperty("tanlCut_combine", m_tanlCut=0.5);
 
    // --- cuts for hits saved in RecMdcTrack
    declareProperty("chi2CutHits", m_chi2CutHits=25);
 
    m_maxFireLayer = -1;
    m_totEvtProcessed=0;
}
 
StatusCode HoughFinder::initialize()
{
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "HoughFinder::initialize()" << endreq;
    StatusCode sc;
 
    // --- RawData
    IRawDataProviderSvc* irawDataProviderSvc;
    sc = service ("RawDataProviderSvc", irawDataProviderSvc);
    m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
    if ( sc.isFailure() ){
        log << MSG::ERROR << name()<<" Could not load RawDataProviderSvc!" << endreq;
        return StatusCode::FAILURE;
    }
 
 
    // --- MDC Geometry
    IMdcGeomSvc*     imdcGeomSvc;
    sc = service ("MdcGeomSvc", imdcGeomSvc);
    m_mdcGeomSvc = dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
    if ( sc.isFailure() ){
        log << MSG::ERROR << "Could not load MdcGeoSvc!" << endreq;
        return StatusCode::FAILURE;
    }
 
    // --- MdcCalibFunSvc
    IMdcCalibFunSvc* imdcCalibSvc; 
    sc = service ("MdcCalibFunSvc", imdcCalibSvc);
    m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
    if ( sc.isFailure() ){
        log << MSG::ERROR << "Could not load MdcCalibFunSvc!" << endreq;
        return StatusCode::FAILURE;
    }
 
    if(m_cgem){
        // --- CGEM geometry
        ICgemGeomSvc*    icgemGeomSvc;
        sc = service ("CgemGeomSvc", icgemGeomSvc);
        m_cgemGeomSvc = dynamic_cast<CgemGeomSvc*> (icgemGeomSvc);
        if ( sc.isFailure() ){
            log << MSG::ERROR << "Could not load CgemGeomSvc!" << endreq;
            return StatusCode::FAILURE;
        }
 
        // --- CgemCalibFunSvc
        ICgemCalibFunSvc* icgemCalibSvc;
        sc = service ("CgemCalibFunSvc", icgemCalibSvc);
        m_cgemCalibFunSvc = dynamic_cast<CgemCalibFunSvc*>(icgemCalibSvc);
        if ( sc.isFailure() ){
            log << MSG::ERROR << "Could not load CgemCalibFunSvc!" << endreq;
            return StatusCode::FAILURE;
        }
 
    }
 
    // --- MagneticFieldSvc 
    sc = service ("MagneticFieldSvc",m_pIMF);
    if( sc.isFailure() ) {
        log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
        return StatusCode::FAILURE;
    }
    m_bfield = new BField(m_pIMF);
    //log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
    m_trkContextEv = new TrkContextEv(m_bfield);
    if(m_debug==5) cout<<"field z = "<<m_bfield->bFieldNominal()<< endl;
 
    m_mdcDetector = new MdcDetector();
    //read pdt
    Pdt::readMCppTable(m_pdtFile);
 
    // --- Mdc Print Service
    //IMdcPrintSvc* imdcPrintSvc;
    //sc = service ("MdcPrintSvc", imdcPrintSvc);
    //m_mdcPrintSvc = dynamic_cast<MdcPrintSvc*> (imdcPrintSvc);
    //if ( sc.isFailure() ){
    //log << MSG::FATAL << "Could not load MdcPrintSvc!" << endreq;
    //return StatusCode::FAILURE;
    //}
 
    // --- Bes Timer Service 
    //sc = service( "BesTimerSvc", m_timersvc);
    //if( sc.isFailure() ) {
    //log << MSG::WARNING << " Unable to locate BesTimerSvc" << endreq;
    //return StatusCode::FAILURE;
    //}
    //m_timer_all = m_timersvc->addItem("Execution");
    //m_timer_all->propName("nExecution");
 
    HoughHit::setMdcGeomSvc(m_mdcGeomSvc);
    HoughHit::setMdcCalibFunSvc(m_mdcCalibFunSvc);
    HoughHit::setCgemGeomSvc(m_cgemGeomSvc);
    HoughHit::setCgemCalibFunSvc(m_cgemCalibFunSvc);
    HoughHit::setMdcDetector(m_mdcDetector);
    HoughTrack::m_clearTrack = m_clearTrack;
    HoughTrack::m_useCgemInGlobalFit = m_useCgemInGlobalFit;
 
    m_roughRhoThetaMap=TH2D("roughRhoThetaMap","roughRhoThetaMap",m_nBinTheta,0.,M_PI, m_nBinRho, -1.*m_rhoRange, m_rhoRange);
    //m_tanlRange = 0.93/sqrt(1-0.93*0.93);//FIXME
    m_tanlRange = 3.0;//FIXME
    m_dzRange = 50;
    m_roughTanlDzMap = TH2D("roughTanlDzMap","roughTanlDzMap",m_nBinTanl,-1.*m_tanlRange,m_tanlRange,m_nBinDz,-1.*m_dzRange,m_dzRange);
 
    const int nPt=12;
    double rho[nPt]         = {0.07, 0.056,0.045,0.038,0.0335,0.03,0.0198,0.0148,0.0074,0.00376,0.00303,0.00157};
    double cut1_Cgem_99[12]={-2.14,-1.36, -0.6 , -0.46, -0.35, -0.59, -0.32, -0.26, -0.22, -0.21, -0.21, -0.211};
    double cut2_Cgem_99[12]={ 0.5 , 1.77,  1.99,  1.94,  1.99,  1.9 ,  0.31,  0.27,  0.24,  0.23,  0.23,  0.222};
    double cut1_ODC1_99[12]={-1.28,-0.71, -0.58, -0.62, -0.64, -0.68, -0.18, -0.14, -0.11, -0.1 , -0.1 ,  -0.11 };
    double cut2_ODC1_99[12]={ 0.9 , 0.74,  0.42,  0.37,  0.32,  0.37,  0.28,  0.25,  0.24,  0.24,  0.24,  0.23};
    double cut1_ODC2_99[12]={ -2.14, -1.25, -1.28, -0.86, -0.47, -0.78, -0.36, -0.44, -0.61, -0.67, -0.64, -0.82 };
    double cut2_ODC2_99[12]={ -1.35, 0.55 , 0.53 , 0.83 , 0.85 , 0.83 , 0.38 , 0.55 , 0.49 , 0.46 , 0.49 , 0.4};
    m_cut1_cgem = new TGraph(nPt, rho, cut1_Cgem_99);
    m_cut2_cgem = new TGraph(nPt, rho, cut2_Cgem_99);
    m_cut1_ODC1 = new TGraph(nPt, rho, cut1_ODC1_99);
    m_cut2_ODC1 = new TGraph(nPt, rho, cut2_ODC1_99);
    m_cut1_ODC2 = new TGraph(nPt, rho, cut1_ODC2_99);
    m_cut2_ODC2 = new TGraph(nPt, rho, cut2_ODC2_99);
 
 
    int s(0),l(0),npt(0);
    double momenta[100] = {0};
    double cuts[2][43][100] = {0};
    string cutFilePath = getenv("HOUGHTRANSALGROOT");
    cutFilePath += "/share/cut.txt";
    ifstream fin(cutFilePath.c_str(), ios::in);
    //cout<<"cutfiledir:"<<cutFilePath.c_str()<<endl;
    //ifstream fin("/workfs/bes/huangzhen/workarea/CgemBoss/6.6.5.c/Reconstruction/HoughFinder/HoughFinder-00-00-13/share/cut.txt",ios::in);
    if (!fin.good())
        cout << "Error in " << __FILE__<<" when opening "<<cutFilePath.c_str()<<endl;
    bool firstline(true);
    string strcom;
    while(std::getline(fin, strcom)){
        if(firstline){
            fin >> npt
                >> momenta[0]  >> momenta[1]  >> momenta[2]  >> momenta[3]  >> momenta[4]  >> momenta[5]  >> momenta[6]  >> momenta[7]  >> momenta[8]  >> momenta[9]
                >> momenta[10] >> momenta[11] >> momenta[12] >> momenta[13] >> momenta[14] >> momenta[15] >> momenta[16] >> momenta[17] >> momenta[18] >> momenta[19]
                >> momenta[20] >> momenta[21] >> momenta[22] >> momenta[23] >> momenta[24] >> momenta[25] >> momenta[26] >> momenta[27] >> momenta[28] >> momenta[29]
                >> momenta[30] >> momenta[31] >> momenta[32] >> momenta[33] >> momenta[34] >> momenta[35] >> momenta[36] >> momenta[37] >> momenta[38] >> momenta[39]
                >> momenta[40] >> momenta[41] >> momenta[42] >> momenta[43] >> momenta[44] >> momenta[45] >> momenta[46] >> momenta[47] >> momenta[48] >> momenta[49]
                >> momenta[50] >> momenta[51] >> momenta[52] >> momenta[53] >> momenta[54] >> momenta[55] >> momenta[56] >> momenta[57] >> momenta[58] >> momenta[59]
                >> momenta[60] >> momenta[61] >> momenta[62] >> momenta[63] >> momenta[64] >> momenta[65] >> momenta[66] >> momenta[67] >> momenta[68] >> momenta[69]
                >> momenta[70] >> momenta[71] >> momenta[72] >> momenta[73] >> momenta[74] >> momenta[75] >> momenta[76] >> momenta[77] >> momenta[78] >> momenta[79]
                >> momenta[80] >> momenta[81] >> momenta[82] >> momenta[83] >> momenta[84] >> momenta[85] >> momenta[86] >> momenta[87] >> momenta[88] >> momenta[89]
                >> momenta[90] >> momenta[91] >> momenta[92] >> momenta[93] >> momenta[94] >> momenta[95] >> momenta[96] >> momenta[97] >> momenta[98] >> momenta[99]
                ;
            firstline = false;
            continue;
        }
        fin >> s >> l;
        fin >> cuts[s][l][0]  >> cuts[s][l][1]  >> cuts[s][l][2]  >> cuts[s][l][3]  >> cuts[s][l][4]  >> cuts[s][l][5]  >> cuts[s][l][6]  >> cuts[s][l][7]  >> cuts[s][l][8]  >> cuts[s][l][9]
            >> cuts[s][l][10] >> cuts[s][l][11] >> cuts[s][l][12] >> cuts[s][l][13] >> cuts[s][l][14] >> cuts[s][l][15] >> cuts[s][l][16] >> cuts[s][l][17] >> cuts[s][l][18] >> cuts[s][l][19]
            >> cuts[s][l][20] >> cuts[s][l][21] >> cuts[s][l][22] >> cuts[s][l][23] >> cuts[s][l][24] >> cuts[s][l][25] >> cuts[s][l][26] >> cuts[s][l][27] >> cuts[s][l][28] >> cuts[s][l][29]
            >> cuts[s][l][30] >> cuts[s][l][31] >> cuts[s][l][32] >> cuts[s][l][33] >> cuts[s][l][34] >> cuts[s][l][35] >> cuts[s][l][36] >> cuts[s][l][37] >> cuts[s][l][38] >> cuts[s][l][39]
            >> cuts[s][l][40] >> cuts[s][l][41] >> cuts[s][l][42] >> cuts[s][l][43] >> cuts[s][l][44] >> cuts[s][l][45] >> cuts[s][l][46] >> cuts[s][l][47] >> cuts[s][l][48] >> cuts[s][l][49]
            >> cuts[s][l][50] >> cuts[s][l][51] >> cuts[s][l][52] >> cuts[s][l][53] >> cuts[s][l][54] >> cuts[s][l][55] >> cuts[s][l][56] >> cuts[s][l][57] >> cuts[s][l][58] >> cuts[s][l][59]
            >> cuts[s][l][60] >> cuts[s][l][61] >> cuts[s][l][62] >> cuts[s][l][63] >> cuts[s][l][64] >> cuts[s][l][65] >> cuts[s][l][66] >> cuts[s][l][67] >> cuts[s][l][68] >> cuts[s][l][69]
            >> cuts[s][l][70] >> cuts[s][l][71] >> cuts[s][l][72] >> cuts[s][l][73] >> cuts[s][l][74] >> cuts[s][l][75] >> cuts[s][l][76] >> cuts[s][l][77] >> cuts[s][l][78] >> cuts[s][l][79]
            >> cuts[s][l][80] >> cuts[s][l][81] >> cuts[s][l][82] >> cuts[s][l][83] >> cuts[s][l][84] >> cuts[s][l][85] >> cuts[s][l][86] >> cuts[s][l][87] >> cuts[s][l][88] >> cuts[s][l][89]
            >> cuts[s][l][90] >> cuts[s][l][91] >> cuts[s][l][92] >> cuts[s][l][93] >> cuts[s][l][94] >> cuts[s][l][95] >> cuts[s][l][96] >> cuts[s][l][97] >> cuts[s][l][98] >> cuts[s][l][99]
            ;
        if(fin.peek()<0)break;
    }
    for(int is=0;is<=1;is++){
        for(int il=0;il<43;il++){
            TGraph* gr = new TGraph();
            int point(0);
            //cout<<il<<"  ";
            for(int ipt=0;ipt<npt;ipt++){
                //cout<<cuts[is][il][ipt]<<"  ";
                if(fabs(cuts[is][il][ipt])<1e-6)continue;
                //cut[is][il]->SetPoint(point,momenta[ipt]/1000.,cuts[is][il][ipt]);
                gr->SetPoint(point,momenta[ipt]/1000.,cuts[is][il][ipt]);
                point++;
            }
            if(point>0)HoughTrack::m_cut[is][il] = gr;
            //cout<<endl;
        }
    }
 
 
    if(m_fillNTuple) sc = bookTuple();
    if ( sc.isFailure() ){
        log << MSG::FATAL << "Could not book Tuple !" << endreq;
        return StatusCode::FAILURE;
    }
 
 
    myDotsHelixFitter.initialize();
 
    return StatusCode::SUCCESS;
}
 
StatusCode HoughFinder::execute()
{
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "HoughFinder::execute()" << endreq;
    cout.precision(6);
 
 
    //cout<<"line "<<__LINE__<<endl;
    // --- Get event header
    SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
    if(!eventHeader){
        log << MSG::ERROR << "Could not find Event Header" << endreq;
        return StatusCode::FAILURE;
    }
    m_run = eventHeader->runNumber();
    m_event = eventHeader->eventNumber();
    m_totEvtProcessed++; //if(m_totEvtProcessed%100==0) cout<<"HoughFinder: "<<m_totEvtProcessed<<" events processed."<<endl;
 
    myDotsHelixFitter.updateBField();
 
    //cout<<"line "<<__LINE__<<endl;
    
    // --- event start time
    SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
    if(!evTimeCol){
        log << MSG::ERROR << "Could not retrieve RecEsTimeCol" << endreq;
        return StatusCode::FAILURE;
    }
    m_bunchT0 = 0;
    RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
    if (iter_evt != evTimeCol->end()){
        m_bunchT0 = (*iter_evt)->getTest()*1.e-9;//m_bunchT0-s, getTest-ns
        //cout<<"this event T0 = "<<m_bunchT0<<endl;
    }
 
    // --- clear, register and store reconstructed track
    registerTrack(m_trackList_tds,m_hitList_tds);
    
 
    //cout<<"line "<<__LINE__<<endl;
    if(m_filter){
        ifstream lowPt_Evt;
        lowPt_Evt.open(m_evtFile.c_str());
        vector<int> evtlist;
        int evtNum;
        while( lowPt_Evt >> evtNum) {
            evtlist.push_back(evtNum);
        }
        vector<int>::iterator iter_evt = find(evtlist.begin(),evtlist.end(),m_event);
        if(iter_evt == evtlist.end()){
            setFilterPassed(false);
            return StatusCode::SUCCESS;
        }
        else{
            setFilterPassed(true);
        }
    }    
    //if(m_event<2000||m_event>3000){
    //setFilterPassed(false);
    //return StatusCode::SUCCESS;
    //}
 
    if(m_debug==5)
    {
        //cout<<endl;
        cout<<"===================================================================================================="<<endl;
        cout<<"run:"<<m_run<<"  , event: "<<m_event<<endl;
        cout<<"field z = "<<m_bfield->bFieldNominal()<< endl;
        //cout<<m_event<<" , ";
        //cout<<endl;
    }
 
    //cout<<"line "<<__LINE__<<endl;
    // --- prepare hits for reconstruction
    int nHit = makeHoughHitList();
 
    //cout<<"line "<<__LINE__<<endl;
    int nMcHit;
    if(m_run<0&&m_mcTruth){
        int nMcHit = getMcHitCol();
        //cout<<"nMcHit "<<nMcHit<<endl;
        if(printFlag)cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ McParticle @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"<<endl;
        getMcParticleCol();
        if(printFlag)cout<<"@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ McParticle @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@"<<endl;
        if(printFlag)cout<<endl;
        if(printFlag)cout<<endl;
        if(printFlag)cout<<endl;
        //halfCircle(m_mcTrackCol,m_mcHitCol);
    }
 
    //cout<<"line "<<__LINE__<<endl;
    searchCircle();
    //cout<<"ntracks:"<<m_houghTrackList.size()<<endl;
    //cout<<"line "<<__LINE__<<endl;
    checkHot(m_houghTrackList); // remove tracks with too many shared hits (if >70%)
    //cout<<"line "<<__LINE__<<endl;
    solveSharedHits(m_XHoughHitList); // solve hits sharing by smaller distance
    //updateXHitsStatistics(m_houghTrackList);
    //checkTrack(m_houghTrackList);
    //cout<<"line "<<__LINE__<<endl;
    associateVHits();
    //cout<<"line "<<__LINE__<<endl;
    if(printFlag)cout<<"################################################## HoughTrack ##################################################"<<endl;
    if(storeFlag==1)storeTracks(m_trackList_tds, m_hitList_tds, m_houghTrackList);
    else
    if(storeFlag==2)storeTrack(m_houghTrackList,m_trackList_tds,m_hitList_tds);
    else
    if(storeFlag==3)storeRecTracks(m_trackList_tds, m_hitList_tds, m_houghTrackList);
    if(printFlag)cout<<"################################################## HoughTrack ##################################################"<<endl;
    if(printFlag)cout<<endl;
    if(printFlag)cout<<endl;
    if(printFlag)cout<<endl;
    //cout<<"line "<<__LINE__<<endl;
 
    if(m_fillNTuple)dumpHit();
    if(m_fillNTuple==1)dumpHoughTrack();
    if(m_fillNTuple==1)dumpHoughEvent();
    if(m_fillNTuple==2)dumpTdsTrack();
    if(m_fillNTuple==2)dumpTdsEvent();
    //cout<<"line "<<__LINE__<<endl;
 
    if(printFlag)cout<<"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$  TdsTrack  $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<endl;
    if(printFlag)printTdsTrack();
    if(printFlag)cout<<"$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$  TdsTrack  $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$"<<endl;
    if(printFlag)cout<<endl;
    if(printFlag)cout<<endl;
    if(printFlag)cout<<endl;
 
    if(printFlag)cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&  HoughHit  &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl;
    if(printFlag)printHitList();
    if(printFlag)cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&  HoughHit  &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<endl;
    if(printFlag)cout<<endl;
    if(printFlag)cout<<endl;
    if(printFlag)cout<<endl;
 
    clearTrack(m_houghTrackList);
    //cout<<"line "<<__LINE__<<endl;
 
    clearMemory();
    //cout<<"line "<<__LINE__<<endl;
    if(m_debug==5) cout<<endl;
 
    return StatusCode::SUCCESS;
}// <-- execute()
 
StatusCode HoughFinder::finalize()
{
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "HoughFinder::finalize()" << endreq;
    return StatusCode::SUCCESS;
}
 
bool largerPt(HoughTrack trk1, HoughTrack trk2)
{
    return fabs(trk1.pt())>fabs(trk2.pt());
}
 
int HoughFinder::makeHoughHitList()
{
    MsgStream log(msgSvc(), name());
    int nHit(0);
    if(!(m_mdcHitCol.empty())) m_mdcHitCol.clear();
    if(!(m_houghHitList.empty()))  m_houghHitList.clear();
    if(!(m_XHoughHitList.empty())) m_XHoughHitList.clear();
    if(!(m_VHoughHitList.empty())) m_VHoughHitList.clear();
    
    if(m_cgem){
        // --- RecCgemCluster 
        SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(), "/Event/Recon/RecCgemClusterCol");
        if(!recCgemClusterCol){
            log << MSG::ERROR << "Could not retrieve RecCgemClusterCol" << endreq;
        }else{
            RecCgemClusterCol::iterator cgemClusterIter = recCgemClusterCol->begin();
            m_recCgemClusterColBegin=cgemClusterIter;
            for (;cgemClusterIter != recCgemClusterCol->end(); cgemClusterIter++,nHit++){
                const RecCgemCluster* recCgemCluster = (*cgemClusterIter);
                //m_recCgemCluster = (*cgemClusterIter);
                HoughHit hit(recCgemCluster, m_bunchT0, nHit);
                m_houghHitList.push_back(hit);
                if(hit.getLayer()>m_maxFireLayer)m_maxFireLayer = hit.getLayer();
            }
        }
    }
 
    uint32_t getDigiFlag = 0;
    if(m_dropHot)getDigiFlag |= MdcRawDataProvider::b_dropHot;
    if(m_keepBadTdc)  getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
    if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
    MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
 
    MdcDigiVec::iterator mdcDigiIter = mdcDigiVec.begin(); 
    for (;mdcDigiIter != mdcDigiVec.end(); mdcDigiIter++,nHit++){
        const MdcDigi* mdcDigi = (*mdcDigiIter);
        HoughHit hit(mdcDigi, m_bunchT0, nHit);
        if(fabs(hit.driftTime())>m_driftTimeUpLimit) continue;
        MdcHit *mdcHit= new MdcHit(mdcDigi,m_mdcDetector);
        hit.setMdcHit(mdcHit);
 
        m_mdcHitCol.push_back(mdcHit);
        m_houghHitList.push_back(hit);
        if(hit.getLayer()>m_maxFireLayer)m_maxFireLayer = hit.getLayer();
    }
 
    for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
        int flag=hitIt->getFlag();
        if(hitIt->getHitType()==HoughHit::CGEMHIT){
            if(flag==0)m_XHoughHitList.push_back(&(*hitIt));
            if(flag==2)m_VHoughHitList.push_back(&(*hitIt));
        }
        if(hitIt->getHitType()==HoughHit::MDCHIT){
            if(flag==0)m_XHoughHitList.push_back(&(*hitIt));
            else m_VHoughHitList.push_back(&(*hitIt));
        }
    }
    //for(vector<HoughHit*>::iterator hitIt = m_XHoughHitList.begin(); hitIt != m_XHoughHitList.end(); hitIt++){
    //(*hitIt)->print();
    //}
    //for(vector<HoughHit*>::iterator hitIt = m_VHoughHitList.begin(); hitIt != m_VHoughHitList.end(); hitIt++){
    //(*hitIt)->print();
    //}
 
    return nHit;
}
 
int HoughFinder::getMcHitCol()
{
    MsgStream log(msgSvc(), name());
    int nMcHit(0);
    if(!(m_mcHitCol.empty()))m_mcHitCol.clear();
 
    // --- get cgemMcHit
    if(m_cgem){
        SmartDataPtr<CgemMcHitCol> cgemMcHitCol(eventSvc(), "/Event/MC/CgemMcHitCol");
        if (!cgemMcHitCol){
            log << MSG::ERROR << "Could not retrieve CgemMcHitCol" << endreq;
        }else{
            CgemMcHitCol::iterator cgemMcHitIt = cgemMcHitCol->begin();
            for(;cgemMcHitIt != cgemMcHitCol->end();cgemMcHitIt++,nMcHit++){
                const CgemMcHit* cgemMcHit = (*cgemMcHitIt);
                HoughHit hit(cgemMcHit,m_bunchT0,nMcHit);
                m_mcHitCol.push_back(hit);
                HoughHit* phit = &(*m_mcHitCol.rbegin());
 
                // --- match fired strip ID between CemgMcHit and RecCgemCluster
                vector<int> xFireStripID;
                vector<int> vFireStripID;
                int sheet;
                TArrayI identifier = phit->getCgemMcHit()->GetIdentifier();
                for(int ii = 0; ii < identifier.GetSize(); ii++){
                    //const Identifier ident = Identifier(identifier.GetAt(ii));
                    const Identifier ident = Identifier(identifier.At(ii));
                    if(CgemID::is_xstrip(ident))xFireStripID.push_back(CgemID::strip(ident));
                    else vFireStripID.push_back(CgemID::strip(ident));
                    sheet = CgemID::sheet(ident);
                    //cout
                    //<< "   hitID:" << phit->getHitID()
                    //<< "   layerID=" << CgemID::layer(ident)
                    //<< "   sheetID=" << CgemID::sheet(ident)
                    //<< "   isXstrip=" << CgemID::is_xstrip(ident)
                    //<< "   stripID=" << CgemID::strip(ident)
                    //<<endl;
                }
                sort(xFireStripID.begin(),xFireStripID.end());
                sort(vFireStripID.begin(),vFireStripID.end());
                //vector<int>::iterator uniqueX = unique(xFireStripID.begin(),xFireStrip     ID.end());
                //vector<int>::iterator uniqueV = unique(vFireStripID.begin(),vFireStrip     ID.end());
                //xFireStripID.erase(uniqueX,xFireStripID.end());
                //vFireStripID.erase(uniqueV,vFireStripID.end());
                //for(vector<int>::iterator i = vFireStripID.begin();i!=vFireStripID.end();i++)cout<<*i<<endl;
                //cout<<endl;
 
                bool findX(false), findV(false);
                const RecCgemCluster* xCluster = NULL;
                const RecCgemCluster* vCluster = NULL;
                for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
                    if(hitIt->getHitType()!=HoughHit::CGEMHIT)continue;
                    const RecCgemCluster* recCgemCluster = hitIt->getCgemCluster();
                    if(recCgemCluster->getlayerid() != phit->getLayer())continue;
                    if(!(recCgemCluster->getflag()==0||recCgemCluster->getflag()==1))continue; 
                    //if(recCgemCluster->getsheetid() != sheet)continue;
                    //if(recCgemCluster->getTrkId() != *(phit->getTrkID().begin())continue;
                    if(!findX){
                        findX = recCgemCluster->getflag() == 0 
                            && xFireStripID.size() > 0 
                            && recCgemCluster->getclusterflagb() == *(xFireStripID.begin()) 
                            && recCgemCluster->getclusterflage() == *(xFireStripID.rbegin());
                        if(findX){
                            xCluster = recCgemCluster;
                            hitIt->setPairHit(phit);
                        }
                    }
                    if(!findV){
                        findV = recCgemCluster->getflag() == 1 
                            && vFireStripID.size() > 0 
                            && recCgemCluster->getclusterflagb() == *(vFireStripID.begin()) 
                            && recCgemCluster->getclusterflage() == *(vFireStripID.rbegin());
                        if(findV){
                            vCluster = recCgemCluster;
                            hitIt->setPairHit(phit);
                        }
                    }
 
                }
 
                /// --- put match RecCgemCluster into CemgMcHit
                if(findX&&findV){
                    for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
                        if(hitIt->getHitType()!=HoughHit::CGEMHIT)continue;
                        const RecCgemCluster* recCgemCluster = hitIt->getCgemCluster();
                        if(recCgemCluster->getlayerid() != phit->getLayer())continue;
                        if(recCgemCluster->getflag()==0||recCgemCluster->getflag()==1)continue;
                        //if(recCgemCluster->getsheetid() != sheet)continue;
                        //if(recCgemCluster->getTrkId() != *(phit->getTrkID().begin())continue;
                        if(recCgemCluster->getclusterflagb() == xCluster->getclusterid()
                        && recCgemCluster->getclusterflage() == vCluster->getclusterid()){
                            phit->setCgemCluster(recCgemCluster);
                            phit->setWire(recCgemCluster->getclusterid());
                            phit->setFlag(recCgemCluster->getflag());
                            phit->setPairHit(&(*hitIt));
                            hitIt->setPairHit(phit);
                        }
                    }
                }
                else if(findX&&!findV){
                    phit->setCgemCluster(xCluster);
                    phit->setWire(xCluster->getclusterid());
                    phit->setFlag(xCluster->getflag());
                }
                else if(!findX&&findV){
                    phit->setCgemCluster(vCluster);
                    phit->setWire(vCluster->getclusterid());
                    phit->setFlag(vCluster->getflag());
                }
            }
        }
    }
 
    // --- get mdcMcHit
    SmartDataPtr<MdcMcHitCol> mdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
    if(!mdcMcHitCol){
        log << MSG::ERROR << "Could not retrieve MdcMcHitCol" << endreq;
    }else{
        MdcMcHitCol::iterator mdcMcHitIt = mdcMcHitCol->begin();
        for(;mdcMcHitIt != mdcMcHitCol->end();mdcMcHitIt++,nMcHit++){
            const MdcMcHit* mdcMcHit = (*mdcMcHitIt);
            HoughHit hit(mdcMcHit,m_bunchT0,nMcHit);
            m_mcHitCol.push_back(hit);
            HoughHit* phit = &(*m_mcHitCol.rbegin());
            for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
                if(hitIt->getHitType()!=HoughHit::MDCHIT)continue;
                const MdcDigi* mdcDigi = hitIt->getDigi();
                if(mdcDigi->identify() == phit->getMdcMcHit()->identify()){
                    phit->setDigi(mdcDigi);
                    phit->setPairHit(&(*hitIt));
                    hitIt->setPairHit(phit);
                }
            }
        }
    }
 
    for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
        for(HitVector_Iterator mcHitIt = m_mcHitCol.begin(); mcHitIt != m_mcHitCol.end();mcHitIt++){
            if(hitIt->getLayer()!=mcHitIt->getLayer())continue;
            if(hitIt->getHitType() == HoughHit::CGEMHIT && mcHitIt->getHitType() == HoughHit::CGEMMCHIT){
                if(mcHitIt->getCgemCluster()==NULL)continue;
                int recClusterID(-1);
                int mcClusterID(-2);
                int recFlag = hitIt->getCgemCluster()->getflag();
                int mcFlag  = mcHitIt->getCgemCluster()->getflag();;
                if(recFlag==mcFlag){
                    recClusterID = hitIt->getCgemCluster()->getclusterid();
                    mcClusterID  = mcHitIt->getCgemCluster()->getclusterid();
                }
 
                if(recFlag==1&&mcFlag==2){
                    recClusterID = hitIt->getCgemCluster()->getclusterid();
                    mcClusterID  = mcHitIt->getCgemCluster()->getclusterflage();
                }
 
                if(recFlag==0&&mcFlag==2){
                    recClusterID = hitIt->getCgemCluster()->getclusterid();
                    mcClusterID  = mcHitIt->getCgemCluster()->getclusterflagb();
                }
 
                if(recFlag==2&&mcFlag==1){
                    recClusterID = hitIt->getCgemCluster()->getclusterflage();
                    mcClusterID  = mcHitIt->getCgemCluster()->getclusterid();
                }
 
                if(recFlag==2&&mcFlag==0){
                    recClusterID = hitIt->getCgemCluster()->getclusterflagb();
                    mcClusterID  = mcHitIt->getCgemCluster()->getclusterid();
                }
 
                if(recClusterID==mcClusterID){
                    hitIt->setPairHit(&(*mcHitIt));
                    hitIt->setCgemMcHit(mcHitIt->getCgemMcHit());
                    //mcHitIt->addPairHit(&(*hitIt));
 
                    //cout<<hitIt->getLayer()<<"  "<<mcHitIt->getLayer()<<"    "<<hitIt->getWire()<<"    "<<mcHitIt->getWire()<<endl;
                    //cout<<hitIt->getDigi()->identify()<<"  "<<mcHitIt->getMdcMcHit()->identify()<<endl;
                    //cout<<hitIt->getLayer()<<"    "<<mcHitIt->getLayer()<<"    ";
                    //cout<<hitIt->getWire()<<"    "<<mcHitIt->getWire()<<endl;
                    //cout<<recFlag<<"  "<<mcFlag<<"    ";
                    //cout<<recClusterID<<"  "<<mcClusterID<<endl;
 
                    //cout<<hitIt->getHitPosition();
                    //cout<<endl;
                    //cout<<mcHitIt->getHitPosition();
                    //cout<<endl;
                    //cout<<endl;
                }
            }
            else if(hitIt->getHitType() == HoughHit::MDCHIT && mcHitIt->getHitType() == HoughHit::MDCMCHIT){
                if(mcHitIt->getDigi()==NULL)continue;
                //if(hitIt->getLayer() == mcHitIt->getLayer() && hitIt->getWire() == mcHitIt->getWire())
                if(hitIt->getDigi()->identify() == mcHitIt->getDigi()->identify())
                {
                    hitIt->setPairHit(&(*mcHitIt));
                    hitIt->setMdcMcHit(mcHitIt->getMdcMcHit());
                    //mcHitIt->addPairHit(&(*hitIt));
                    //cout<<hitIt->getLayer()<<"    "<<hitIt->getWire()<<endl;
                    //cout<<mcHitIt->getLayer()<<"    "<<mcHitIt->getWire()<<endl;
                    //cout<<endl;
                    //cout<<hitIt->getDigi()->identify()<<"  "<<mcHitIt->getMdcMcHit()->identify()<<endl;
                }
            }
        }
    }
    //for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
        //if(hitIt->getHitType()!=HoughHit::MDCHIT)continue;
        //cout<<hitIt->getLayer()<<"  "<<hitIt->getWire()<<"  ";
        //if(hitIt->getPairHit()!=NULL){
            //cout<<hitIt->getPairHit()->getHitID()<<"  ";
            //cout<<hitIt->getPairHit()->getLayer()<<"  "<<hitIt->getPairHit()->getWire()<<"  ";
        //}
        //cout<<endl;
    //}
 
    return nMcHit;
}
 
int HoughFinder::getMcParticleCol(){
    const int maxCell[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 };
    m_mcTrackCol.clear();
    MsgStream log(msgSvc(), name());
    SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
    if (!mcParticleCol) {
        log << MSG::ERROR << "Could not find McParticle" << endreq;
    }else{
        bool psipDecay(false);
        McParticleCol::iterator iter_mc = mcParticleCol->begin();
        for (;iter_mc != mcParticleCol->end(); iter_mc++){
            int trackIndex  = (*iter_mc)->trackIndex();
            int pid         = (*iter_mc)->particleProperty();
            bool primaryParticle    = (*iter_mc)->primaryParticle();
            bool leafParticle       = (*iter_mc)->leafParticle();
            bool decayFromGenerator = (*iter_mc)->decayFromGenerator();
            bool decayInFlight      = (*iter_mc)->decayInFlight();
            HepLorentzVector initialPosition        = (*iter_mc)->initialPosition();
            HepLorentzVector initialFourMomentum    = (*iter_mc)->initialFourMomentum();
            //HepLorentzVector finalPosition          = (*iter_mc)->finalPosition();
            //HepLorentzVector finalFourMomentum      = (*iter_mc)->finalFourMomentum();
            unsigned int statusFlags = (*iter_mc)->statusFlags();
            //string process = (*iter_mc)->getProcess();
            //if(primaryParticle)continue;
            //if(!decayFromGenerator)continue;
            //if(pid==100443)psipDecay = true;
            //if(!psipDecay) continue;
            //if(!(fabs(pid)==11||fabs(pid)==211))continue;
 
            string name;
            int charge(0);
            if( pid == 0 ) { 
                //cout << "Wrong particle id" <<endl;
                continue;
            }else{
                IPartPropSvc* p_PartPropSvc;
                static const bool CREATEIFNOTTHERE(true);
                StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
                if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
                    log << MSG::ERROR << "Could not initialize Particle Properties Service" << endreq;
                }
                HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT(); 
                if( p_particleTable->particle(pid) ){
                    name = p_particleTable->particle(pid)->name();
                    charge = p_particleTable->particle(pid)->charge();
                }else if( p_particleTable->particle(-pid) ){
                    name = "anti " + p_particleTable->particle(-pid)->name();
                    charge = (-1)*p_particleTable->particle(-pid)->charge();
                }
            }
 
            //cout
            //<<setw(8)<<trackIndex
            //<<setw(12)<<name
            //<<setw(8)<<pid
            //<<setw(8)<<primaryParticle
            //<<setw(8)<<leafParticle
            //<<setw(8)<<decayFromGenerator
            //<<setw(8)<<decayInFlight
            //<<setw(8)<<statusFlags
            //<<endl;
 
            HoughTrack track(charge,initialPosition.v(),initialFourMomentum.v(),trackIndex);
            //HoughTrack* track = new HoughTrack(charge,initialPosition.v(),initialFourMomentum.v(),trackIndex);
            HepPoint3D pivot(0,0,0);
            track.pivot(pivot);
 
            int nHot(0);
            vector<HoughHit*> hotList;
            vector<HoughHit>& hitList = m_mcHitCol;
            for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
                vector<int> trkID = hitIt->getTrkID();
                if(trkID[0]==track.getTrkID()){
                    if(hitIt->getHitType()==HoughHit::CGEMMCHIT){
                        track.addHitPnt(&(*hitIt));
                        track.addVecStereoHitPnt(&(*hitIt));
                        hotList.push_back(&(*hitIt));
                    }
                    if(hitIt->getHitType()==HoughHit::MDCMCHIT){
                        if(hitIt->getFlag()==0)track.addHitPnt(&(*hitIt));
                        else track.addVecStereoHitPnt(&(*hitIt));
                        hotList.push_back(&(*hitIt));
                    }
                    nHot++;
                }
            }
            if(nHot>0)m_mcTrackCol.push_back(track);
 
            ///*
            bool classifyHotByHotLayer(false);
            bool classifyHotByTrackCenter(true);
            if(classifyHotByHotLayer){
                int lastHitLayer(0);
                int deltaLayer(0);
                int halfCircle(0);
                HoughHit* lastLayerHit;
                vector<HoughHit*> hitOnLayer;
                int nMdcHot(0);
                for(vector<HoughHit*>::iterator hitIt = hotList.begin(); hitIt != hotList.end(); hitIt++){
                    if((*hitIt)->getHitType()!=HoughHit::MDCHIT)continue;
                    if(nMdcHot==0){//the first hit on track
                        hitOnLayer.clear();
                        hitOnLayer.push_back(*hitIt);
                        lastLayerHit = *hitIt;
                        lastHitLayer = (*hitIt)->getLayer();
                        deltaLayer = 0;
                        halfCircle = 1;
                    }else{
                        if((*hitIt)->getLayer()==lastHitLayer){//find hits on the same layer
                            hitOnLayer.push_back(*hitIt);
                            lastHitLayer = (*hitIt)->getLayer();
                        }else{//another layer hit
                            if(deltaLayer*((*hitIt)->getLayer()-lastHitLayer)>=0){//not the turn layer
                                for(vector<HoughHit*>::iterator iter = hitOnLayer.begin(); iter != hitOnLayer.end(); iter++){
                                    (*iter)->setHalfCircle(halfCircle);
                                    lastLayerHit = *iter;
                                }
                            }else{//the turn layer
                                halfCircle++;
                                int maxWireGap(0);
                                int turnPoint(0);
                                int deltaWire(0);
                                for(vector<HoughHit*>::iterator iter = hitOnLayer.begin(); iter != hitOnLayer.end(); iter++){
                                    deltaWire = fabs((*iter)->getWire()-lastLayerHit->getWire());
                                    if(deltaWire>(maxCell[(*iter)->getLayer()])/2) deltaWire = maxCell[(*iter)->getLayer()] - deltaWire;
                                    if(deltaWire>maxWireGap){
                                        maxWireGap = deltaWire;
                                        turnPoint = iter - hitOnLayer.begin();
                                    }
                                    lastLayerHit = *iter;
                                }
                                deltaWire = fabs((*hitIt)->getWire()-lastLayerHit->getWire());
                                if(deltaWire>(maxCell[(*hitIt)->getLayer()])/2) deltaWire = maxCell[(*hitIt)->getLayer()] - deltaWire;
                                if(deltaWire>maxWireGap){
                                    maxWireGap = deltaWire;
                                    turnPoint = hitOnLayer.size();
                                }
 
                                if(maxWireGap>2){
                                    for(vector<HoughHit*>::iterator iter = hitOnLayer.begin(); iter != hitOnLayer.end(); iter++){
                                        int shift = iter-hitOnLayer.begin();
                                        if(shift<turnPoint)(*iter)->setHalfCircle(halfCircle-1);
                                        else (*iter)->setHalfCircle(halfCircle);
                                    }
                                }else{
                                    for(vector<HoughHit*>::iterator iter = hitOnLayer.begin(); iter != hitOnLayer.end(); iter++){
                                        int shift = iter-hitOnLayer.begin();
                                        if(shift<hitOnLayer.size()/2)(*iter)->setHalfCircle(halfCircle-1);
                                        else (*iter)->setHalfCircle(halfCircle);
                                    }
                                }
                            }
                            hitOnLayer.clear();
                            hitOnLayer.push_back(*hitIt);
                            lastHitLayer = (*hitIt)->getLayer();
                            deltaLayer = (*hitIt)->getLayer()-lastHitLayer;
                        }
 
                    }
                    nMdcHot++;
                }
            }
            //*/
            ///*
            if(classifyHotByTrackCenter){
                int cgemHalf(0), mdcHalf(0);
                int cgemSign(0), mdcSign(0); 
                int half(0), sign(0);
                for(vector<HoughHit*>::iterator hitIt = hotList.begin(); hitIt != hotList.end(); hitIt++){
                    //(*hitIt)->print();
                    if((*hitIt)->getHitType() == HoughHit::CGEMMCHIT){
                        if(track.judgeHalf(*hitIt) != cgemSign){
                            cgemHalf++;
                            cgemSign = track.judgeHalf(*hitIt);
                            half = cgemHalf;
                            sign = cgemSign;
                        }
                    }
                    else if((*hitIt)->getHitType() == HoughHit::MDCMCHIT){
                        if(track.judgeHalf(*hitIt) != mdcSign){
                            mdcHalf++;
                            mdcSign = track.judgeHalf(*hitIt);
                            half = mdcHalf;
                            sign = mdcSign;
                        }
                    }
                    (*hitIt)->setHalfCircle(half);
                }
            }
            if(hotList.size()>0){
                if(printFlag)track.print();
                //track.printHot();
                //cout<<endl;
            }
            for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++){
                HoughHit* hitIt = *hotIt;
                if(printFlag)hitIt->print();
                //cout<<hitIt->getHitID();
                //cout<<"  ("<<hitIt->getLayer()<<", "<<hitIt->getWire()<<")    "<<hitIt->getFlag();
                //cout<<endl;
            }
            if(printFlag)if(hotList.size())cout<<endl;
        }
        sort(m_mcTrackCol.begin(),m_mcTrackCol.end(),largerPt);
        //cout<<"N MCtrack: "<<m_mcTrackCol.size()<<endl;
        for(vector<HoughTrack>::iterator trkIter=m_mcTrackCol.begin();trkIter!=m_mcTrackCol.end();trkIter++){
            //trkIter->print();
        }
    }
    return m_mcTrackCol.size();
}
 
int HoughFinder::fillHistogram(HoughHit* hit, TH2D* hitMap, int charge, int vote)
{
    //cout<<"fillHistogram "<<__LINE__<<endl;
    int    nBin = hitMap->GetXaxis()->GetNbins();
    double xMin = hitMap->GetXaxis()->GetXmin();
    double xMax = hitMap->GetXaxis()->GetXmax();
    double yMin = hitMap->GetYaxis()->GetXmin();
    double yMax = hitMap->GetYaxis()->GetXmax();
    double dx   = (xMax-xMin)/(nBin*vote)>1e-6?(xMax-xMin)/(nBin*vote):1e-6;
    double x    = xMin + dx/2;
 
    int    flag = hit->getFlag();
    int    nPoint(0);
    if(flag==0){
        double xHit = hit->getHitPosition().x();   
        double yHit = hit->getHitPosition().y();   
        double rHit = hit->getDriftDist();//drift distance
        double denominator = xHit*xHit+yHit*yHit-rHit*rHit;
 
        //double X1(0),Y1(0),X2(0),Y2(0),R(0);
        double X1 = 2*xHit/denominator;
        double Y1 = 2*yHit/denominator;
        double X2 = 2*xHit/denominator;
        double Y2 = 2*yHit/denominator;
        double R  = 2*rHit/denominator;
 
        while(x<xMax){
            double y1 = X1*cos(x) + Y1*sin(x) + R;
            double y2 = X2*cos(x) + Y2*sin(x) - R;
            double slope1 = -X1*sin(x) + Y1*cos(x);
            double slope2 = -X2*sin(x) + Y2*cos(x);
            double hitOnCircle1 = charge*slope1*y1;
            double hitOnCircle2 = charge*slope2*y2;
 
            bool cut(true);
            cut = yMin<y1 && y1<yMax;
            cut = cut && hitOnCircle1 <=0;
            if(cut){
                hitMap->Fill(x,y1);
                nPoint++;
            }
 
        if(hit->getHitType()!=HoughHit::CGEMHIT) {
            cut = yMin<y2 && y2<yMax;
            cut = cut && hitOnCircle2 <= 0;
            if(cut){
                //if(hit.getHitType() == HoughHit::CGEMHIT)continue;// comment by llwang
                //int bin = hitMap->FindBin(x,y2);
                //if(hitMap->GetBinContent(bin)>0)continue;
                hitMap->Fill(x,y2);
                nPoint++;
            }
        }
 
        x += dx;
    }
    }else{
        vector<HoughHit::S_Z> sz = hit->getSZ();
        vector<HoughHit::S_Z>::iterator iter = sz.begin();
        for(;iter!=sz.end();iter++){
            double S = iter->first;
            double Z = iter->second;
            double dy = -S*dx;
            double y  = -dy + Z;
            while(x<xMax){
                x += dx;
                y += dy;
                //y = S*cos(x) + Z*sin(x);
                bool cut(true);
                cut = yMin<y && y<yMax;
                if(cut){
                    hitMap->Fill(x,y);
                    nPoint++;
                }
            }
        }
    }
    //cout<<"fillHistogram "<<__LINE__<<endl;
    return nPoint;
}
 
int HoughFinder::fillHistogram(HoughHit* hit, TH2D* hitMap, HoughTrack* trkCandi, int vote)
{
    int    charge = trkCandi->getCharge();
    int    nBin = hitMap->GetXaxis()->GetNbins();
    double xMin = hitMap->GetXaxis()->GetXmin();
    double xMax = hitMap->GetXaxis()->GetXmax();
    double yMin = hitMap->GetYaxis()->GetXmin();
    double yMax = hitMap->GetYaxis()->GetXmax();
    double dx   = (xMax-xMin)/(nBin*vote)>1e-6?(xMax-xMin)/(nBin*vote):1e-6;// step size
    double x    = xMin + dx/2;// first step
 
    int    flag = hit->getFlag();
    int    nPoint(0);
    if(flag==0) // X-view hits
    {
        bool   firstHalf = trkCandi->judgeHalf(hit)==1;
        if(!firstHalf) return 0;
        double xHit = hit->getHitPosition().x();   
        double yHit = hit->getHitPosition().y();   
        double rHit = hit->getDriftDist();//drift distance
        double denominator = xHit*xHit+yHit*yHit-rHit*rHit;
 
        //double X1(0),Y1(0),X2(0),Y2(0),R(0);
        double X1 = 2*xHit/denominator;
        double Y1 = 2*yHit/denominator;
        double X2 = 2*xHit/denominator;
        double Y2 = 2*yHit/denominator;
        double R  = 2*rHit/denominator;
 
        while(x<xMax){
            double y1 = X1*cos(x) + Y1*sin(x) + R;
            double y2 = X2*cos(x) + Y2*sin(x) - R;
            double slope1 = -X1*sin(x) + Y1*cos(x);
            double slope2 = -X2*sin(x) + Y2*cos(x);
            double hitOnCircle1 = charge*slope1*y1;
            double hitOnCircle2 = charge*slope2*y2;
 
            bool cut(true);
            cut = yMin<y1 && y1<yMax;
            //cut = cut && hitOnCircle1 <=0;
            if(cut){
                hitMap->Fill(x,y1);
                nPoint++;
            }
 
            cut = yMin<y2 && y2<yMax;
            //cut = cut && hitOnCircle2 <= 0;
            if(cut){
                //if(hit.getHitType() == HoughHit::CGEMHIT)continue;// comment by llwang
                //int bin = hitMap->FindBin(x,y2);
                //if(hitMap->GetBinContent(bin)>0)continue;
                hitMap->Fill(x,y2);
                nPoint++;
            }
            x += dx;
        }
    }else{ // V-view hits
        //if(hit.getPairHit()==NULL)return nPoint;// FIXME
        //if(hit.getPairHit()->getHalfCircle()!=1)return nPoint;// FIXME
        int nTangency = trkCandi->calculateZ_S(hit);
        if(nTangency==0){
            //hit.setUse(-1);// why? FIXME
        }
        else
        {
            vector<HoughHit::S_Z> sz = hit->getSZ();
            vector<HoughHit::S_Z>::iterator iter = sz.begin();
            for(;iter!=sz.end();iter++)
            {
                double S = iter->first;
                double Z = iter->second;
                while(x<xMax)
                {
                    double y  = -S*x + Z;
                    //y = S*cos(x) + Z*sin(x);
                    bool cut(true);
                    cut = yMin<y && y<yMax;
                    if(cut)
                    {
                        hitMap->Fill(x,y);
                        nPoint++;
                    }
                    x += dx;
                }
            }
        }
        //cout<<setw(8)<<""<<"nTangency="<<nTangency<<", nPointFilled="<<nPoint<<endl;
    }
    return nPoint;
}
 
int HoughFinder::fillHistogram(vector<HoughHit*> &hitList, TH2D* hitMap, int charge, int vote, HoughTrack* trkCandi )
{
    bool nearStereoHits= false;
    if(trkCandi!=NULL) nearStereoHits = trkCandi->nearStereoHits();
 
    int nHitsFilled = 0;
    m_nVClusterOnSZmap=0;
    m_VHoughHitListOnSZmap.clear();
    std::vector<HoughHit*>::iterator iter = hitList.begin();
    for(; iter!= hitList.end(); iter++)
    {
        int used = (*iter)->getUse();
        if(used==0) // too tight? FIXME
        {
            int nPoint = 0;
            //(*iter)->print();
            if(trkCandi!=NULL){
                //if((*iter)->getPairHit()==NULL) continue;// FIXME
                //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
                if(!nearStereoHits&&(*iter)->getHitType()==HoughHit::MDCHIT&&(*iter)->getFlag()!=0) continue;// if no axial hits around, the stereo hits should not be filled
 
                nPoint=fillHistogram((*iter),hitMap,trkCandi,vote);
            }
            else{
                //if((*iter)->getPairHit()==NULL) continue;// FIXME
                //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
                nPoint=fillHistogram(*iter,hitMap,charge,vote);
            }
            if(nPoint>0) {
                nHitsFilled++;
                if(trkCandi!=NULL)
                {
                    if((*iter)->getHitType()==HoughHit::CGEMHIT) m_nVClusterOnSZmap++;
                    m_VHoughHitListOnSZmap.push_back(*iter);
                }
            }
            //cout<<nPoint<<"        nPoint filled"<<(*iter)->print();
        }
    }
    return nHitsFilled;
}
 
//////////////////////////////////////////////////---Hough Transform---//////////////////////////////////////////////////
/*
//Depth First Search Peak
int HoughFinder::findTrack(vector<HoughHit> &hitList, vector<HoughTrack*>& trackVector, double x, double y, double dx, double dy, int charge, int mapType, int step)
{
step++;
if(mapType==0){
m_nStep    = m_nStep_xy;
m_xBin     = m_xBin_xy;     
m_yBin     = m_yBin_xy;     
m_xExtend  = m_xExtend_xy; 
m_yExtend  = m_yExtend_xy; 
m_nVote    = m_nVote_xy;    
m_nRMS     = m_nRMS_xy;
m_xInclude = m_xInclude_xy;
m_yInclude = m_yInclude_xy;
}else{
m_nStep    = m_nStep_sz;
m_xBin     = m_xBin_sz;     
m_yBin     = m_yBin_sz;     
m_xExtend  = m_xExtend_sz; 
m_yExtend  = m_yExtend_sz; 
m_nVote    = m_nVote_sz;    
m_nRMS     = m_nRMS_sz;
m_xInclude = m_xInclude_sz;
m_yInclude = m_yInclude_sz;
}    
 
// --- build maps
map<int,TH2D*> hitMapList;
TH2D* houghMap = buildMap(hitList,hitMapList,x,y,dx,dy,charge,mapType,step);
//cout<<"hitMapList.size()+++"<<hitMapList.size()<<"  "<<charge<<endl;
// --- find peaks
findPeak(hitList,hitMapList,houghMap,trackVector,charge,mapType,step);
// --- clear maps
clearMap(houghMap,hitMapList);
//cout<<"hitMapList.size()---"<<hitMapList.size()<<endl;
//cout<<endl<<endl;
return trackVector.size();
}
 
 
TH2D* HoughFinder::buildMap(vector<HoughHit> &hitList, map<int,TH2D*>& hitMapList, double x, double y, double dx, double dy, int charge, int mapType, int step)
{
int vote = m_nVote[step-1];
int xBin = m_xBin[step-1];
int yBin = m_yBin[step-1];
double xExtend = 0.5+m_xExtend[step-1];
double yExtend = 0.5+m_yExtend[step-1];
if(m_findPeakMethod==1){
vote = m_nVote[step-m_nStep-1];
xBin = m_xBin[step-m_nStep-1];
yBin = m_yBin[step-m_nStep-1];
xExtend = 0.5+m_xExtend[step-m_nStep-1];
yExtend = 0.5+m_yExtend[step-m_nStep-1];
}
double xMin = x-dx*xExtend;
double xMax = x+dx*xExtend;
double yMin = y-dy*yExtend;
double yMax = y+dy*yExtend;
stringstream ssname;
ssname<<"step"<<step;
string sname = ssname.str();
const char * name = sname.c_str();
TH2D* houghMap = new TH2D(name, name, xBin, xMin, xMax, yBin, yMin, yMax);
HitVector_Iterator hitIt = hitList.begin();
for(; hitIt != hitList.end();hitIt++){
//hitIt->print();
if(hitIt->getFlag()==1 && hitIt->getHitType() == HoughHit::CGEMHIT)continue;
if(mapType==0 && hitIt->getFlag() != 0)continue;
if(mapType!=0 && hitIt->getFlag() == 0)continue;
int hitID = hitIt->getHitID();
stringstream ssname;
ssname<<"step"<<step<<",hitID"<<hitID;
string sname = ssname.str();
const char * name = sname.c_str();
TH2D* hitMap = new TH2D(name, name, xBin, xMin, xMax, yBin, yMin, yMax);
int nPoint = fillHistogram((*hitIt),hitMap,charge,vote);
if(nPoint>0){
    houghMap->Add(hitMap);
    hitMapList[hitID] = hitMap;
}else delete hitMap;
}
return houghMap;
}
 
int HoughFinder::findPeak(vector<HoughHit> &hitList, map<int,TH2D*>& hitMapList, TH2D* houghMap, vector<HoughTrack*>& trackVector, int charge, int mapType, int step)
{
    int trkID = trackVector.size();
    int binx(0),biny(0),binz(0);
    houghMap->GetMaximumBin(binx,biny,binz);
    int peakHeight = houghMap->GetBinContent(binx,biny);
    int nRMS = m_nRMS[step-1];
    if(m_findPeakMethod==1)nRMS = m_nRMS[step-m_nStep-1];
    double peakCut = mapDev(houghMap,nRMS);
    int nPeaks(0);
    int nRemove(0);
    while(peakHeight>peakCut){
        houghMap->GetMaximumBin(binx,biny,binz);
        double x  = houghMap->GetXaxis()->GetBinCenter(binx);
        double y  = houghMap->GetYaxis()->GetBinCenter(biny);
        double dx = houghMap->GetXaxis()->GetBinWidth(binx);
        double dy = houghMap->GetYaxis()->GetBinWidth(biny);
        //cout<<x<<"  "<<y<<endl;
        centroid(houghMap,m_xInclude,m_yInclude,x,y);
        //cout<<x<<"  "<<y<<endl;
        if(step<m_nStep){
            findTrack(hitList,trackVector,x,y,dx,dy,charge,mapType,step);
            nRemove = removeHitMap(hitMapList,houghMap,trackVector);
            if(nRemove==0)nRemove = removeHitMap(hitMapList,houghMap,binx,biny);
        }else{
            if(mapType == 0){
                HoughTrack* track = new HoughTrack(charge,x,y,dx,dy,trkID);
                //track->print();
                int nHot = track->findHot(hitList,mapType);
                if(nHot>=3){
                    trackVector.push_back(track);
                    nRemove = removeHitMap(hitMapList,houghMap,track);
                    if(nRemove==0)nRemove = removeHitMap(hitMapList,houghMap,binx,biny);
                    trkID++;
                }else{
                    nRemove = removeHitMap(hitMapList,houghMap,binx,biny);
                }
 
            }else if(trackVector.size()>0){
                HoughTrack* track = *(trackVector.begin());
                double tanl = track->tanl();
                double dz   = track->dz();
                double dTanl    = track->getDTanl();
                double dDz      = track->getDDz();
                vector<HoughHit> hot;
                int nHot1 = track->findHot(hitList,hot,mapType);
                track->setTanl(x);
                track->setDz(y);
                track->setDTanl(dx);
                track->setDDz(dy);
                track->updateHelix();
                int nHot2 = track->findHot(hitList,hot,mapType);
                if(nHot1>nHot2){
                    track->setTanl(tanl);
                    track->setDz(dz);
                    track->setDTanl(dTanl);
                    track->setDDz(dDz);
                    track->updateHelix();
                }
                nRemove = removeHitMap(hitMapList,houghMap,binx,biny);
            }
        }
        //double xKurtosis = houghMap->GetKurtosis(1);
        //double yKurtosis = houghMap->GetKurtosis(2);
        houghMap->GetMaximumBin(binx,biny,binz);
        peakHeight = houghMap->GetBinContent(binx,biny);
        peakCut = mapDev(houghMap,nRMS);
        nPeaks++;
        //if(mapType!=0 && nPeaks>0)break;
        if(nRemove==0)break;
    }
    return nPeaks;
}
 
//remove map of hits that on the peak cell
int HoughFinder::removeHitMap(map<int,TH2D*>& hitMapList, TH2D* houghMap,int binx, int biny)
{
    bool debug(0);
    int nRemove(0);
    //if(debug){
    //cout<<endl;
    //for(map<int,TH2D*>::iterator iter = hitMapList.begin();iter != hitMapList.end();iter++){
    //cout<<"hit->"<<iter->first<<"    ";
    //if(iter->second==NULL)cout<<"erased"<<endl;
    //else cout<<iter->second->GetBinContent(binx,biny)<<endl;
    //}
    //cout<<endl;
    //}
 
    //int binx(0),biny(0),binz(0);
    //houghMap->GetMaximumBin(binx,biny,binz);
    map<int,TH2D*>::iterator iter = hitMapList.begin();
    for(;iter != hitMapList.end();iter++){
        if(debug)cout<<"hit->"<<iter->first<<"    ";
        if(iter->second!=NULL){
            if(iter->second->GetBinContent(binx,biny)>0){
                houghMap->Add(iter->second,-1);
                delete iter->second;
                iter->second = NULL;
                //hitMapList.erase(iter);
                //hitMapList.erase(iter->first);
                nRemove++;
                if(debug)cout<<"    erasing";
            }
        }
        else{
            if(debug)cout<<"    erased";
        }
        if(debug)cout<<endl;
    }
 
    //if(debug){
    //cout<<endl;
    //for(map<int,TH2D*>::iterator iter = hitMapList.begin();iter != hitMapList.end();iter++){
    //cout<<"hit->"<<iter->first<<"    ";
    //if(iter->second==NULL)cout<<"erased"<<endl;
    //else cout<<iter->second->GetBinContent(binx,biny)<<endl;
    //}
    //cout<<endl;
    //}
 
    return nRemove;
}
 
//remove map of hits that on the track
int HoughFinder::removeHitMap(map<int,TH2D*>& hitMapList, TH2D* houghMap, HoughTrack* track)
{
    bool debug(1);
    int nRemove(0);
    //if(debug){
    //cout<<endl;
    //for(map<int,TH2D*>::iterator iter = hitMapList.begin();iter != hitMapList.end();iter++){
    //cout<<"hit->"<<iter->first<<"    ";
    //if(iter->second==NULL)cout<<"erased"<<endl;
    //else{
    //int bin = iter->second->FindBin(track->getAngle(),track->getRho());
    //cout<<iter->second->GetBinContent(bin)<<endl;
    //}
    //}
    //cout<<endl;
    //}
 
    vector<HoughHit>* hot = track->getHot();
    HitVector_Iterator hotIt = hot->begin();
    for(;hotIt!=hot->end();hotIt++){
        int hitID = hotIt->getHitID();
        if(debug)cout<<"hit->"<<hitID<<"    ";
        map<int,TH2D*>::iterator iter = hitMapList.find(hitID);
        if(iter->second!=NULL){
            if(iter != hitMapList.end()){
                houghMap->Add(iter->second,-1);
                delete iter->second;
                iter->second = NULL;
                //hitMapList.erase(iter);
                //hitMapList.erase(iter->first);
                nRemove++;
                if(debug)cout<<"erasing"<<endl;
            }
        }else{
            if(debug)cout<<"erased"<<endl;
        }
    }
 
    //if(debug){
    //cout<<endl;
    //for(map<int,TH2D*>::iterator iter = hitMapList.begin();iter != hitMapList.end();iter++){
    //cout<<"hit->"<<iter->first<<"    ";
    //if(iter->second==NULL)cout<<"erased"<<endl;
    //else{
    //int bin = iter->second->FindBin(track->getAngle(),track->getRho());
    //cout<<iter->second->GetBinContent(bin)<<endl;
    //}
    //}
    //cout<<endl;
    //}
    return nRemove;
}
 
//remove map of hits that on all tracks
int HoughFinder::removeHitMap(map<int,TH2D*>& hitMapList, TH2D* houghMap, vector<HoughTrack*>& trackVector)
{
    int nRemove(0);
    TrackVector_Iterator trkIter = trackVector.begin();
    for(;trkIter!=trackVector.end();trkIter++){
        if((*trkIter)->getFlag() == 1)continue;
        vector<HoughHit>* hot = (*trkIter)->getHot();
        HitVector_Iterator iter = hot->begin();
        for(;iter!=hot->end();iter++){
            int hitID = iter->getHitID();
            map<int,TH2D*>::iterator iter = hitMapList.find(hitID);
            if(iter != hitMapList.end()){
                houghMap->Add(iter->second,-1);
                delete iter->second;
                //iter->second = NULL;
                hitMapList.erase(iter);
                nRemove++;
            }
        }
    }
    //for(map<int,TH2D*>::iterator iter = hitMapList.begin();iter != hitMapList.end();iter++){
    //if(iter->second!=NULL)cout<<iter->first<<"  "<<iter->second->GetBinContent(binx,biny)<<endl;
    //else cout<<iter->first<<"  "<<"  erase"<<endl;
    //}
    //cout<<endl;
    return nRemove;
}
 
int HoughFinder::clearMap(TH2D* houghMap, map<int,TH2D*>& hitMapList)
{
    int nRemove(0);
    map<int,TH2D*>::iterator iter = hitMapList.begin();
    for(;iter != hitMapList.end();iter++){
        if(iter->second != NULL){
            delete iter->second;
            iter->second = NULL;
            //hitMapList.erase(iter);
            nRemove++;
        }
    }
    hitMapList.clear();
 
    delete houghMap;
    houghMap = NULL;
 
    return nRemove;
}
 
 
int HoughFinder::centroid(TH2D* houghMap, int xInclude, int yInclude, double &xc, double &yc)
{
    double X(0),Y(0),W(0);
    int binx(0),biny(0),binz(0);
    houghMap->GetMaximumBin(binx,biny,binz);
    int xMin = binx-xInclude>0?binx-xInclude:0;
    int xMax = binx+xInclude<houghMap->GetNbinsX()?binx+xInclude:houghMap->GetNbinsX();
    int yMin = biny-yInclude>0?biny-yInclude:0;
    int yMax = biny+yInclude<houghMap->GetNbinsY()?biny+yInclude:houghMap->GetNbinsY();
    //cout<<houghMap->GetNbinsX()<<"    "<<houghMap->GetNbinsY()<<endl;
    for(int i=binx-xInclude;i<=binx+xInclude;i++){
        for(int j=biny-yInclude;j<=biny+yInclude;j++){
            double x = houghMap->GetXaxis()->GetBinCenter(i);
            double y = houghMap->GetYaxis()->GetBinCenter(j);
            int    w = houghMap->GetBinContent(i,j);
            X += w*x;
            Y += w*y;
            W += w;
        }
    }
    if(W>0){
        xc = X/W;
        yc = Y/W;
    }
    return W;
}
 
double HoughFinder::mapDev(TH2D* houghMap,double nRMS)
{
    int xBins = houghMap->GetXaxis()->GetNbins();
    int yBins = houghMap->GetYaxis()->GetNbins();
    int count_use=0;
    double sum=0;
    double sum2=0;
    for(int i=0;i<xBins;i++){
        for(int j=0;j<yBins;j++){
            int count=houghMap->GetBinContent(i+1,j+1);
            if(count!=0){
                count_use++;
                sum+=count;
                sum2+=count*count;
            }
        }
    }
    double aver(0),rms(0);
    if(count_use != 0){
        aver=sum/count_use;
        if(count_use==1) rms = sqrt(sum2/count_use-aver*aver);
        else rms = sqrt((sum2-count_use*aver*aver)/(count_use-1));
    }
    double cut = aver + nRMS*rms;
    return cut;
}
//*/
/*
//Breadth First Search Peak
int HoughFinder::findTrack(vector<HoughTrack>& trackVector, int charge, int mapType, int step)
{
step++;
vector<HoughTrack>& trkVector;
TrackVector_Iterator trkIter = trackVector.begin();
for(;trkIter!=trackVector.end();trkIter++){
int    xBin = 10;
int    yBin = 10;
double xPeak = trkIter->getAlpha();
double yPeak = trkIter->getRho();
double dx(0),dy(0);
if(mapType==0){
dx = trkIter->getDAngle();
dy = trkIter->getDRho();
}else{
dx = trkIter->getDTanl();
dy = trkIter->getDDz();
}
findTrack(trackVector,xPeak-dx,xPeak+dx,yPeak-dy,yPeak+dy,charge,mapType,m_nStep+step);
}
trackVector.swap(trkVector);
if(step<m_nStep)findTrack(trackVector);
}
 
//Breadth First Search Peak
int HoughFinder::findTrack(vector<HoughTrack>& trackVector, int charge, int mapType, int step)
{
while(step<m_nStep){
step++;
vector<HoughTrack>& trkVector;
trackVector.swap(trkVector);
trackVector.clear();
TrackVector_Iterator trkIter = trkVector.begin();
for(;trkIter!=trkVector.end();trkIter++){
int    xBin = 10;
int    yBin = 10;
double xPeak = trkIter->getAlpha();
double yPeak = trkIter->getRho();
double dx(0),dy(0);
if(mapType==0){
dx = trkIter->getDAngle();
dy = trkIter->getDRho();
}else{
dx = trkIter->getDTanl();
dy = trkIter->getDDz();
}
findTrack(trackVector,xPeak-dx,xPeak+dx,yPeak-dy,yPeak+dy,charge,mapType,m_nStep+step);
}
}
//trackVector.swap(trkVector);
if(step<m_nStep)findTrack(trackVector);
}
//*/
 
bool moreHot(HoughTrack trk1, HoughTrack trk2)
{
    int nHot1 = trk1.getHotList().size();
    int nHot2 = trk2.getHotList().size();
    return nHot1>nHot2;
    //return trk1.getHot().size() > trk2.getHot().size();
}
 
//check track by hots, remove the shorter track with too many shared hits
// --- v1 ---
/*
   int HoughFinder::checkHot(vector<HoughTrack>& trackVector)
   {
   int nDelete(0);
   std::sort(trackVector.begin(),trackVector.end(),moreHot);// sort track candidates from more hits to less
   for(vector<HoughTrack>::iterator trkIT1 = trackVector.begin(); trkIT1!=trackVector.end(); trkIT1++){
   if((trkIT1)->getFlag() == 1)continue;
   vector<HoughHit>* hot1 = (trkIT1)->getHot();
   trkIT1->dropRedundentCgemXHits();
//cout<<"hit size "<<hot1->size()<<endl;
if(hot1->size()<3){
(trkIT1)->setFlag(1);
trkIT1->releaseSelHits();
continue;
}
 
for(vector<HoughTrack>::iterator trkIT2 = trkIT1+1; trkIT2!=trackVector.end();trkIT2++){
if((trkIT2)->getFlag() == 1)continue;
vector<HoughHit>* hot2 = (trkIT2)->getHot();
if(hot2->size()<3){
(trkIT2)->setFlag(1);
trkIT2->releaseSelHits();
continue;
}
 
int shareHit(0);
for(HitVector_Iterator hitIt1 = hot1->begin(); hitIt1 != hot1->end(); hitIt1++){
for(HitVector_Iterator hitIt2 = hot2->begin(); hitIt2 != hot2->end(); hitIt2++){
//cout<<"id2, id1 = "<<hitIt2->getHitID()<<", "<<hitIt1->getHitID()<<endl;
if(hitIt2->getHitID()==hitIt1->getHitID())shareHit++;
}
}
//cout<<"shareHit = "<<shareHit<<endl;
if(shareHit>m_shareHitRate*hot2->size()){
//delete *trkIT2;
//trackVector.erase(trkIT2);
//trkIT2--;
(trkIT2)->setFlag(1);
trkIT2->releaseSelHits();
nDelete++;
cout<<"too many shared hits ("<<shareHit<<"/"<<hot2->size()
<<") in track "<<(trkIT2)->getTrkID()
<<" with track "<<trkIT1->getTrkID()<<endl;
}
}
}
return nDelete;
}
//*/
 
// --- v2 ---
int HoughFinder::checkHot(vector<HoughTrack>& trackVector)
{
    //m_debug=true;
    int nDelete(0);
    if(trackVector.size()==0)return 0;
    std::sort(trackVector.begin(),trackVector.end(),moreHot);// sort track candidates from more hits to less
 
    vector<HoughTrack>::iterator trkIT1 = trackVector.end(); 
    //vector<HoughTrack>::iterator trkIT1 = trackVector.rbegin(); 
    bool loop = true;
 
    while(loop)
    {
        trkIT1--;
        if(trkIT1==trackVector.begin()) loop=false;
        if((trkIT1)->getFlag() == 1) continue;
        //trkIT1->dropRedundentCgemXHits();
        //cout<<"hit size "<<hot1->size()<<endl;
        int nHits = trkIT1->getVecHitPnt().size();
        if(nHits<3){
            (trkIT1)->setFlag(1);
            trkIT1->releaseSelHits();
            if(m_debug==1) {
                cout<<"too less hits ("<<nHits<<") in track "<<(trkIT1)->getTrkID()<<endl;
            }
            continue;
        }
        int nShared = trkIT1->getNHitsShared();
 
        //cout<<"shareHit = "<<shareHit<<endl;
        if(nShared>m_shareHitRate*nHits){
            trkIT1->setFlag(1);
            trkIT1->releaseSelHits();
            nDelete++;
            if(m_debug==1) {
                cout<<"too many shared hits ("<<nShared<<"/"<<nHits
                    <<") in track "<<(trkIT1)->getTrkID()
                    <<endl;
            }
        }
    }
    //m_debug=false;
    return nDelete;
}
 
int HoughFinder::checkTrack(vector<HoughTrack>& trackVector)
{
    int nDelete(0);
    double dDr(999);
    double dPhi0(999);
    double dKappa(999);
    double dDz(999);
    double dTanl(999);
 
    //m_debug=true;
 
    std::sort(trackVector.begin(),trackVector.end(),moreHot);// sort track candidates from more hits to less
    for(vector<HoughTrack>::iterator trkIT1 = trackVector.begin(); trkIT1!=trackVector.end(); trkIT1++){
        if((trkIT1)->getFlag() == 1)continue;
        int charge = (trkIT1)->getCharge();
        double dr      = (trkIT1)->dr();
        double phi0    = (trkIT1)->phi0();
        double kappa   = (trkIT1)->kappa();
        double dz      = (trkIT1)->dz();
        double tanl    = (trkIT1)->tanl();
        for(vector<HoughTrack>::iterator trkIT2 = trkIT1+1; trkIT2!=trackVector.end();trkIT2++){
            if((trkIT2)->getFlag() == 1)continue;
            bool sameTrack(false);
            dDr      = fabs(dr - (trkIT2)->dr());
            dPhi0    = fabs(phi0 - (trkIT2)->phi0());
            if((trkIT2)->getCharge() == charge){
                dKappa   = fabs(kappa - (trkIT2)->kappa());
                dTanl    = fabs(tanl - (trkIT2)->tanl());
            }else{
                dKappa = fabs(fabs(kappa) - fabs((trkIT2)->kappa()));
                dTanl  = fabs(fabs(tanl) - fabs((trkIT2)->tanl()));
            }
            double r = fabs(dz)<fabs((trkIT2)->dz())?(trkIT1)->radius():(trkIT2)->radius();
            double k = dz<fabs((trkIT2)->dz())?tanl:(trkIT2)->tanl();
            dDz      = fabs(dz - (trkIT2)->dz());
            int nCircle = fabs(dDz)/(2*M_PI*r*k);
            dDz = fabs(dDz - nCircle*(2*M_PI*r*k));
            if(dDz>M_PI*r*k)dDz = fabs(dDz - 2*M_PI*r*k);
            if(dPhi0<m_phi0Cut&&dKappa<m_kappaCut&&dTanl<m_tanlCut) sameTrack=true;
            //sameTrack = sameTrack&&(dDr<m_drCut);//FIXME
            //sameTrack = sameTrack&&(dPhi0<m_phi0Cut);//FIXME
            //sameTrack = sameTrack&&(dKappa<m_kappaCut);//FIXME
            //sameTrack = sameTrack&&(dDz<m_dzCut);//FIXME
            //sameTrack = sameTrack&&(dTanl<m_tanlCut);//FIXME
            if(sameTrack){
                //delete *trkIT2;
                //trackVector.erase(trkIT2);
                //trkIT2--;
                (trkIT2)->setFlag(1);
                if(m_debug==5) {
                    cout<<"dDr, dPhi0, dKappa, dDz, dTanl "
                        <<setw(10)<<dDr
                        <<setw(10)<<dPhi0
                        <<setw(10)<<dKappa
                        <<setw(10)<<dDz
                        <<setw(10)<<dTanl
                        <<endl;
                    cout<<"similar track "<<trkIT2->getTrkID()<<" is removed from hough track list"<<endl;
                }
 
                // --- associate the hits of track 2 to track 1 if within pi/2 // to-do
                nDelete++;
            }
        }
    }
    if(m_debug==5) cout<<"checkTrack(): "<<nDelete<<" similar tracks are removed from hough track list"<<endl;
    //m_debug=false;
    return nDelete;
}
 
/*
//fit tracks on xy plane, calculate S and Z of hits as well as add them to tracks, then calculate tanl and dz by Hough Transform and fit tracks again
int HoughFinder::completeTrack(vector<HoughTrack*>& trackVector)
{
int nTrk(0);
for(TrackVector_Iterator trackIT = trackVector.begin();trackIT!=trackVector.end();trackIT++){
if((*trackIT)->getFlag() == 1)continue;
 
//cout<<"begin circlr fit "<<endl;
TrkErrCode trkErrCode = (*trackIT)->fitCircle(m_mdcDetector,m_trkContextEv,m_bunchT0);
//cout<<"finish circlr fit "<<endl;
if(trkErrCode.failure()){
(*trackIT)->setFlag(1);
continue;
}
//cout<<"finish circlr fit "<<endl;
 
vector<HoughHit> stereoHot;
for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++){
HoughHit hit(*hitIt);
if((*trackIT)->calculateZ_S(hit)>0)stereoHot.push_back(hit);
}
if(stereoHot.size()<2){//FIXME
(*trackIT)->setFlag(1);
continue;
}
//cout<<"finish stereo hot finding "<<endl;
 
double tanl  = 0;
double dz    = 0;
//double dTanl = 0.93/sqrt(1-0.93*0.93);
double dTanl = 3.0;
double dDz   = 50.0;
int  mapType = 1;
int     step = 0;
int   charge = 0;
vector<HoughTrack*> trkVector;
trkVector.push_back(*trackIT);
//cout<<"before trk finding "<<endl;
findTrack(stereoHot,trkVector,tanl,dz,dTanl,dDz,charge,mapType,step);
//cout<<"end trk finding "<<endl;
TrackVector_Iterator trkIT = trkVector.begin();
HoughTrack* trk = (*trkIT);
int nHot = trk->findHot(stereoHot,mapType);
if(nHot<2){//FIXME
(*trackIT)->setFlag(1);
continue;
}
trkErrCode = trk->fitHelix(m_mdcDetector,m_trkContextEv,m_bunchT0,m_mdcHitCol);
if(trkErrCode.failure()){
(*trackIT)->setFlag(1);
continue;
}
nTrk++;
}
return nTrk;
}
//*/
StatusCode HoughFinder::registerTrack(RecMdcTrackCol*& recMdcTrackCol ,RecMdcHitCol*& recMdcHitCol)
{
    MsgStream log(msgSvc(), name());
    StatusCode sc;
    IDataProviderSvc* eventSvc = NULL;
    Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
    if (eventSvc) {
        log << MSG::INFO << "makeTds:event Svc has been found" << endreq;
    } else {
        log << MSG::FATAL << "makeTds:Could not find eventSvc" << endreq;
        return StatusCode::FAILURE;
        //return StatusCode::SUCCESS;
    }
    IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*>(eventSvc);;
 
 
    // --- register RecMdcTrack
    recMdcTrackCol = new RecMdcTrackCol;
    DataObject *aRecMdcTrackCol;
    eventSvc->findObject("/Event/Recon/RecMdcTrackCol",aRecMdcTrackCol);
    if(aRecMdcTrackCol != NULL) {
        dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
        eventSvc->unregisterObject("/Event/Recon/RecMdcTrackCol");
    }
 
    sc =  eventSvc->registerObject("/Event/Recon/RecMdcTrackCol", recMdcTrackCol);
    if(sc.isFailure()) {
        log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
        return StatusCode::FAILURE;
    }
    log << MSG::INFO << "RecMdcTrackCol registered successfully!" <<endreq;
 
 
 
    recMdcHitCol = new RecMdcHitCol;
    DataObject *aRecMdcHitCol;
    eventSvc->findObject("/Event/Recon/RecMdcHitCol",aRecMdcHitCol);
    if(aRecMdcHitCol != NULL) {
        dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
        eventSvc->unregisterObject("/Event/Recon/RecMdcHitCol");
    }
 
    sc =  eventSvc->registerObject("/Event/Recon/RecMdcHitCol", recMdcHitCol);
    if(sc.isFailure()) {
        log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
        return StatusCode::FAILURE;
    }
    log << MSG::INFO << "RecMdcHitCol registered successfully!" <<endreq;
 
 
    return StatusCode::SUCCESS;
}//end of stroeTracks
 
int HoughFinder::storeTrack(vector<HoughTrack>& trackVector, RecMdcTrackCol*& recMdcTrackCol ,RecMdcHitCol*& recMdcHitCol)
{
    MdcTrackParams mdcTrackParams;
    //mdcTrackParams.lPrint=8;
    mdcTrackParams.lRemoveInActive=1;
    mdcTrackParams.maxGapLength=99;
    mdcTrackParams.nHitDeleted=99;
    MdcTrackList mdcTrackList(mdcTrackParams);
    mdcTrackList.setNoInner(true);
    vector<MdcTrack*> mdcTrackVector;
    for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
        if((trkIT)->getFlag() == 1)continue;
        if((trkIT)->getFlag() == -1)continue;
        if((trkIT)->getFlag() == -2)continue;
        //(trkIT)->sortHot();
        MdcTrack* mdcTrack = new MdcTrack((trkIT)->getTrkRecoTrk());
        mdcTrackList.append(mdcTrack);
        mdcTrackVector.push_back(mdcTrack);
    }
    int nDeleted(0);
    mdcTrackList.setNoInner(true);
    if(m_checkHits)nDeleted = mdcTrackList.arbitrateHits();
    if(nDeleted>0)cout<<"nDeleted "<<nDeleted<<endl;
 
    int trkID=0;
    for(int l=0;l<mdcTrackList.length();l++){
        mdcTrackList[l]->storeTrack(trkID,recMdcTrackCol,recMdcHitCol,4);
        trkID++;
    }
 
    vector<MdcTrack*>::iterator iter = mdcTrackVector.begin();
    for(;iter!=mdcTrackVector.end();iter++){
        if(m_clearTrack)delete *iter;
    }
    return trkID;
}
 
bool sortCluster(const RecCgemCluster* clusterA , const RecCgemCluster* clusterB){
    return clusterA->getlayerid()<clusterB->getlayerid();
}
 
int HoughFinder::storeRecTracks(RecMdcTrackCol* trackList, RecMdcHitCol* hitList, vector<HoughTrack>& trackVector)
{
    bool debug=false; //debug=true;
    int trackId(0);
    int tkStat =4;
    if(debug) cout<<"N Rectrack: "<<trackVector.size()<<endl;
    sort(trackVector.begin(),trackVector.end(),largerPt);
 
    for(vector<HoughTrack>::iterator trkIT = trackVector.begin(); trkIT!=trackVector.end(); trkIT++)
    {
 
        if(debug) cout<<"trk flag: "<<(trkIT)->getFlag()<<endl;
 
        //if((trkIT)->getFlag() == 1)continue;
        //if((trkIT)->getFlag() == -1)continue;
        //if((trkIT)->getFlag() == -2)continue;
        if((trkIT)->getFlag()!=0) continue;
 
        if(debug) 
        {
            cout<<"the following track will be saved: "<<endl;
            trkIT->print();
        }
 
        //new a Rec. Track MdcTrack 
        RecMdcTrack* recMdcTrack = new RecMdcTrack();
 
        recMdcTrack->setTrackId(trackId);
 
        double helixPar[5];
        helixPar[0]=trkIT->getDr();
        helixPar[1]=trkIT->getPhi0();
        helixPar[2]=trkIT->getKappa();
        helixPar[3]=trkIT->getDz();
        helixPar[4]=trkIT->getTanl();//helixPar[4]+=0.1;
        recMdcTrack->setHelix(helixPar);
 
        int    q     = trkIT->getCharge();
        double pxy   = trkIT->pt();
        double px    = trkIT->momentum(0).x();
        double py    = trkIT->momentum(0).y();
        double pz    = trkIT->momentum(0).z();
        double p     = trkIT->momentum(0).mag();
        double theta = trkIT->direction(0).theta();
        double phi   = trkIT->direction(0).phi();
        HepPoint3D poca  = trkIT->x(0);
        HepPoint3D pivot = trkIT->pivot();
        double r = poca.perp();
        HepSymMatrix Ea = trkIT->Ea();
        //cout<<"Ea="<<Ea<<endl;
        double errorMat[15];
        int k = 0;
        for (int ie = 0 ; ie < 5 ; ie ++){
            for (int je = ie ; je < 5 ; je ++){
                errorMat[k] = Ea[ie][je];
                k++;
            }
        }
        // ---> a Helix for DotsHelixFitter
        HepVector a(5,0);
        a(1) =trkIT->getDr();
        a(2) =trkIT->getPhi0();
        a(3) =trkIT->getKappa();
        a(4) =trkIT->getDz();
        a(5) =trkIT->getTanl();
        KalmanFit::Helix aHelix(pivot,a,Ea); // <-- a Helix for DotsHelixFitter ---
        //myDotsHelixFitter.setInitialHelix(aHelix);// set a Helix for DotsHelixFitter
        double chisq    = trkIT->getChi2();
        recMdcTrack->setCharge(q);
        recMdcTrack->setPxy(pxy);
        recMdcTrack->setPx(px);
        recMdcTrack->setPy(py);
        recMdcTrack->setPz(pz);
        recMdcTrack->setP(p);
        recMdcTrack->setTheta(theta);
        recMdcTrack->setPhi(phi);
        recMdcTrack->setPoca(poca);
        recMdcTrack->setX(poca.x());//poca
        recMdcTrack->setY(poca.y());
        recMdcTrack->setZ(poca.z());
        recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
        HepPoint3D apivot(0.,0.,0.);
        recMdcTrack->setPivot(apivot);
        recMdcTrack->setVX0(0.);//pivot 
        recMdcTrack->setVY0(0.);
        recMdcTrack->setVZ0(0.);
        recMdcTrack->setError(errorMat);
        recMdcTrack->setError(Ea);
        recMdcTrack->setChi2(chisq);
        recMdcTrack->setStat(tkStat);
 
        //-----------hit list-------------
        int hitId = 0;
        HitRefVec  hitRefVec;
        ClusterRefVec clusterRefVec;
        //vector<int> vecHits;
        map<int,int> clusterFitStat;
        //int maxLayer = 0;
        int maxLayerId = -1;
        int minLayerId = 43;
        double fiTerm = 999.;
        HoughHit* fiTermHot = NULL;
 
        //vector<HoughHit*> hotList = trkIT->getVecHitPnt();
        vector<HoughHit*> hotList = trkIT->getHotList(2);
        for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++)
        {
            //cout<<"handle hit "<<hitId<<endl;
            HoughHit* hot(*hotIt);
            //int layer = (**hotIt).getLayer();
            //int wire  = (**hotIt).getWire();
            //cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
 
            /*
               if(hot->getHitType()==HoughHit::MDCHIT)
               {
            //if(hot->getUse()!=1) continue;
            //cout<<"create a new RecMdcHit"<<endl;
 
            // --- make a RecMdcHit from digi, aHelix using DotsHelixFitter
            RecMdcHit* recMdcHit = new RecMdcHit(myDotsHelixFitter.makeRecMdcHit(hot->getDigi(), aHelix));
            recMdcHit->setId(hitId);
            recMdcHit->setTrkId(trackId);
 
            // 
            // --> make a RecMdcHit from a HoughHit (parameters from Hough)
            RecMdcHit* recMdcHit = new RecMdcHit;
            recMdcHit->setId(hitId);
            recMdcHit->setTrkId(trackId);
            HepPoint3D hitPoint = hot->getHitPosition();
            double dPhi = trkIT->dPhi(hitPoint);
            HepPoint3D position  = trkIT->x(dPhi);
            Hep3Vector direction = trkIT->direction(dPhi);
            //cout<<"get 1"<<endl;
 
            double ddl      = -1*hot->getDriftDist();
            double ddr      = hot->getDriftDist();
            //double erddl    = hot->hitSigma();
            //double erddr    = hot->hitSigma();
            double pChisq   = hot->getResidual()*hot->getResidual();//FIXME
            int    lr       = 0;
            int    stat     = hot->getUse();
            stat = 1; // llwang
            //cout<<"hit stat = "<<stat<<endl;
            //cout<<"get 2"<<endl;
            Identifier mdcId = hot->getDigi()->identify();
            //cout<<"get 2.1"<<endl;
            double tdc      = hot->getDigi()->getTimeChannel();
            //cout<<"get 2.2"<<endl;
            double adc      = hot->getDigi()->getChargeChannel();
            //cout<<"get 2.3"<<endl;
            double driftT   = hot->driftTime();
            //cout<<"get 2.4"<<endl;
            double doca     = hot->getDriftDist() +hot->getResidual();
            //cout<<"get 2.5"<<endl;
            double entra    = direction.phi()-position.phi();//FIXME
            //cout<<"get 2.6"<<endl;
            while(entra<-0.5*M_PI)entra+= M_PI;
            while(entra> 0.5*M_PI)entra -= M_PI;
            //cout<<"get 3"<<endl;
            double zhit     = hitPoint.z();
            double fltLen   = trkIT->flightLength(hitPoint);
            //cout<<"after get "<<hitId<<endl;
 
            recMdcHit->setDriftDistLeft(ddl);
            recMdcHit->setDriftDistRight(ddr);
            //recMdcHit->setErrDriftDistLeft(erddl);
            //recMdcHit->setErrDriftDistRight(erddr);
            recMdcHit->setChisqAdd(pChisq);
            recMdcHit->setFlagLR(lr);
            recMdcHit->setStat(stat);
            recMdcHit->setMdcId(mdcId);
            recMdcHit->setTdc(tdc);
            recMdcHit->setAdc(adc);
            recMdcHit->setDriftT(driftT);
            recMdcHit->setDoca(doca);
            recMdcHit->setEntra(entra);
            recMdcHit->setZhit(zhit);
            recMdcHit->setFltLen(fltLen);
            //cout<<"after set "<<hitId<<endl;
            // <--- make a RecMdcHit from a HoughHit  
 
            hitList->push_back(recMdcHit);
            SmartRef<RecMdcHit> refHit(recMdcHit);
            hitRefVec.push_back(refHit);
 
            if(hot->getLayer()>maxLayerId){
                maxLayerId = hot->getLayer();
                fiTermHot = hot;
            }
            if(hot->getLayer()<minLayerId){
                minLayerId = hot->getLayer();
            }
            //if(hot->getUse()==1){
            recMdcHit->setStat(1);
            //}else{
            //  recMdcHit->setStat(0);
            //}
            //cout<<"finish a DC hit "<<hitId<<endl;
            hitId++;
        }*/
        if(hot->getHitType()==HoughHit::CGEMHIT)
        {
            //hot->print();
            //cout<<endl;
            const RecCgemCluster* recCgemCluster = hot->getCgemCluster();
            int clusterid = hot->getCgemCluster()->getclusterid();
            //int stat = hot->getUse();
            int stat = 1;
            clusterFitStat[clusterid] = stat;
            clusterRefVec.push_back(recCgemCluster);
            if(debug) cout<<" cluster "<<clusterid<<" kept! "<<endl;
            if(hot->getLayer()>maxLayerId)
            {
                maxLayerId = hot->getLayer();
                fiTermHot = hot;
            }
            if(hot->getLayer()<minLayerId)
            {
                minLayerId = hot->getLayer();
            }
        }
        }// loop hotList (vector<HoughHit*>) of a hough track
 
        // --- loop RecMdcHitVec in Hough track
        double fltLen = -1;
        vector<RecMdcHit>* RecMdcHitVecPnt = trkIT->getRecMdcHitVec();
        int nMdcHits=RecMdcHitVecPnt->size();
        int nMdcHitsKept=0;
        vector<RecMdcHit>::iterator iter_recMdcHit = RecMdcHitVecPnt->begin();
        for(; iter_recMdcHit!=RecMdcHitVecPnt->end(); iter_recMdcHit++)
        {
            RecMdcHit* recMdcHit = new RecMdcHit(*iter_recMdcHit);
            if(recMdcHit->getChisqAdd()>m_chi2CutHits) // skip hit with chi2>m_chi2CutHits
            {
                //cout<<"skip hit with chi2="<<recMdcHit->getChisqAdd()<<endl;
                delete recMdcHit;
                continue;
            }
            recMdcHit->setId(hitId);
            recMdcHit->setTrkId(trackId);
            recMdcHit->setStat(1);
            hitList->push_back(recMdcHit);
            SmartRef<RecMdcHit> refHit(recMdcHit);
            hitRefVec.push_back(refHit);
            nMdcHitsKept++;
 
            Identifier mdcid = recMdcHit->getMdcId();
            int layer = MdcID::layer(mdcid);
            int wire  = MdcID::wire(mdcid);
            if(debug) cout<<" hit at layer "<<layer<<" wire "<<wire<<" kept! "<<endl;
            if(layer>maxLayerId)
            {
                maxLayerId = layer;
                fltLen=recMdcHit->getFltLen();
                //fiTermHot = hot;
            }
            if(layer<minLayerId)
            {
                minLayerId = layer;
            }
            hitId++;
        }
        if(debug) cout<<"track "<<trackId<<", "<<nMdcHitsKept<<"/"<<nMdcHits<<" hits kept"<<endl;
 
        //fi terminal angle
        if (fiTermHot!=NULL){
            //fiTerm=(1./sqrt(1.+tanDip*tanDip))*fiTermHot->fltLen()*helix.omega();
            HepPoint3D point = fiTermHot->getHitPosition();
            fiTerm = trkIT->flightArc(point)/trkIT->radius();
        }
        if(fltLen>0) fiTerm=-fltLen*sin(theta)/trkIT->radius();
        recMdcTrack->setFiTerm(fiTerm);
 
        // --- setN* functions called in setVec* functions
        recMdcTrack->setVecHits(hitRefVec);
        std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
        recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
        trackList->push_back(recMdcTrack);
        trackId++;
    }// loop HoughTrack list
    return trackList->size();
 
}
 
int HoughFinder::storeTracks(RecMdcTrackCol* trackList, RecMdcHitCol* hitList, vector<HoughTrack>& trackVector){
    bool debug=false; //debug=true;
    int trackId(0);
    int tkStat =4;
    if(debug) cout<<"N Rectrack: "<<trackVector.size()<<endl;
    sort(trackVector.begin(),trackVector.end(),largerPt);
 
    for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
        if(debug) cout<<"trk flag: "<<(trkIT)->getFlag()<<endl;
        if((trkIT)->getFlag() == 1)continue;
        if((trkIT)->getFlag() == -1)continue;
        if((trkIT)->getFlag() == -2)continue;
        if(trkIT->getTrkRecoTrk()==NULL)
        {
            if(debug) cout<<"trk getTrkRecoTrk NULL!"<<endl;
            continue;
        }
        //get the results of the fit to this track
        const TrkFit* fitresult = trkIT->getTrkRecoTrk()->fitResult();
        //check the fit worked
        if (fitresult == 0) 
        {
            if(debug) cout<<"no fit result!"<<endl;
            continue;
        }
        if(debug) 
        {
            cout<<"the following track will be saved: "<<endl;
            trkIT->print();
        }
        if(printFlag)trkIT->print();
        if(printFlag)trkIT->printHot();
        //cout<<endl;
 
        //new a Rec. Track MdcTrack 
        RecMdcTrack* recMdcTrack = new RecMdcTrack();
        const TrkHitList* aList = trkIT->getTrkRecoTrk()->hits();
        //some track info
        const BField& theField = trkIT->getTrkRecoTrk()->bField();
        double Bz = theField.bFieldZ();
        //Fit info
        int hitId = 0;
        int nHits=0;
        int nSt=0;
        //int nAct=0; int nSeg=0; 
        //int maxlay = 0; int minlay = 43; int seg[11];/* 11 slayer */ int segme;
        //****************************
        //* active flag open this
        //****************************
        //* if (allHit>0){
        //*   nHits = aList->nHit();//add inActive hit also
        //* }else{
        //*   nHits = fitresult->nMdc();// = nActive()
        //* }
        //****************************
        //for 5.1.0 use all hits
        nHits = aList->nHit();
        //****************************
        // use active only
        // nHits = fitresult->nMdc();// = nActive()
        //****************************
 
        int    q        = fitresult->charge();
        double chisq    = fitresult->chisq();
        int    nDof     = fitresult->nDof();
 
        //if(nHits-nDof!=5) cout<<"HoughFinder::storeTracks: nHits-nDof!=5, nHits="<<nHits<<" , nDof="<<nDof<<", chi2="<<chisq<<" , evt="<<m_event<<endl;
 
        //track parameters at closest approach to beamline
        double fltLenPoca = 0.0;
        TrkExchangePar helix = fitresult->helix(fltLenPoca);
        //std::cout<< __FILE__ << " phi0  " << helix.phi0() << " omega  " <<helix.omega()<<std::endl;
        double phi0 = helix.phi0();
        double tanDip = helix.tanDip();
        //double z0 = helix.z0();
        double d0 = helix.d0();
        //momenta and position
        //CLHEP::Hep3Vector mom = fitresult->momentum(fltLenPoca);
        HepPoint3D poca = fitresult->position(fltLenPoca);
 
        double helixPar[5];
        //dro =-d0
        helixPar[0] = -d0;
        //phi0 = phi0 - pi/2   [0,2pi)
        double tphi0 = phi0 - Constants::pi/2.;
        if (tphi0 < 0. ) tphi0 += Constants::pi * 2.;
        helixPar[1] = tphi0;
        //kappa = q/pxy
        double pxy = fitresult->pt();
        if (pxy == 0.)  helixPar[2] = 9999.;
        else helixPar[2] = q/fabs(pxy);
        if(pxy>9999.) helixPar[2] = 0.00001;
        //dz = z0
        helixPar[3] = helix.z0();
        //tanl 
        helixPar[4] = tanDip;
        //error
        //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , helix.covariance() = V(X)
        HepSymMatrix mS(helix.params().num_row(),0);
        mS[0][0]=-1.;//dr0=-d0
        mS[1][1]=1.;
        mS[2][2]=-333.567/Bz;//pxy -> cpa
        mS[3][3]=1.;
        mS[4][4]=1.;
        HepSymMatrix mVy = helix.covariance().similarity(mS);
        double errorMat[15];
        int k = 0;
        for (int ie = 0 ; ie < 5 ; ie ++){
            for (int je = ie ; je < 5 ; je ++) {
                errorMat[k] = mVy[ie][je];
                k++;
            }
        }
        double p,px,py,pz;
        px = pxy * (-sin(helixPar[1]));
        py = pxy * cos(helixPar[1]);
        pz = pxy * helixPar[4];
        p = sqrt(pxy*pxy + pz*pz);
        //theta, phi
        double theta = acos(pz/p);
        double phi = atan2(py,px);
        recMdcTrack->setTrackId(trackId);
        recMdcTrack->setHelix(helixPar);
        recMdcTrack->setCharge(q);
        recMdcTrack->setPxy(pxy);
        recMdcTrack->setPx(px);
        recMdcTrack->setPy(py);
        recMdcTrack->setPz(pz);
        recMdcTrack->setP(p);
        recMdcTrack->setTheta(theta);
        recMdcTrack->setPhi(phi);
        recMdcTrack->setPoca(poca);
        recMdcTrack->setX(poca.x());//poca
        recMdcTrack->setY(poca.y());
        recMdcTrack->setZ(poca.z());
        recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
        HepPoint3D pivot(0.,0.,0.);
        recMdcTrack->setPivot(pivot);
        recMdcTrack->setVX0(0.);//pivot 
        recMdcTrack->setVY0(0.);
        recMdcTrack->setVZ0(0.);
        recMdcTrack->setError(errorMat);
        recMdcTrack->setError(mVy);
        recMdcTrack->setChi2(chisq);
        recMdcTrack->setStat(tkStat);
        //recMdcTrack->setNhits(nHits);
        //-----------hit list-------------
        HitRefVec  hitRefVec;
        ClusterRefVec clusterRefVec;
        vector<int> vecHits;
        map<int,int> clusterFitStat;
        //for fiterm
        int maxLayer = 0;
        //cout<<__FILE__ <<"  "<<__LINE__ <<endl;
        for(TrkHitList::hot_iterator hot = aList->begin();hot!=aList->end();hot++){
            if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
                const MdcRecoHitOnTrack* recoHot;             
                recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
                int layerId = recoHot->mdcHit()->layernumber();
                int wireId = recoHot->mdcHit()->wirenumber();
                if(maxLayer < layerId)maxLayer = layerId;
            }
        }
        //cout<<"maxLayer:"<<maxLayer<<endl;
        int maxLayerId = -1;
        int minLayerId = 43;
        double fiTerm = 999.;
        const MdcRecoHitOnTrack* fiTermHot = NULL;
        TrkHitList::hot_iterator hot = aList->begin();
        for (;hot!=aList->end();hot++){
            //cout<<__FILE__ <<"  "<<__LINE__ <<"  "<< typeid(*hot).name()<<endl;
            if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
                const MdcRecoHitOnTrack* recoHot;             
                recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
                if(!(recoHot->isActive())) continue;// added 2020-05-26 wangll
                int layerId = recoHot->mdcHit()->layernumber();
                if(layerId>(maxLayer - m_removeNOuterHits))continue;
                //if (!recoHot->mdcHit()) return;
                RecMdcHit* recMdcHit = new RecMdcHit;
                //id
                recMdcHit->setId(hitId);
                //---------BESIII----------
                //phiWire - phiHit <0 flagLR=0 left driftleft<0 ;
                //phiWire - phiHit >0 flagLR=1 right driftright>0;
                //flag = 2 ambig;
                //-----Babar wireAmbig----
                //-1 -> right, 0 -> don't know, +1 -> left
                //+1 phiWire-phiHit<0 Left,
                //-1 phiWire-phiHit>0 Right, 
                //0  don't know   
                //ambig() w.r.t track direction
                //wireAmbig() w.r.t. wire location 
                double hotWireAmbig = recoHot->wireAmbig();
                double driftDist = fabs(recoHot->drift());
                double sigma = recoHot->hitRms();       
                double doca = fabs(recoHot->dcaToWire());
                int flagLR=2;
                if ( hotWireAmbig == 1){
                    flagLR = 0; //left driftDist <0
                    doca *= -1.;//2012-07-19 
                }else if( hotWireAmbig == -1){
                    flagLR = 1; //right driftDist >0
                }else if( hotWireAmbig == 0){
                    flagLR = 2; //don't know 
                }
                recMdcHit->setFlagLR(flagLR);
                recMdcHit->setDriftDistLeft((-1)*driftDist);
                recMdcHit->setErrDriftDistLeft(sigma);
                recMdcHit->setDriftDistRight(driftDist);
                recMdcHit->setErrDriftDistRight(sigma);
                //Mdc Id
                Identifier mdcId = recoHot->mdcHit()->digi()->identify();
                recMdcHit->setMdcId(mdcId);
                //corrected ADC
 
                //contribution to chisquare
                //fill chisq
                double res=999.,rese=999.;
                if (recoHot->resid(res,rese,false)){
                }else{}
                double deltaChi=0;
                recoHot->getFitStuff(deltaChi);//yzhang open 2010-09-20 
                recMdcHit->setChisqAdd( deltaChi * deltaChi );
                //set tracks containing this hit
                recMdcHit->setTrkId(trackId);
                //doca
                recMdcHit->setDoca(doca);//doca sign left<0
                //entrance angle
 
                recMdcHit->setEntra(recoHot->entranceAngle());
                //z of hit 
                HepPoint3D pos; Hep3Vector dir;
                double fltLen = recoHot->fltLen(); 
                recMdcHit->setFltLen(fltLen);
                recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
                recMdcHit->setZhit(pos.z());
                double tof = recoHot->getParentRep()->arrivalTime(fltLen);
                recMdcHit->setTdc(recoHot->mdcHit()->tdcIndex());
                recMdcHit->setAdc(recoHot->mdcHit()->adcIndex()); 
                //drift time      
                recMdcHit->setDriftT(recoHot->mdcHit()->driftTime(tof,pos.z()));
                //for fiterm
                int wireId = recoHot->mdcHit()->wirenumber();
                //std::cout<<" MdcTrack::store() ("<<layerId<<","<<wireId<<") fltLen "<<fltLen<<" wireAmbig  "<<hotWireAmbig<<"    y "<<pos.y()<<std::endl;
                //  <<recoHot->entranceAngle()<<std::endl;
                if (layerId >= maxLayerId){
                    maxLayerId = layerId;
                    fiTermHot = recoHot;
                } 
                if (layerId < minLayerId){
                    minLayerId = layerId;
                }
                // status flag
                if (recoHot->isActive()) {
                    recMdcHit->setStat(1);
                }else{
                    recMdcHit->setStat(0);    
                }
                // for 5.1.0 use all hits
                if (recoHot->layer()->view()) {
                    ++nSt;
                }
                hitList->push_back(recMdcHit);
                SmartRef<RecMdcHit> refHit(recMdcHit);
                hitRefVec.push_back(refHit);        
                vecHits.push_back(mdcId.get_value());
                ++hitId;          
            }else if(typeid(*hot)==typeid(CgemHitOnTrack)){
                const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));            
                const RecCgemCluster* recCgemCluster = recoHot->baseHit();
                int clusterid = recCgemCluster->getclusterid();
                int stat = recoHot->isActive();
                if(stat==0) continue;// added 2020-05-26 wangll
                //cout<<"stat:"<<stat<<endl;
                clusterFitStat[clusterid] = stat;
                //cout<<clusterid<<"  "<<stat<<endl;
                //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
                clusterRefVec.push_back(recCgemCluster);
            }
        }
        std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
        //fi terminal angle
        if (fiTermHot!=NULL){
            fiTerm=(1./sqrt(1.+tanDip*tanDip))*fiTermHot->fltLen()*helix.omega();
        }
        recMdcTrack->setFiTerm(fiTerm);
        // number of  stereo hits contained
        //recMdcTrack->setNhits(hitId);
        //recMdcTrack->setNster(nSt);  
        //recMdcTrack->setNdof(nDof);
        //recMdcTrack->setFirstLayer(maxLayerId);  
        //recMdcTrack->setLastLayer(minLayerId);  
        //recMdcTrack->setFirstLayer(minLayerId);  
        //recMdcTrack->setLastLayer(maxLayerId);  
        //recMdcTrack->setClusterFitStat(clusterFitStat);
        //recMdcTrack->setVecClusters(clusterRefVec);
        //recMdcTrack->setNcluster(clusterRefVec.size());
 
        // --- setN* functions called in setVec* functions
        recMdcTrack->setVecHits(hitRefVec);
        recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
        trackList->push_back(recMdcTrack);
        trackId++;
    }
    //cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
    /*
       int trackId(0);
       for(vector<HoughTrack>::iterator trkIT = trackVector.begin();trkIT!=trackVector.end();trkIT++){
//TrkErrCode trkErrCode = trkIT->fitCircle(m_mdcDetector,m_trkContextEv,m_bunchT0);
HoughTrack* houghTrack = &(*trkIT);
if(houghTrack->getFlag() == 1)continue;
if(houghTrack->getFlag() < 0)continue;
//houghTrack->print();
//houghTrack->printHot();
trackId=trackList->size();
//cout<<"start store trk "<<trackId;
RecMdcTrack* recMdcTrack = new RecMdcTrack();
//vector<HoughHit>* hotList = houghTrack->getHot();
//vector<HoughHit*> hotList = houghTrack->getVecHitPnt();
vector<HoughHit*> hotList = houghTrack->getHotList();
//cout<<"  with "<<hotList.size()<<" hits";//<<endl;
double helixPar[5];
helixPar[0]  = houghTrack->dr();
helixPar[1]  = houghTrack->phi0();
helixPar[2]  = houghTrack->kappa();
helixPar[3]  = houghTrack->dz();
helixPar[4]  = houghTrack->tanl();
int    q     = houghTrack->getCharge();
double pxy   = houghTrack->pt();
double px    = houghTrack->momentum(0).x();
double py    = houghTrack->momentum(0).y();
double pz    = houghTrack->momentum(0).z();
double p     = houghTrack->momentum(0).mag();
double theta = houghTrack->direction(0).theta();
double phi   = houghTrack->direction(0).phi();
HepPoint3D poca  = houghTrack->x(0);
HepPoint3D pivot = houghTrack->pivot();
double r = poca.perp();
HepSymMatrix Ea = houghTrack->Ea();
//cout<<"  phi, pxy, pz = "<<phi<<", "<<pxy<<", "<<pz<<endl;
double errorMat[15];
int k = 0;
for (int ie = 0 ; ie < 5 ; ie ++){
for (int je = ie ; je < 5 ; je ++){
errorMat[k] = Ea[ie][je];
k++;
}
}
double chisq = houghTrack->getChi2();
int    tkStat = 4+houghTrack->getCircleFitStat();
//if(tkStat==4) cout<<"    track "<<trackId<<" has no good circle fit"<<endl;
int    nHits = hotList.size();
 
recMdcTrack->setTrackId(trackId);
recMdcTrack->setHelix(helixPar); 
recMdcTrack->setCharge(q);
recMdcTrack->setPxy(pxy);
recMdcTrack->setPx(px);
recMdcTrack->setPy(py);
recMdcTrack->setPz(pz);
recMdcTrack->setP(p);
recMdcTrack->setTheta(theta);
recMdcTrack->setPhi(phi);
recMdcTrack->setPoca(poca);
recMdcTrack->setX(poca.x());//poca
recMdcTrack->setY(poca.y());
recMdcTrack->setZ(poca.z());
recMdcTrack->setR(r);
recMdcTrack->setPivot(pivot);
recMdcTrack->setVX0(0.);//pivot 
recMdcTrack->setVY0(0.);
recMdcTrack->setVZ0(0.);
recMdcTrack->setError(errorMat);
recMdcTrack->setError(Ea);
recMdcTrack->setChi2(chisq);
recMdcTrack->setStat(tkStat);
recMdcTrack->setNhits(nHits);
//-----------hit list-------------
int hitId = 0;
int nSt(0);
int maxLayerId = -1;
int minLayerId = 43;
double fiTerm = 999.;
HoughHit fiTermHot;
 
HitRefVec  hitRefVec;
ClusterRefVec clusterRefVec;
map<int,int> clusterFitStat;
for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end();hotIt++)
{
    //cout<<"handle hit "<<hitId<<endl;
    HoughHit* hot(*hotIt);
    //int layer = (**hotIt).getLayer();
    //int wire  = (**hotIt).getWire();
    //cout<<"("<<setw(2)<<layer<<","<<setw(3)<<wire<<") ";
    if(hot->getHitType()==HoughHit::MDCHIT){
        //cout<<"create a new RecMdcHit"<<endl;
        RecMdcHit* recMdcHit = new RecMdcHit;
 
        recMdcHit->setId(hitId);
        recMdcHit->setTrkId(trackId);
        HepPoint3D hitPoint = hot->getHitPosition();
        double dPhi = houghTrack->dPhi(hitPoint);
        HepPoint3D position  = houghTrack->x(dPhi);
        Hep3Vector direction = houghTrack->direction(dPhi);
        //cout<<"get 1"<<endl;
 
        double ddl      = -1*hot->getDriftDist();
        double ddr      = hot->getDriftDist();
        //double erddl    = hot->hitSigma();
        //double erddr    = hot->hitSigma();
        double pChisq   = hot->getResidual()*hot->getResidual();//FIXME
        int    lr       = 0;
        int    stat     = hot->getUse();
        stat = 1; // llwang
        //cout<<"hit stat = "<<stat<<endl;
        //cout<<"get 2"<<endl;
        Identifier mdcId = hot->getDigi()->identify();
        //cout<<"get 2.1"<<endl;
        double tdc      = hot->getDigi()->getTimeChannel();
        //cout<<"get 2.2"<<endl;
        double adc      = hot->getDigi()->getChargeChannel();
        //cout<<"get 2.3"<<endl;
        double driftT   = hot->driftTime();
        //cout<<"get 2.4"<<endl;
        double doca     = hot->getDriftDist() +hot->getResidual();
        //cout<<"get 2.5"<<endl;
        double entra    = direction.phi()-position.phi();//FIXME
        //cout<<"get 2.6"<<endl;
        while(entra<-0.5*M_PI)entra+= M_PI;
        while(entra> 0.5*M_PI)entra -= M_PI;
        //cout<<"get 3"<<endl;
        double zhit     = hitPoint.z();
        double fltLen   = houghTrack->flightLength(hitPoint);
        //cout<<"after get "<<hitId<<endl;
 
        recMdcHit->setDriftDistLeft(ddl);
        recMdcHit->setDriftDistRight(ddr);
        //recMdcHit->setErrDriftDistLeft(erddl);
        //recMdcHit->setErrDriftDistRight(erddr);
        recMdcHit->setChisqAdd(pChisq);
        recMdcHit->setFlagLR(lr);
        recMdcHit->setStat(stat);
        recMdcHit->setMdcId(mdcId);
        recMdcHit->setTdc(tdc);
        recMdcHit->setAdc(adc);
        recMdcHit->setDriftT(driftT);
        recMdcHit->setDoca(doca);
        recMdcHit->setEntra(entra);
        recMdcHit->setZhit(zhit);
        recMdcHit->setFltLen(fltLen);
        //cout<<"after set "<<hitId<<endl;
 
        hitList->push_back(recMdcHit);
        SmartRef<RecMdcHit> refHit(recMdcHit);
        hitRefVec.push_back(refHit);
 
        if(hot->getLayer()>maxLayerId){
            maxLayerId = hot->getLayer();
            fiTermHot = *hot;
        }
        if(hot->getLayer()<minLayerId){
            minLayerId = hot->getLayer();
        }
        if(hot->getUse()==1){
            recMdcHit->setStat(1);
        }else{
            recMdcHit->setStat(0);
        }
        if(hot->getFlag()!=0)++nSt;
        //cout<<"finish a DC hit "<<hitId<<endl;
    }else if(hot->getHitType()==HoughHit::CGEMHIT){
        //hot->print();
        //cout<<endl;
        const RecCgemCluster* recCgemCluster = hot->getCgemCluster();
        int clusterid = hot->getCgemCluster()->getclusterid();
        //int stat = hot->getUse();
        int stat = 1;
        clusterFitStat[clusterid] = stat;
        clusterRefVec.push_back(recCgemCluster);
    }
    hitId++;
}
 
//HepPoint3D point = fiTermHot.getHitPosition();
//fiTerm = houghTrack->flightArc(point)/houghTrack->radius();
//recMdcTrack->setFiTerm(fiTerm);
recMdcTrack->setNdof(nHits-5);
recMdcTrack->setNster(nSt);
recMdcTrack->setFirstLayer(minLayerId);
recMdcTrack->setLastLayer(maxLayerId);  
//cout<<"push back a RecMdcTrack with "<<recMdcTrack->getNhits()<<" hits, "<<recMdcTrack->getNcluster()<<" clusters"<<endl;
recMdcTrack->setVecHits(hitRefVec);
//cout<<"before setVecClusters"<<endl;
//recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
//recMdcTrack->setVecClusters(clusterRefVec);
recMdcTrack->setVecClusters(clusterRefVec,clusterFitStat);
//cout<<"clusterRefVec.size = "<<clusterRefVec.size()<<endl;
//cout<<"end setVecClusters"<<endl;
recMdcTrack->setNcluster(clusterRefVec.size());
//cout<<"push back a RecMdcTrack with "<<recMdcTrack->getNhits()<<" hits, "<<recMdcTrack->getNcluster()<<" clusters"<<endl;
trackList->push_back(recMdcTrack);
//trackId++;
//cout<<helixPar[2]<<endl;
}
return trackList->size();
// */
return trackList->size();
}
 
void HoughFinder::printHitList()
{
    //cout<<"nHit ==========> "<<m_houghHitList.size()<<endl;
    cout<<"========== truth =========="<<endl;
    for(HitVector_Iterator hitIt=m_mcHitCol.begin();hitIt!=m_mcHitCol.end();hitIt++){
        hitIt->print();
    }
    cout<<endl;
    cout<<"========== Digi =========="<<endl;
    for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end();hitIt++)
    {
        hitIt->print();
    }
    cout<<endl;
    /*
       {
       vector<HoughHit>& hitList = m_mcHitCol;
       for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
       cout<<setw(4)<<hitIt->getHitID();
       cout<<"("<<setw(2)<<hitIt->getLayer()<<", "<<setw(3)<<hitIt->getWire()<<")    ";
       cout<<"("<<setw(10)<<hitIt->getHitPosition().x()<<", "<<setw(10)<<hitIt->getHitPosition().y()<<", "<<setw(10)<<hitIt->getHitPosition().z()<<")    ";
       vector<int> trkID = hitIt->getTrkID();
       cout<<"trkID:"<<setw(4)<<trkID[0];
       cout<<"half:"<<setw(4)<<hitIt->getHalfCircle();
       if(hitIt->getHitType()==HoughHit::CGEMMCHIT){
       if(hitIt->getCgemCluster()!=NULL){
    //cout<<hitIt->getFlag()<<"    ";
    //if(hitIt->getFlag()==2)
    cout<<setw(5)<<hitIt->getCgemCluster()->getclusterid();
    cout<<"["<<setw(3)<<hitIt->getCgemCluster()->getclusterflagb()<<", "<<setw(3)<<hitIt->getCgemCluster()->getclusterflage()<<"]    ";
    }
    //else cout<<"NULL  ";
    }
    else{
    if(hitIt->getDigi()!=NULL){
    cout<<setw(5)<<hitIt->getDigi()->getTrackIndex();
    cout<<hitIt->getDigi()->identify();
    }
    //else cout<<"NULL  ";
 
    }
    cout<<endl;
    }
    }
    cout<<endl;
    //*/
    /*
       {
       vector<HoughHit>& hitList = m_houghHitList;
       cout<<hitList.size()<<endl;
       for(HitVector_Iterator hitIt = hitList.begin(); hitIt != hitList.end();hitIt++){
       cout<<setw(4)<<hitIt->getHitID();
       cout<<"("<<setw(2)<<hitIt->getLayer()<<", "<<setw(3)<<hitIt->getWire()<<")    ";
       cout<<"("<<setw(10)<<hitIt->getHitPosition().x()<<", "<<setw(10)<<hitIt->getHitPosition().y()<<", "<<setw(10)<<hitIt->getHitPosition().z()<<")    ";
       cout<<setw(4)<<hitIt->getFlag();
       if(hitIt->getHitType()==HoughHit::CGEMHIT){
       cout<<setw(4)<<hitIt->getCgemCluster()->getclusterid();
       cout<<"["<<setw(4)<<hitIt->getCgemCluster()->getclusterflagb()<<", "<<setw(4)<<hitIt->getCgemCluster()->getclusterflage()<<"]    ";
       }
       else{
       cout<<setw(5)<<hitIt->getDigi()->getTrackIndex();
       cout<<hitIt->getDigi()->identify()<<"     ";
    //cout<<"("<<hitIt->getPairHit()->getLayer()<<", "<<hitIt->getPairHit()->getWire()<<")    ";
    //if(hitIt->getPairHit()!=NULL)cout<<hitIt->getDriftDist()<<"     "<<hitIt->getPairHit()->getDriftDist()<<"     "<<hitIt->getDriftDist()-hitIt->getPairHit()->getDriftDist()<<endl;
    }
    if(hitIt->getPairHit()!=NULL)cout<<hitIt->getPairHit()->getHitID();
    //else cout<<"NULL";
    cout<<endl;
    }
    }
    cout<<endl;
    // */
}
 
 
int HoughFinder::searchCircle()
{
    //m_debug=false; 
    //m_debug=true;
    if(m_debug==5) {
        cout<<"***************************************"<<endl;
        cout<<"start circle search : "<<endl;
    }
    m_triedCellInRoughMap.clear();
    bool tryCgemLeftOnly=false;
    while(true) 
    {
        int charge=0;
        // --- fill rough map
        int nXhitsLeft=0;
        m_roughRhoThetaMap.Reset();
        //int charge=0; 
        int vote=1;
 
        //cout<<"searchCircle "<<__LINE__<<endl;
        nXhitsLeft = fillHistogram(m_XHoughHitList,&m_roughRhoThetaMap,charge,vote);
        //cout<<"searchCircle "<<__LINE__<<endl;
        if(m_debug==5) {
            cout<<"------------------------------------------------------------------------------------ "<<endl;
            cout<<"nXhitsLeft = "<<nXhitsLeft<<endl;
        }
 
        if(nXhitsLeft<3) {
            if(!tryCgemLeftOnly) {
                int nCgemHitsLeft = activeUnusedCgemHitsOnly(m_XHoughHitList);
                tryCgemLeftOnly=true;
                if(nCgemHitsLeft>=3) {
                    if(m_debug==5) cout<<"Now try CGEM left only"<<endl;
                    continue;
                }
            }
            break;
        }
 
        // --- find the peak position on rough Hough map
        double x_peak(0.), y_peak(0.);
        double x_weight(0.), y_weight(0.);
        getWeightedPeak(m_roughRhoThetaMap, x_peak, y_peak, x_weight, y_weight);
        if(m_debug==5) cout<<"rough x_weight, y_weight = "<<x_weight<<", "<<y_weight<<endl;
        //cout<<endl;
 
        // --- range for fine map
        x_peak = x_weight;
        y_peak = y_weight;
        double x_min=x_peak-0.5*m_ExtPeak_theta*M_PI/m_nBinTheta;
        double x_max=x_peak+0.5*m_ExtPeak_theta*M_PI/m_nBinTheta;
        double y_min=y_peak-0.5*m_ExtPeak_rho*0.2/m_nBinRho;
        double y_max=y_peak+0.5*m_ExtPeak_rho*0.2/m_nBinRho;
 
        // --- binning for fine map
        int nFineBin_Theta = nFineBinTheta(y_peak);
        int nFineBin_Rho   = nFineBinRho(y_peak);
        if(m_debug==5) cout<<"nThetaFine, nRhoFine = "<<nFineBin_Theta<<", "<<nFineBin_Rho<<endl;
 
        int trkFoundThisTry = 0;
        int NHitTried = 0;
        // --- loop -/+ hypotheses
        for(charge=-1; charge<2; charge+=2)
        {
 
            if(m_debug==5) cout<<"/********    Try iCharge = "<<charge<<"  **********/"<<endl;
 
            // create a rough track candidate
            int trkID = m_houghTrackList.size();
            double binTheta(0.), binRho(0.);
            HoughTrack track(charge, x_peak, y_peak, binTheta, binRho, trkID);
            if(m_debug==5) cout<<"circle center = "<<track.center()<<endl;
            // --- fill the fine map
            m_fineRhoThetaMap.Reset();
            m_fineRhoThetaMap.SetBins(nFineBin_Theta,x_min,x_max,nFineBin_Rho,y_min,y_max);
            vote=1;
            int nHitFilled(0);
            //track.print();
            if(fabs(y_weight)<0.02) // > 300 MeV
            {
                //cout<<"fill high pt histogram"<<endl;
                nHitFilled = fillHistogram(m_XHoughHitList,&m_fineRhoThetaMap,0,vote);
            }
            else
            {
                //cout<<"fill low  pt histogram"<<endl;
                nHitFilled = fillHistogram(m_XHoughHitList,&m_fineRhoThetaMap,charge,vote);
            }
            if(nHitFilled<3) {
                if(m_debug==5) cout<<"nHitFilled<3 -> continue"<<endl;
                //delete track;
                continue;
            }
            double theta_peak, rho_peak, theta_weight,  rho_weight;
            getWeightedPeak(m_fineRhoThetaMap, theta_peak, rho_peak, theta_weight,  rho_weight);
 
            // update the track candidate
            track.update(theta_weight,rho_weight);
            //cout<<" ============== before fitting ================="<<endl;
            //cout<<"fine  x_weight, y_weight = "<<theta_weight<<", "<<rho_weight<<endl;
            //track.print();
 
            // --- select X-view hits
            //vector<HoughHit*> XHitList_sel;
            int flag=0;// only X-view hits
            track.resetNhitHalf();
            int nXsel(0);
            if(fabs(y_weight)<0.02) // > 300 MeV
            { 
                track.findXHot(m_XHoughHitList,0);
                //if(track->getNhitFirstHalf()<track->getNhitSecondHalf()) 
                if(charge==-1&&track.getNhitUnusedFirstHalf()<track.getNhitUnusedSecondHalf()) 
                {
                    if(m_debug==5) cout<<"--- the second half has more hits, try next charge: "<<endl;
                    //delete track;
                    continue;
                }
                else 
                {
                    track.clearHits();
                    //XHitList_sel.clear();
                    if(m_debug==5) cout<<" --- the first  half has more hits, go ahead:  "<<endl;
                }
            }
            nXsel = track.findXHot(m_XHoughHitList, charge);
            //track.printHot();
            NHitTried = track.getNTried();
            if(m_debug==5) {
                cout<<"sel "<<nXsel<<" hits, "<<NHitTried<<" new hits tried"<<endl;
                track.printHot(0);
            }
 
            //if(track.isAGoodCircle()) 
            if(nXsel>=3) 
            {
                //track.printHot();
                //track.print();
                //TrkContextEv* trkContextEv = new TrkContextEv(m_bfield);
                //if(m_fitFlag>0)TrkErrCode trkErrCode = track.fitCircle(m_mdcDetector,trkContextEv,m_bunchT0);
                int fitFlag_circle = 0;
                if(m_fitFlag>0) 
                {
                    //TrkErrCode trkErrCode = track.fitCircle(m_mdcDetector,m_trkContextEv,m_bunchT0);
                    //if(m_debug) 
                    //cout<<"circle helix from TrkFitter"<<track.a()<<", chi2="<<track.getChi2()<<endl;
                    //cout<<"m_bunchT0 = "<<m_bunchT0<<endl;
                    //track.update(theta_weight,rho_weight);
                    //cout<<"to start circle fit"<<endl;
                    fitFlag_circle = track.fitCircle(&myDotsHelixFitter,m_bunchT0*1e9, m_circle_chiCut, m_debug);// FIXME, try track from IP first, track with ODC hits only?
                    //cout<<"circle fit ends"<<endl;
                    if(fitFlag_circle==0) 
                    {
                        HepVector a_circle_fit = myDotsHelixFitter.getHelix();
                        track.updateCirclePar(a_circle_fit[0], a_circle_fit[1], a_circle_fit[2]);
                        track.clearHits();
                        int nXsel_new = track.findXHot(m_XHoughHitList, charge);// too loose? FIXME
                        int nDrop = track.dropRedundentCgemXHits();// or tag them only? FIXME
                        if(m_debug==5) {
                            cout<<"   sel "<<nXsel_new<<" hits after fitCircle, fitted pt = "<<track.pt();
                            if(nDrop>0) cout<<", drop "<<nDrop<<" redundent CgemXHits";
                            cout<<endl;
                            track.printHot(0);
                        }
                        NHitTried = track.getNTried();
                        if(nXsel_new==0||NHitTried==0) {
                            track.update(theta_weight,rho_weight);
                            nXsel = track.findXHot(m_XHoughHitList, charge);
                            NHitTried = track.getNTried();
                        }
                        //cout<<" ============== after fitting ================="<<endl;
                        //track.printHot();
                        //track.print();
                    }// good circle fit
                    //else if(m_debug) cout<<"fitCircle failed "<<endl;
                }// if do circle fit
                //track.print();
 
                //*/
                //track.markUsedHot(XHitList_sel);
                if(track.isAGoodCircle()) 
                {
                    track.markUsedHot(1);
                    //track.setMdcHit( &m_mdcHitCol);
                    m_houghTrackList.push_back(track);
                    if(m_debug==5) cout<<"save this circle trk"<<endl;
                    //track.print();
                    //track.printHot();
                    trkFoundThisTry++;
                    break;// a good circle found, search next peak
                }
                //else 
                //{
                //        if(m_debug) cout<<"drop this circle trk"<<endl;
                //        track.markUsedHot(-1);
                //}
            }// >= 3 hits
            //else 
            //{
            // --- drop the circle track
            if(m_debug==5) cout<<"drop this circle trk"<<endl;
            //track.markUsedHot(XHitList_sel, -1);
            track.markUsedHot(-1);// drop too many hits? FIXME
            //delete track;
            //}
        }// loop charge
        if(NHitTried==0) 
        {
            if(m_debug==5) cout<<"no new hits tried in this try"<<endl;
            if(!tryCgemLeftOnly) {
                int nCgemHitsLeft = activeUnusedCgemHitsOnly(m_XHoughHitList);
                tryCgemLeftOnly=true;
                if(m_debug==5) cout<<"check nCgemHitsLeft="<<nCgemHitsLeft<<endl;
                if(nCgemHitsLeft>=3) {
                    if(m_debug==5) cout<<"Now try CGEM left only"<<endl;
                    continue;
                }
            }
            if(m_debug==5) cout<<"stop circle trk finding "<<endl;
            break;
        }
    }// while
 
    //m_debug=false;
    return 1;
}
 
int HoughFinder::searchCircle2()
{
    //m_debug=true;
    if(m_debug) {
        cout<<"***************************************"<<endl;
        cout<<"start circle search : "<<endl;
    }
    m_roughRhoThetaMap.Reset();
    int charge=0;
    int nXhitsLeft=0;
    int vote=1;
    nXhitsLeft = fillHistogram(m_XHoughHitList,&m_roughRhoThetaMap,charge,vote);
 
    if(nXhitsLeft<=3) return 0;
 
    // --- find local peaks
 
}
 
void HoughFinder::getWeightedPeak(TH2D& h, double& x_peak, double& y_peak, double& x_weight, double& y_weight, int x_ext, int y_ext)
{
    int ix_max, iy_max, iz_max;
    h.GetMaximumBin(ix_max, iy_max, iz_max);
    int nx=h.GetXaxis()->GetNbins();
    int ny=h.GetYaxis()->GetNbins();
    //cout<<"peak at "<<ix_max<<", "<<iy_max<<" in a map "<<nx<<" * "<<ny<<endl;
 
    if(ix_max==0&&iy_max==0) {
        //cout<<"nx, ny = "<<nx<<", "<<ny<<endl;
        for(int i=0; i<nx; i++)
        {
            //cout<<endl;
            for(int j=0; j<ny; j++)
            {
                double n_tmp = h.GetBinContent(i,j);
                //cout<<n_tmp<<"  ";
            }
        }
    }
 
    x_peak=h.GetXaxis()->GetBinCenter(ix_max);
    y_peak=h.GetYaxis()->GetBinCenter(iy_max);
    x_weight=0.; y_weight=0.;
    double weight(0.);
    if(x_ext>=0&&y_ext>=0) {
        for(int i1=ix_max-x_ext; i1<=ix_max+x_ext; i1++)
        {
            double x_tmp = h.GetXaxis()->GetBinCenter(i1);
            int i1_tmp = i1;
            //if(i1<1) i1_tmp=m_nBinTheta+i1;
            //if(i1>1) i1_tmp=m_nBinTheta+i1;
            //if(i1==m_nBinTheta+1) i1_tmp=1;
            for(int i2=iy_max-y_ext; i2<=iy_max+y_ext; i2++)
            {
                double n_tmp = h.GetBinContent(i1_tmp,i2);
                if(i2<1||i2>ny) n_tmp=0.;
                if(i1<1||i1>nx) n_tmp=0.;
                weight+=n_tmp;
                double y_tmp = h.GetYaxis()->GetBinCenter(i2);
                x_weight+=x_tmp*n_tmp;
                y_weight+=y_tmp*n_tmp;
                //cout<<"i1,i2="<<i1<<","<<i2<<endl;
            }
        }
        x_weight=x_weight/weight;
        y_weight=y_weight/weight;
    }
    else {
        x_weight=x_peak;
        y_weight=y_peak;
    }
}
 
void HoughFinder::getLocalPeaks(TH2D& h, int minCut)
{
    int nx=h.GetXaxis()->GetNbins();
    int ny=h.GetYaxis()->GetNbins();
    vector<int> cell_locMax(nx*ny,0);
    for(int i=0; i<nx; i++) 
    {
        for(int j=0; j<ny; j++) 
        {
            int index=i+j*nx;
            if(h.GetBinContent(i+1,j+1)<minCut) continue;
            if(cell_locMax[index]<0) continue;
            for(int i1=max(i-1,0); i1<min(i+2,nx); i1++)
            {
                for(int j1=max(j-1,0); j1<min(j+2,ny); j1++)
                {
                    if(i==i1&&j==j1) continue;
                    if(h.GetBinContent(i+1,j+1)<=h.GetBinContent(i1+1,j1+1)) {
                        cell_locMax[index]=-1;
                        break;
                    }
                    else //if(h.GetBinContent(i+1,j+1)>h.GetBinContent(i1+1,j1+1))
                    {
                        cell_locMax[i1+j1*nx]=-1;
                        //cell_locMax[index]++;
                    }
                }
                if(cell_locMax[index]<0) break;
            }
            if(cell_locMax[index]==0) {
                cout<<"a local maximum found ("<<i<<", "<<j<<")"<<endl;
            }
        }
    }
}
 
int HoughFinder::nFineBinTheta(double rho)
{
    double rho_rough = fabs(rho);
    double k_theta = (600.-1200.)/0.0335;
    int nThetaFine = int(k_theta*rho_rough+1200);
    if(nThetaFine<500) nThetaFine=500;
    //cout<<"rho_rough, nThetaFine = "<<rho_rough<<", "<<nThetaFine<<endl;
    //cout<<"nThetaFine = "<<nThetaFine<<endl;
    nThetaFine=ceil(1.0*nThetaFine/m_nBinTheta*m_ExtPeak_theta);
    return nThetaFine;
}
 
int HoughFinder::nFineBinRho(double rho)
{
    double rho_rough = fabs(rho);
    int nRhoFine = 50;
    if(rho_rough<0.0198)
    {
        double k1_rho = (850.-1200.)/0.0198;
        nRhoFine=int(k1_rho*rho_rough+1200);
    }
    else if(rho_rough<0.03)
    {
        double k2_rho = (300.-850.)/(0.03-0.0198);
        nRhoFine=int(k2_rho*(rho_rough-0.0198)+850.);
    }
    else
    {
        double k3_rho = (100.-300.)/(0.056-0.03);
        nRhoFine=int(k3_rho*(rho_rough-0.03)+300.);
        if(nRhoFine<50) nRhoFine=50;
    }
    //cout<<"nRhoFine = "<<nRhoFine<<endl;
    nRhoFine=ceil(1.0*nRhoFine/m_nBinRho*m_ExtPeak_rho);
 
    return nRhoFine;
}
 
void HoughFinder::XhitCutWindow(double rho, int ilayer, double charge, double& cut1, double &cut2)
{
    rho=fabs(rho);
    if(rho>0.07) rho=0.07;
    TGraph* g_cut1, *g_cut2;
    if(ilayer<=3) {
        g_cut1=m_cut1_cgem;
        g_cut2=m_cut2_cgem;
    }
    else if(ilayer<=19) {
        g_cut1=m_cut1_ODC1;
        g_cut2=m_cut2_ODC1;
    }
    else if(ilayer<=42) {
        g_cut1=m_cut1_ODC2;
        g_cut2=m_cut2_ODC2;
    }
    // charge < 0
    cut1=g_cut1->Eval(rho);
    cut2=g_cut2->Eval(rho);
    if(charge>0) {// charge > 0
        double cut=cut1;
        cut1=-1*cut2;
        cut2=-1*cut;
    }
    else {// charge = 0
        double cut=max(fabs(cut1),fabs(cut2));
        cut1=-1*cut;
        cut2=cut;
    }
    cut1=1.1*cut1;
    cut2=1.1*cut2;
}
 
void HoughFinder::clearTrack(vector<HoughTrack>& trackVector)
{
    /*
       vector<HoughTrack*>::iterator it = trackVector.begin();
       for(; it!=trackVector.end(); it++)
       {
       HoughTrack* trk = (*it);
       delete trk;
       }
       */
    trackVector.clear();
}
 
void HoughFinder::clearMemory()
{
    //cout<<"m_mdcHitCol.size()"<<m_mdcHitCol.size()<<endl;
    vector<MdcHit*>::iterator imdcHit = m_mdcHitCol.begin();
    for(;imdcHit != m_mdcHitCol.end();imdcHit++){
        delete (*imdcHit);
    }
 
    for(vector<HoughTrack>::iterator trkIT = m_houghTrackList.begin();trkIT!=m_houghTrackList.end();trkIT++){
        trkIT->clearMemory();
        //TrkRecoTrk* trkRecoTrk = trkIT->getTrkRecoTrk();
        //delete trkRecoTrk;
    }
}
 
 
int HoughFinder::activeUnusedCgemHitsOnly(vector<HoughHit*>& hitPntList)
{
    int nCgemLeft(0);
    vector<HoughHit*>::iterator iter = hitPntList.begin();
    for(; iter!=hitPntList.end(); iter++)
    {
        //if((*iter)->getPairHit()==NULL) continue;// FIXME
        //if((*iter)->getPairHit()->getHalfCircle()!=1)continue;// FIXME
        int useOld = (*iter)->getUse();
        int iLayer = (*iter)->getLayer();
        if(iLayer<3) {
            if(useOld==-1||useOld==0) {
                (*iter)->setUse(0);
                nCgemLeft++;
            }
        }
        else 
            if(useOld==0)  (*iter)->setUse(-1);
    }
 
    return nCgemLeft;
}
 
// re-associate shard hits, llwang 2018-12-26
void HoughFinder::solveSharedHits(vector<HoughHit*>& hitList)
{
    vector<HoughHit*>::iterator itHit = hitList.begin();
    for(; itHit!=hitList.end(); itHit++)
    {
        //vector<HoughTrack*> trkPntList = it->getTrkPntVec();
        //int nTrkSharing = trkPntList.size();
        vector<int> trkIdList = (*itHit)->getTrkID();
        int nTrkSharing = trkIdList.size();
        //cout<<"nTrkSharing = "<<nTrkSharing<<endl;
        if(nTrkSharing>1) {
            //cout<<"found one hit shared by "<<nTrkSharing<<" trks: layer "<<(*itHit)->getLayer()
            //<<", pos("<<(*itHit)->getHitPosition().x()<<","<<(*itHit)->getHitPosition().y()<<") with ";
            int trkId_minDelD=-1;
            double minResid = 99.0;
            //vector<int> trkIdList_toDel; 
            vector<double> vecRes = (*itHit)->getVecResid();
            int nTrk = vecRes.size(); 
            //cout<<nTrk<<" residuals: ";
            for(int i=0; i<nTrk; i++)// find out the minimum residual
            {
                //cout<<"   trk "<<trkIdList[i]<<":"<<vecRes[i];
                if(minResid>fabs(vecRes[i])) {
                    minResid=fabs(vecRes[i]);
                    trkId_minDelD=trkIdList[i];
                }
            }
            //cout<<endl;
            //cout<<"    trkId_minDelD, minResid = "<<trkId_minDelD<<", "<<minResid<<endl;
            for(int i=0; i<nTrk; i++)// remove other associations
            {
                if(trkIdList[i]!=trkId_minDelD) {
                    (*itHit)->dropTrkID(trkIdList[i]);// remove the trk from the hit
                    vector<HoughTrack>::iterator itTrk = getHoughTrkIt(m_houghTrackList,trkIdList[i]);
                    if(itTrk!=m_houghTrackList.end()) itTrk->dropHitPnt(&(*(*itHit)));// remove the hit from the trk
                }
            }
            trkIdList = (*itHit)->getTrkID();
            nTrkSharing = trkIdList.size();
            //cout<<"new nTrkSharing = "<<nTrkSharing<<endl;
        }// if(nTrkSharing>1)
    }// loop hits
}
 
 
vector<HoughTrack>::iterator HoughFinder::getHoughTrkIt(vector<HoughTrack>& houghTrkList, int trkId)
{
    vector<HoughTrack>::iterator it = houghTrkList.begin();
    for(; it!=houghTrkList.end(); it++) 
    {
        if(it->getTrkID()==trkId) break;
    }
    return it;
}
 
//void HoughFinder::updateXHitsStatistics(vector<HoughTrack>& trackVector)
//{
//}
 
int HoughFinder::associateVHits()
{
    //m_debug=true;
    //int trkSize = m_houghTrackList.size();
    //cout<<"found "<<trkSize<<" circles on x-y plane"<<endl;
    if(m_debug==5) 
    {
        cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
        cout<<" start V hits association: "<<endl;
    }
    //cout<<__FILE__<<"   "<<__LINE__<<endl;
    for(vector<HoughTrack>::iterator trkIT = m_houghTrackList.begin(); trkIT!=m_houghTrackList.end(); trkIT++)// loop circles
    {
        // --- check circle track
        //trkIT->print();
        //trkIT->printHot();
        //HoughTrack* trkCandi = m_houghTrackList[i];
        if(trkIT->getFlag()!=0) continue; // skip bad circle candidates
        trkIT->updateLayerStat();
        if(m_debug==5) {
            cout<<"try V-hits association for circle: "<<trkIT->getTrkID()<<endl;
            cout<<"if near stereo hits: "<<trkIT->nearStereoHits()<<endl;
        }
        trkIT->clearXVhitsCgem();
        
        // --- fill tanl-dz map
        int nXVhit = 0;
        m_roughTanlDzMap.Reset();
        int vote = 3;
        //nXVhit = fillHistogram(m_VHoughHitList,&m_roughTanlDzMap,0,vote,trkCandi);
        nXVhit = fillHistogram(m_VHoughHitList,&m_roughTanlDzMap,0,vote,&(*trkIT));// fill tanLambda dz map
        if(m_debug==5) {
            cout<<nXVhit<<" nXVhits filled in the sz Hough map"<<endl;
        }
        double x_peak(0.), y_peak(0.);
        double x_weight(0.), y_weight(0.);
        if(nXVhit<3) 
        {
            trkIT->setFlag(3); // not enough V hits
            if(m_debug==5) {
                cout<<"not enough V hits "<<endl;
                cout<<endl;
            }
            continue;// skip circle with Vhits candidates <3
        }
        else if(m_nVClusterOnSZmap>0.5*nXVhit)// need to self V-cluster selectionn
        {
            // FIXME, to be implemented
        }
        //else {
        getWeightedPeak(m_roughTanlDzMap, x_peak, y_peak, x_weight, y_weight);
        //cout<<"rough x_peak, y_peak = "<<x_peak<<", "<<y_peak<<endl;
        //cout<<"rough x_weight, y_weight = "<<x_weight<<", "<<y_weight<<endl;
        //int binx(0),biny(0),binz(0);
        //m_roughTanlDzMap.GetMaximumBin(binx,biny,binz);
        //int peakHeight = m_roughTanlDzMap.GetBinContent(binx,biny);
        //cout<<"peakHeight:"<<peakHeight<<endl;
        //cout<<endl;
        //cout<<endl;
        //x_peak = x_weight;
        //y_peak = y_weight;
        //-----------------------------------------------------------------
        //binning study
        //}
        
        double cosTh=x_weight/sqrt(x_weight*x_weight+1);
        if(m_debug==5) {
            cout<<" cosTh="<<cosTh<<", dz="<<y_weight<<endl;
        }
        trkIT->setTanl(x_weight);
        trkIT->setDz(y_weight);
        trkIT->updateHelix();// update pivot and a
        int flag = 1;
        int charge = trkIT->getCharge();
 
        //int nHot = trkIT->findHot(m_VHoughHitList, flag, charge,3,m_event);
        //trkIT->dropRedundentCgemVHits();
        //trkIT->setMdcHit( &m_mdcHitCol);
        //trkIT->print();
        //int l(20);
        //int l = Layer; l = 35;
        //while(l<36){
        //        //while(1)
        //        //cout<<"----------------------------------------------------------------"<<endl;
        //        //cout<<"----"<<endl;
        //        double dr = trkIT->dr();
        //        double phi0 = trkIT->phi0();
        //        double kappa = trkIT->kappa();
        //        double dz = trkIT->dz();
        //        double tanl = trkIT->tanl();
        //        double cutFactor=1.0;
        //        if(nXVhit<4) cutFactor=3;
        //        //int nHot = trkIT->findVHot(m_VHoughHitList,charge,l,cutFactor);
        //        int nHot = trkIT->findVHot(m_VHoughHitListOnSZmap,charge,l,cutFactor);
        //        trkIT->dropRedundentCgemVHits();
        //        if(m_debug) 
        //        {
        //          cout<<"track id "<<trkIT->getTrkID()<<", maxLayer="<<l<<", nVhits="<<nHot<<endl;
        //          trkIT->print();
        //        }
        //        //cout<<nHot<<endl;
        //        //trkIT->print();
        //        vector<HoughHit*> hot = trkIT->getHotList();
        //        //cout<<"nHot: "<<hot.size()<<endl;
        //        //cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
        //        //trkIT->print();
        //        //trkIT->printHot();
        // 
        //        TrkErrCode trkErrCode;
        //        if(m_fitFlag>1) {
        //          trkErrCode = trkIT->fitHelix(m_mdcDetector,m_bfield,m_bunchT0,hot,l);
        //          //trkIT->fitHelix(&myDotsHelixFitter,m_bunchT0*1e9,m_recCgemClusterColBegin);
        //        }
        //        //trkIT->print();
        //        //cout<<"----------------------------------------------------------------"<<endl;
        //        bool isFitDivergency(false);
        //        isFitDivergency = isFitDivergency||fabs(trkIT->dr()    -dr   )>999.;
        //        isFitDivergency = isFitDivergency||fabs(trkIT->phi0()  -phi0 )>999.;
        //        isFitDivergency = isFitDivergency||fabs(trkIT->kappa() -kappa)>999.;
        //        isFitDivergency = isFitDivergency||fabs(trkIT->dz()    -dz   )>999.;
        //        isFitDivergency = isFitDivergency||fabs(trkIT->tanl()  -tanl )>999.;
        //        if(isFitDivergency){
        //          trkIT->setFlag(-2);
        //          if(m_debug)cout<<"helix fit diverged!"<<endl;
        //          break;
        //        }
        //        if(trkErrCode.failure()){
        //          //if(m_debug)cout<<__FILE__<<"  line"<<__LINE__<<":helix fit failure!"<<endl;
        //          trkIT->setFlag(-1);
        //          if(m_debug){
        //              trkErrCode.print(std::cout);cout<<endl;//cout<<endl;
        //              cout<<"helix fit failed!"<<endl;
        //          }
        //        }
        //        else{
        //          trkIT->setFlag(0);
        //          if(m_mcTruth)findMcTrack(&(*trkIT));
        //          if(m_debug)cout<<"helix fit good!"<<endl;
        //        }
        //        if(m_debug) 
        //        {
        //          cout<<"track info after fit: "<<endl;
        //          trkIT->print();
        //        }
        // 
        //        l++;
        //        if(l>m_maxFireLayer)break;
        //        //break;
        //        //if(nHot==0)break;
        //}// while
        // 
        //if(trkIT->getFlag()!=0)
        //{
        //        //cout<<"trk flag !=0 !"<<endl;"
        //        continue;
        //}
        ////cout<<"nHot: "<<nHot<<endl;
        ////trkIT->print();
        ////trkIT->printHot();
        //double dr = trkIT->dr();
        //double phi0 = trkIT->phi0();
        //double kappa = trkIT->kappa();
        //double dz = trkIT->dz();
        //double tanl = trkIT->tanl();
        // 
        //double cutFactor=1.0;
        //if(nXVhit<4) cutFactor=3;
        ////int nHot = trkIT->findVHot(m_VHoughHitList,charge,l,cutFactor);
        //int nHot = trkIT->findVHot(m_VHoughHitListOnSZmap,charge,l,cutFactor);
        //trkIT->dropRedundentCgemVHits();
        //if(m_debug) cout<<"after helix fit, nVhits="<<nHot<<endl;
        //vector<HoughHit*> hot = trkIT->getHotList();
        //TrkErrCode trkErrCode = trkIT->fitHelix(m_mdcDetector,m_trkContextEv,m_bunchT0,m_mdcHitCol,hot);
        // 
        //bool isFitDivergency(false);
        //isFitDivergency = isFitDivergency||fabs(trkIT->dr()    -dr   )>999.;
        //isFitDivergency = isFitDivergency||fabs(trkIT->phi0()  -phi0 )>999.;
        //isFitDivergency = isFitDivergency||fabs(trkIT->kappa() -kappa)>999.;
        //isFitDivergency = isFitDivergency||fabs(trkIT->dz()    -dz   )>999.;
        //isFitDivergency = isFitDivergency||fabs(trkIT->tanl()  -tanl )>999.;
        //if(isFitDivergency){
        //        trkIT->setFlag(-2);
        //        if(m_debug)cout<<"helix fit diverged!"<<endl;
        //        break;
        //}
        //if(trkErrCode.failure()){
        //        //if(m_debug)cout<<__FILE__<<"  line"<<__LINE__<<":helix fit failure!"<<endl;
        //        trkIT->setFlag(-1);
        //        if(m_debug)cout<<"helix fit failed!"<<endl;
        //}
        //else{
        //        trkIT->setFlag(0);
        //        if(m_mcTruth)findMcTrack(&(*trkIT));
        //        if(m_debug)cout<<"helix fit good!"<<endl;
        //}
        //if(trkIT->getFlag()!=0)
        //{
        //        if(m_debug) cout<<"helix flag !=0 !"<<endl;
        //}
 
        // --- V-hits selection and helix fit
        int helixFitFlag(0);
        int nHot_last=0;
        int lmax=35;
        //while(1)
        for(int j=0; j<2; j++)// fit Helix twice
        {
            double cutFactor=1.0;
            if(nXVhit<4) cutFactor=3;
            if(fabs(cosTh)>0.8) cutFactor=5;
            int nHot = trkIT->findVHot(m_VHoughHitListOnSZmap,charge,lmax,cutFactor); // includes a sz line fit, to be done FIXME
            trkIT->dropRedundentCgemVHits();// at most one V-cluster on one layer
            // gap check, and remove hits after a big distance // to be done FIXME; could be solved by checking if any axial hit around
            if(m_debug==5) trkIT->printHot(3);
            helixFitFlag = trkIT->fitHelix(&myDotsHelixFitter,m_bunchT0*1e9,m_recCgemClusterColBegin, m_chi2CutHits);// trkIT's m_dr etc. updated already
            if(helixFitFlag==0)
            {
                trkIT->updateHelix();// update pivot and a
                //if(nHot_last!=nHot) 
                //{
                //nHot_last=nHot;
                //continue;
                //}
                //else break;
            }
            else break;
        }
        if(helixFitFlag!=0) // if fit fails
        {
            trkIT->setFlag(helixFitFlag);
            if(m_debug==5) {
                cout<<"Helix fit failed "<<endl;
                cout<<endl;
            }
            continue; // skip bad helix fit
        }
        else // if fit successes
        {
            vector<RecMdcHit> aRecMdcHitVec=myDotsHelixFitter.makeRecMdcHitVec(1);
            trkIT->setRecMdcHitVec(aRecMdcHitVec);
            //trkIT->setRecMdcHitVec(myDotsHelixFitter.makeRecMdcHitVec(1));
            if(m_debug==5) {
                cout<<"Helix fit success and aRecMdcHitVec filled "<<endl;
                trkIT->printRecMdcHitVec();
                cout<<endl;
                //trkIT->printHot();
            }
        }
 
        //cout<<"getId() "<<m_trkContextEv->getId()<<endl;
        //trkIT->print();
        //cout<<trkIT->dr()    -dr   <<endl;
        //cout<<trkIT->phi0()  -phi0 <<endl;
        //cout<<trkIT->kappa() -kappa<<endl;
        //cout<<trkIT->dz()    -dz   <<endl;
        //cout<<trkIT->tanl()  -tanl <<endl;
        //trkIT->printHot();
    }// loop circles
    //cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
    //m_debug=false;
    return 0;
}
 
void HoughFinder::findMcTrack(HoughTrack* track)
{
    TH1I trkID("track index","",100,0,100);
    vector<HoughHit*> hotList = track->getHotList();
    for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end(); hotIt++){
        if((*hotIt)->getPairHit()!=NULL)trkID.Fill(((*hotIt)->getPairHit()->getTrkID())[0]);
    }
    int trkid = trkID.GetMaximumBin()-1;
    for(vector<HoughTrack>::iterator trkIter = m_mcTrackCol.begin();trkIter!=m_mcTrackCol.end();trkIter++){
        if(trkIter->getTrkID()==trkid){
            if(m_run<0&&m_mcTruth)track->setMcTrack(&(*trkIter));
            break;
        }
    }
}
 
int HoughFinder::nFineBinTanl(double tanl)
{
}
 
int HoughFinder::nFineBinDz(double tanl)
{
}
 
void HoughFinder::XVhitCutWindow(double tanl, int ilayer, double charge, double& cut1, double &cut2)
{
}
 
//if(m_hist)ntuple_hit->write();
StatusCode HoughFinder::bookTuple(){
    MsgStream log(msgSvc(), name());
 
    NTuplePtr nt1(ntupleSvc(), "mdcHoughFinder/hit");
    if ( nt1 ){
        ntuple_hit= nt1;
    } else {
        ntuple_hit= ntupleSvc()->book("mdcHoughFinder/hit", CLID_ColumnWiseTuple, "hit");
        if(ntuple_hit){
            ntuple_hit->addItem("hit_run",          m_hit_run);
            ntuple_hit->addItem("hit_event",        m_hit_event);
 
            ntuple_hit->addItem("hit_nhit",         m_hit_nhit, 0, 10000);
            ntuple_hit->addItem("hit_hitID",        m_hit_nhit, m_hit_hitID);
            ntuple_hit->addItem("hit_hitType",      m_hit_nhit, m_hit_hitType);
            ntuple_hit->addItem("hit_layer",        m_hit_nhit, m_hit_layer);
            ntuple_hit->addItem("hit_wire",         m_hit_nhit, m_hit_wire);
            ntuple_hit->addItem("hit_flag",         m_hit_nhit, m_hit_flag);
            ntuple_hit->addItem("hit_halfCircle",   m_hit_nhit, m_hit_halfCircle);
            ntuple_hit->addItem("hit_x",            m_hit_nhit, m_hit_x);
            ntuple_hit->addItem("hit_y",            m_hit_nhit, m_hit_y);
            ntuple_hit->addItem("hit_z",            m_hit_nhit, m_hit_z);
            ntuple_hit->addItem("hit_drift",        m_hit_nhit, m_hit_drift);
 
            ntuple_hit->addItem("mcHit_hitID",      m_hit_nhit, m_mcHit_hitID);
            ntuple_hit->addItem("mcHit_hitType",    m_hit_nhit, m_mcHit_hitType);
            ntuple_hit->addItem("mcHit_layer",      m_hit_nhit, m_mcHit_layer);
            ntuple_hit->addItem("mcHit_wire",       m_hit_nhit, m_mcHit_wire);
            ntuple_hit->addItem("mcHit_flag",       m_hit_nhit, m_mcHit_flag);
            ntuple_hit->addItem("mcHit_halfCircle", m_hit_nhit, m_mcHit_halfCircle);
            ntuple_hit->addItem("mcHit_x",          m_hit_nhit, m_mcHit_x);
            ntuple_hit->addItem("mcHit_y",          m_hit_nhit, m_mcHit_y);
            ntuple_hit->addItem("mcHit_z",          m_hit_nhit, m_mcHit_z);
            ntuple_hit->addItem("mcHit_drift",      m_hit_nhit, m_mcHit_drift);
        } else { log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/hit" <<endmsg;
            return StatusCode::FAILURE;
        }
    }
 
    NTuplePtr nt2(ntupleSvc(), "mdcHoughFinder/track");
    if ( nt2 ){
        ntuple_track = nt2;
    } else {
        ntuple_track = ntupleSvc()->book("mdcHoughFinder/track", CLID_ColumnWiseTuple, "track");
        if(ntuple_track){
            ntuple_track->addItem("trk_run",            m_trk_run);
            ntuple_track->addItem("trk_event",          m_trk_event);
            ntuple_track->addItem("trk_nTrack",         m_trk_nTrack);
            ntuple_track->addItem("trk_trackID",        m_trk_trackID);
            ntuple_track->addItem("trk_charge",         m_trk_charge);
            ntuple_track->addItem("trk_flag",           m_trk_flag);
            ntuple_track->addItem("trk_angle",          m_trk_angle);
            ntuple_track->addItem("trk_rho",            m_trk_rho);
            ntuple_track->addItem("trk_dAngle",         m_trk_dAngle);
            ntuple_track->addItem("trk_dRho",           m_trk_dRho);
            ntuple_track->addItem("trk_dTanl",          m_trk_dTanl);
            ntuple_track->addItem("trk_dDz",            m_trk_dDz);
            ntuple_track->addItem("trk_Xc",             m_trk_Xc);
            ntuple_track->addItem("trk_Yc",             m_trk_Yc);
            ntuple_track->addItem("trk_R",              m_trk_R);
            ntuple_track->addItem("trk_dr",             m_trk_dr);
            ntuple_track->addItem("trk_phi0",           m_trk_phi0);
            ntuple_track->addItem("trk_kappa",          m_trk_kappa);
            ntuple_track->addItem("trk_dz",             m_trk_dz);
            ntuple_track->addItem("trk_tanl",           m_trk_tanl);
            ntuple_track->addItem("trk_pxy",            m_trk_pxy);
            ntuple_track->addItem("trk_px",             m_trk_px);
            ntuple_track->addItem("trk_py",             m_trk_py);
            ntuple_track->addItem("trk_pz",             m_trk_pz);
            ntuple_track->addItem("trk_p",              m_trk_p);
            ntuple_track->addItem("trk_phi",            m_trk_phi);
            ntuple_track->addItem("trk_theta",          m_trk_theta);
            ntuple_track->addItem("trk_cosTheta",       m_trk_cosTheta);
            ntuple_track->addItem("trk_vx",             m_trk_vx);
            ntuple_track->addItem("trk_vy",             m_trk_vy);
            ntuple_track->addItem("trk_vz",             m_trk_vz);
            ntuple_track->addItem("trk_vr",             m_trk_vr);
            ntuple_track->addItem("trk_chi2",           m_trk_chi2);
            ntuple_track->addItem("trk_fiTerm",         m_trk_fiTerm);
            ntuple_track->addItem("trk_nhit",           m_trk_nhit);
            ntuple_track->addItem("trk_ncluster",       m_trk_ncluster);
            ntuple_track->addItem("trk_stat",           m_trk_stat);
            ntuple_track->addItem("trk_ndof",           m_trk_ndof);
            ntuple_track->addItem("trk_nster",          m_trk_nster);
            ntuple_track->addItem("trk_nlayer",         m_trk_nlayer);
            ntuple_track->addItem("trk_firstLayer",     m_trk_firstLayer);
            ntuple_track->addItem("trk_lastLayer",      m_trk_lastLayer);
            ntuple_track->addItem("trk_nCgemXClusters", m_trk_nCgemXClusters);
            ntuple_track->addItem("trk_nCgemVClusters", m_trk_nCgemVClusters);
 
            ntuple_track->addItem("trk_nHot",           m_trk_nHot, 0, 10000);
            ntuple_track->addItem("hot_hitID",          m_trk_nHot, m_hot_hitID);
            ntuple_track->addItem("hot_hitType",        m_trk_nHot, m_hot_hitType);
            ntuple_track->addItem("hot_layer",          m_trk_nHot, m_hot_layer);
            ntuple_track->addItem("hot_wire",           m_trk_nHot, m_hot_wire);
            ntuple_track->addItem("hot_flag",           m_trk_nHot, m_hot_flag);
            ntuple_track->addItem("hot_halfCircle",     m_trk_nHot, m_hot_halfCircle);
            ntuple_track->addItem("hot_x",              m_trk_nHot, m_hot_x);
            ntuple_track->addItem("hot_y",              m_trk_nHot, m_hot_y);
            ntuple_track->addItem("hot_z",              m_trk_nHot, m_hot_z);
            ntuple_track->addItem("hot_drift",          m_trk_nHot, m_hot_drift);
            ntuple_track->addItem("hot_residual",       m_trk_nHot, m_hot_residual);
 
            ntuple_track->addItem("mcHot_hitID",        m_trk_nHot, m_mcHot_hitID);
            ntuple_track->addItem("mcHot_hitType",      m_trk_nHot, m_mcHot_hitType);
            ntuple_track->addItem("mcHot_layer",        m_trk_nHot, m_mcHot_layer);
            ntuple_track->addItem("mcHot_wire",         m_trk_nHot, m_mcHot_wire);
            ntuple_track->addItem("mcHot_flag",         m_trk_nHot, m_mcHot_flag);
            ntuple_track->addItem("mcHot_halfCircle",   m_trk_nHot, m_mcHot_halfCircle);
            ntuple_track->addItem("mcHot_x",            m_trk_nHot, m_mcHot_x);
            ntuple_track->addItem("mcHot_y",            m_trk_nHot, m_mcHot_y);
            ntuple_track->addItem("mcHot_z",            m_trk_nHot, m_mcHot_z);
            ntuple_track->addItem("mcHot_drift",        m_trk_nHot, m_mcHot_drift);
 
            ntuple_track->addItem("mcTrk_trackID",        m_mcTrk_trackID);
            ntuple_track->addItem("mcTrk_charge",         m_mcTrk_charge);
            ntuple_track->addItem("mcTrk_flag",           m_mcTrk_flag);
            ntuple_track->addItem("mcTrk_angle",          m_mcTrk_angle);
            ntuple_track->addItem("mcTrk_rho",            m_mcTrk_rho);
            ntuple_track->addItem("mcTrk_dAngle",         m_mcTrk_dAngle);
            ntuple_track->addItem("mcTrk_dRho",           m_mcTrk_dRho);
            ntuple_track->addItem("mcTrk_dTanl",          m_mcTrk_dTanl);
            ntuple_track->addItem("mcTrk_dDz",            m_mcTrk_dDz);
            ntuple_track->addItem("mcTrk_Xc",             m_mcTrk_Xc);
            ntuple_track->addItem("mcTrk_Yc",             m_mcTrk_Yc);
            ntuple_track->addItem("mcTrk_R",              m_mcTrk_R);
            ntuple_track->addItem("mcTrk_dr",             m_mcTrk_dr);
            ntuple_track->addItem("mcTrk_phi0",           m_mcTrk_phi0);
            ntuple_track->addItem("mcTrk_kappa",          m_mcTrk_kappa);
            ntuple_track->addItem("mcTrk_dz",             m_mcTrk_dz);
            ntuple_track->addItem("mcTrk_tanl",           m_mcTrk_tanl);
            ntuple_track->addItem("mcTrk_pxy",            m_mcTrk_pxy);
            ntuple_track->addItem("mcTrk_px",             m_mcTrk_px);
            ntuple_track->addItem("mcTrk_py",             m_mcTrk_py);
            ntuple_track->addItem("mcTrk_pz",             m_mcTrk_pz);
            ntuple_track->addItem("mcTrk_p",              m_mcTrk_p);
            ntuple_track->addItem("mcTrk_phi",            m_mcTrk_phi);
            ntuple_track->addItem("mcTrk_theta",          m_mcTrk_theta);
            ntuple_track->addItem("mcTrk_cosTheta",       m_mcTrk_cosTheta);
            ntuple_track->addItem("mcTrk_vx",             m_mcTrk_vx);
            ntuple_track->addItem("mcTrk_vy",             m_mcTrk_vy);
            ntuple_track->addItem("mcTrk_vz",             m_mcTrk_vz);
            ntuple_track->addItem("mcTrk_vr",             m_mcTrk_vr);
            ntuple_track->addItem("mcTrk_chi2",           m_mcTrk_chi2);
            ntuple_track->addItem("mcTrk_fiTerm",         m_mcTrk_fiTerm);
            ntuple_track->addItem("mcTrk_nhit",           m_mcTrk_nhit);
            ntuple_track->addItem("mcTrk_ncluster",       m_mcTrk_ncluster);
            ntuple_track->addItem("mcTrk_stat",           m_mcTrk_stat);
            ntuple_track->addItem("mcTrk_ndof",           m_mcTrk_ndof);
            ntuple_track->addItem("mcTrk_nster",          m_mcTrk_nster);
            ntuple_track->addItem("mcTrk_nlayer",         m_mcTrk_nlayer);
            ntuple_track->addItem("mcTrk_firstLayer",     m_mcTrk_firstLayer);
            ntuple_track->addItem("mcTrk_lastLayer",      m_mcTrk_lastLayer);
            ntuple_track->addItem("mcTrk_nCgemXClusters", m_mcTrk_nCgemXClusters);
            ntuple_track->addItem("mcTrk_nCgemVClusters", m_mcTrk_nCgemVClusters);
 
            ntuple_track->addItem("mcTrk_nHot",           m_mcTrk_nHot, 0, 10000);
            ntuple_track->addItem("mcTrkHot_hitID",       m_mcTrk_nHot, m_mcTrkHot_hitID);
            ntuple_track->addItem("mcTrkHot_hitType",     m_mcTrk_nHot, m_mcTrkHot_hitType);
            ntuple_track->addItem("mcTrkHot_layer",       m_mcTrk_nHot, m_mcTrkHot_layer);
            ntuple_track->addItem("mcTrkHot_wire",        m_mcTrk_nHot, m_mcTrkHot_wire);
            ntuple_track->addItem("mcTrkHot_flag",        m_mcTrk_nHot, m_mcTrkHot_flag);
            ntuple_track->addItem("mcTrkHot_halfCircle",  m_mcTrk_nHot, m_mcTrkHot_halfCircle);
            ntuple_track->addItem("mcTrkHot_x",           m_mcTrk_nHot, m_mcTrkHot_x);
            ntuple_track->addItem("mcTrkHot_y",           m_mcTrk_nHot, m_mcTrkHot_y);
            ntuple_track->addItem("mcTrkHot_z",           m_mcTrk_nHot, m_mcTrkHot_z);
            ntuple_track->addItem("mcTrkHot_drift",       m_mcTrk_nHot, m_mcTrkHot_drift);
        } else {
            log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/track" <<endmsg;
            return StatusCode::FAILURE;
        }
    }
 
    NTuplePtr nt3(ntupleSvc(), "mdcHoughFinder/event");
    if ( nt3 ){
        ntuple_event = nt3;
    } else {
        ntuple_event = ntupleSvc()->book("mdcHoughFinder/event", CLID_ColumnWiseTuple, "event");
        if(ntuple_event){
            ntuple_event->addItem("evt_run",                  m_evt_run);
            ntuple_event->addItem("evt_event",                m_evt_event);
            ntuple_event->addItem("evt_nXCluster",            m_evt_nXCluster);
            ntuple_event->addItem("evt_nVCluster",            m_evt_nVCluster);
            ntuple_event->addItem("evt_nXVCluster",           m_evt_nXVCluster);
            ntuple_event->addItem("evt_nCgemCluster",         m_evt_nCGEMCluster);
            ntuple_event->addItem("evt_nAxialHit",            m_evt_nAxialHit);
            ntuple_event->addItem("evt_nStereoHit",           m_evt_nStereoHit);
            ntuple_event->addItem("evt_nODCHit",              m_evt_nODCHit);
            ntuple_event->addItem("evt_nHit",                 m_evt_nHit);
 
            ntuple_event->addItem("evt_nTrack",                  m_evt_nTrack, 0, 100);
            ntuple_event->addItem("evtTrk_trackID",              m_evt_nTrack, m_evtTrk_trackID);
            ntuple_event->addItem("evtTrk_charge",               m_evt_nTrack, m_evtTrk_charge);
            ntuple_event->addItem("evtTrk_flag",                 m_evt_nTrack, m_evtTrk_flag);
            ntuple_event->addItem("evtTrk_angle",                m_evt_nTrack, m_evtTrk_angle);
            ntuple_event->addItem("evtTrk_rho",                  m_evt_nTrack, m_evtTrk_rho);
            ntuple_event->addItem("evtTrk_dAngle",               m_evt_nTrack, m_evtTrk_dAngle);
            ntuple_event->addItem("evtTrk_dRho",                 m_evt_nTrack, m_evtTrk_dRho);
            ntuple_event->addItem("evtTrk_dTanl",                m_evt_nTrack, m_evtTrk_dTanl);
            ntuple_event->addItem("evtTrk_dDz",                  m_evt_nTrack, m_evtTrk_dDz);
            ntuple_event->addItem("evtTrk_Xc",                   m_evt_nTrack, m_evtTrk_Xc);
            ntuple_event->addItem("evtTrk_Yc",                   m_evt_nTrack, m_evtTrk_Yc);
            ntuple_event->addItem("evtTrk_R",                    m_evt_nTrack, m_evtTrk_R);
            ntuple_event->addItem("evtTrk_dr",                   m_evt_nTrack, m_evtTrk_dr);
            ntuple_event->addItem("evtTrk_phi0",                 m_evt_nTrack, m_evtTrk_phi0);
            ntuple_event->addItem("evtTrk_kappa",                m_evt_nTrack, m_evtTrk_kappa);
            ntuple_event->addItem("evtTrk_dz",                   m_evt_nTrack, m_evtTrk_dz);
            ntuple_event->addItem("evtTrk_tanl",                 m_evt_nTrack, m_evtTrk_tanl);
            ntuple_event->addItem("evtTrk_pxy",                  m_evt_nTrack, m_evtTrk_pxy);
            ntuple_event->addItem("evtTrk_px",                   m_evt_nTrack, m_evtTrk_px);
            ntuple_event->addItem("evtTrk_py",                   m_evt_nTrack, m_evtTrk_py);
            ntuple_event->addItem("evtTrk_pz",                   m_evt_nTrack, m_evtTrk_pz);
            ntuple_event->addItem("evtTrk_p",                    m_evt_nTrack, m_evtTrk_p);
            ntuple_event->addItem("evtTrk_phi",                  m_evt_nTrack, m_evtTrk_phi);
            ntuple_event->addItem("evtTrk_theta",                m_evt_nTrack, m_evtTrk_theta);
            ntuple_event->addItem("evtTrk_cosTheta",             m_evt_nTrack, m_evtTrk_cosTheta);
            ntuple_event->addItem("evtTrk_vx",                   m_evt_nTrack, m_evtTrk_vx);
            ntuple_event->addItem("evtTrk_vy",                   m_evt_nTrack, m_evtTrk_vy);
            ntuple_event->addItem("evtTrk_vz",                   m_evt_nTrack, m_evtTrk_vz);
            ntuple_event->addItem("evtTrk_vr",                   m_evt_nTrack, m_evtTrk_vr);
            ntuple_event->addItem("evtTrk_chi2",                 m_evt_nTrack, m_evtTrk_chi2);
            ntuple_event->addItem("evtTrk_fiTerm",               m_evt_nTrack, m_evtTrk_fiTerm);
            ntuple_event->addItem("evtTrk_nhit",                 m_evt_nTrack, m_evtTrk_nhit);
            ntuple_event->addItem("evtTrk_ncluster",             m_evt_nTrack, m_evtTrk_ncluster);
            ntuple_event->addItem("evtTrk_stat",                 m_evt_nTrack, m_evtTrk_stat);
            ntuple_event->addItem("evtTrk_ndof",                 m_evt_nTrack, m_evtTrk_ndof);
            ntuple_event->addItem("evtTrk_nster",                m_evt_nTrack, m_evtTrk_nster);
            ntuple_event->addItem("evtTrk_nlayer",               m_evt_nTrack, m_evtTrk_nlayer);
            ntuple_event->addItem("evtTrk_firstLayer",           m_evt_nTrack, m_evtTrk_firstLayer);
            ntuple_event->addItem("evtTrk_lastLayer",            m_evt_nTrack, m_evtTrk_lastLayer);
            ntuple_event->addItem("evtTrk_nCgemXClusters",       m_evt_nTrack, m_evtTrk_nCgemXClusters);
            ntuple_event->addItem("evtTrk_nCgemVClusters",       m_evt_nTrack, m_evtTrk_nCgemVClusters);
 
            ntuple_event->addItem("mcEvtTrk_trackID",              m_evt_nTrack, m_mcEvtTrk_trackID);
            ntuple_event->addItem("mcEvtTrk_charge",               m_evt_nTrack, m_mcEvtTrk_charge);
            ntuple_event->addItem("mcEvtTrk_flag",                 m_evt_nTrack, m_mcEvtTrk_flag);
            ntuple_event->addItem("mcEvtTrk_angle",                m_evt_nTrack, m_mcEvtTrk_angle);
            ntuple_event->addItem("mcEvtTrk_rho",                  m_evt_nTrack, m_mcEvtTrk_rho);
            ntuple_event->addItem("mcEvtTrk_dAngle",               m_evt_nTrack, m_mcEvtTrk_dAngle);
            ntuple_event->addItem("mcEvtTrk_dRho",                 m_evt_nTrack, m_mcEvtTrk_dRho);
            ntuple_event->addItem("mcEvtTrk_dTanl",                m_evt_nTrack, m_mcEvtTrk_dTanl);
            ntuple_event->addItem("mcEvtTrk_dDz",                  m_evt_nTrack, m_mcEvtTrk_dDz);
            ntuple_event->addItem("mcEvtTrk_Xc",                   m_evt_nTrack, m_mcEvtTrk_Xc);
            ntuple_event->addItem("mcEvtTrk_Yc",                   m_evt_nTrack, m_mcEvtTrk_Yc);
            ntuple_event->addItem("mcEvtTrk_R",                    m_evt_nTrack, m_mcEvtTrk_R);
            ntuple_event->addItem("mcEvtTrk_dr",                   m_evt_nTrack, m_mcEvtTrk_dr);
            ntuple_event->addItem("mcEvtTrk_phi0",                 m_evt_nTrack, m_mcEvtTrk_phi0);
            ntuple_event->addItem("mcEvtTrk_kappa",                m_evt_nTrack, m_mcEvtTrk_kappa);
            ntuple_event->addItem("mcEvtTrk_dz",                   m_evt_nTrack, m_mcEvtTrk_dz);
            ntuple_event->addItem("mcEvtTrk_tanl",                 m_evt_nTrack, m_mcEvtTrk_tanl);
            ntuple_event->addItem("mcEvtTrk_pxy",                  m_evt_nTrack, m_mcEvtTrk_pxy);
            ntuple_event->addItem("mcEvtTrk_px",                   m_evt_nTrack, m_mcEvtTrk_px);
            ntuple_event->addItem("mcEvtTrk_py",                   m_evt_nTrack, m_mcEvtTrk_py);
            ntuple_event->addItem("mcEvtTrk_pz",                   m_evt_nTrack, m_mcEvtTrk_pz);
            ntuple_event->addItem("mcEvtTrk_p",                    m_evt_nTrack, m_mcEvtTrk_p);
            ntuple_event->addItem("mcEvtTrk_phi",                  m_evt_nTrack, m_mcEvtTrk_phi);
            ntuple_event->addItem("mcEvtTrk_theta",                m_evt_nTrack, m_mcEvtTrk_theta);
            ntuple_event->addItem("mcEvtTrk_cosTheta",             m_evt_nTrack, m_mcEvtTrk_cosTheta);
            ntuple_event->addItem("mcEvtTrk_vx",                   m_evt_nTrack, m_mcEvtTrk_vx);
            ntuple_event->addItem("mcEvtTrk_vy",                   m_evt_nTrack, m_mcEvtTrk_vy);
            ntuple_event->addItem("mcEvtTrk_vz",                   m_evt_nTrack, m_mcEvtTrk_vz);
            ntuple_event->addItem("mcEvtTrk_vr",                   m_evt_nTrack, m_mcEvtTrk_vr);
            ntuple_event->addItem("mcEvtTrk_chi2",                 m_evt_nTrack, m_mcEvtTrk_chi2);
            ntuple_event->addItem("mcEvtTrk_fiTerm",               m_evt_nTrack, m_mcEvtTrk_fiTerm);
            ntuple_event->addItem("mcEvtTrk_nhit",                 m_evt_nTrack, m_mcEvtTrk_nhit);
            ntuple_event->addItem("mcEvtTrk_ncluster",             m_evt_nTrack, m_mcEvtTrk_ncluster);
            ntuple_event->addItem("mcEvtTrk_stat",                 m_evt_nTrack, m_mcEvtTrk_stat);
            ntuple_event->addItem("mcEvtTrk_ndof",                 m_evt_nTrack, m_mcEvtTrk_ndof);
            ntuple_event->addItem("mcEvtTrk_nster",                m_evt_nTrack, m_mcEvtTrk_nster);
            ntuple_event->addItem("mcEvtTrk_nlayer",               m_evt_nTrack, m_mcEvtTrk_nlayer);
            ntuple_event->addItem("mcEvtTrk_firstLayer",           m_evt_nTrack, m_mcEvtTrk_firstLayer);
            ntuple_event->addItem("mcEvtTrk_lastLayer",            m_evt_nTrack, m_mcEvtTrk_lastLayer);
            ntuple_event->addItem("mcEvtTrk_nCgemXClusters",       m_evt_nTrack, m_mcEvtTrk_nCgemXClusters);
            ntuple_event->addItem("mcEvtTrk_nCgemVClusters",       m_evt_nTrack, m_mcEvtTrk_nCgemVClusters);
        } else {
            log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/event" <<endmsg;
            return StatusCode::FAILURE;
        }
    }
    return StatusCode::SUCCESS;
}
 
int HoughFinder::dumpHit()
{
    m_hit_run   = m_run;
    m_hit_event = m_event;
    m_hit_nhit  = m_houghHitList.size();
    int i(0);
    for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end(); hitIt++){
        m_hit_hitID[i]      = hitIt->getHitID();
        m_hit_hitType[i]    = hitIt->getHitType();
        m_hit_layer[i]      = hitIt->getLayer();
        m_hit_wire[i]       = hitIt->getWire();
        m_hit_flag[i]       = hitIt->getFlag();
        m_hit_halfCircle[i] = hitIt->getHalfCircle();
        m_hit_x[i]          = hitIt->getHitPosition().x();
        m_hit_y[i]          = hitIt->getHitPosition().y();
        m_hit_z[i]          = hitIt->getHitPosition().z();
        m_hit_drift[i]      = hitIt->getDriftDist();
 
        if(hitIt->getPairHit()!=NULL){
            m_mcHit_hitID[i]      = hitIt->getPairHit()->getHitID();
            m_mcHit_hitType[i]    = hitIt->getPairHit()->getHitType();
            m_mcHit_layer[i]      = hitIt->getPairHit()->getLayer();
            m_mcHit_wire[i]       = hitIt->getPairHit()->getWire();
            m_mcHit_flag[i]       = hitIt->getPairHit()->getFlag();
            m_mcHit_halfCircle[i] = hitIt->getPairHit()->getHalfCircle();
            m_mcHit_x[i]          = hitIt->getPairHit()->getHitPosition().x();
            m_mcHit_y[i]          = hitIt->getPairHit()->getHitPosition().y();
            m_mcHit_z[i]          = hitIt->getPairHit()->getHitPosition().z();
            m_mcHit_drift[i]      = hitIt->getPairHit()->getDriftDist();
        }
        i++;
    }
    ntuple_hit->write();
    return m_houghHitList.size();
}
 
int HoughFinder::dumpHoughTrack()
{
    if(m_houghTrackList.size()==0)return 0;
    m_trk_run               = m_run;
    m_trk_event             = m_event;
    m_trk_nTrack            = m_houghTrackList.size();
    HepPoint3D positionOntrack(999,999,999), positionOnDetector(999,999,999);
    for(vector<HoughTrack>::iterator trkIt = m_houghTrackList.begin(); trkIt != m_houghTrackList.end(); trkIt++){
        m_trk_trackID           = trkIt->getTrkID();
        m_trk_charge            = trkIt->getCharge();
        m_trk_flag              = trkIt->getFlag();
        m_trk_angle             = trkIt->getAngle();
        m_trk_rho               = trkIt->getRho();
        m_trk_dAngle            = trkIt->getDAngle();
        m_trk_dRho              = trkIt->getDRho();
        m_trk_dTanl             = trkIt->getDTanl();
        m_trk_dDz               = trkIt->getDDz();
        m_trk_Xc                = trkIt->center().x();
        m_trk_Yc                = trkIt->center().y();
        m_trk_R                 = trkIt->radius();
        m_trk_dr                = trkIt->dr();
        m_trk_phi0              = trkIt->phi0();
        m_trk_kappa             = trkIt->kappa();
        m_trk_dz                = trkIt->dz();
        m_trk_tanl              = trkIt->tanl();
        m_trk_pxy               = trkIt->pt();
        m_trk_px                = trkIt->momentum(0).x();
        m_trk_py                = trkIt->momentum(0).y();
        m_trk_pz                = trkIt->momentum(0).z();
        m_trk_p                 = trkIt->momentum(0).mag();
        m_trk_phi               = trkIt->direction(0).phi();
        m_trk_theta             = trkIt->direction(0).theta();
        m_trk_cosTheta          = cos(trkIt->direction(0).theta());
        m_trk_vx                = trkIt->x(0).x();
        m_trk_vy                = trkIt->x(0).y();
        m_trk_vz                = trkIt->x(0).z();
        m_trk_vr                = trkIt->x(0).perp();
        m_trk_chi2              = trkIt->getChi2();
        //m_trk_chi2[i]           = trkIt->
        //m_trk_fiTerm[i]         = trkIt->
        //m_trk_nhit[i]           = trkIt->
        //m_trk_ncluster[i]       = trkIt->
        //m_trk_stat[i]           = trkIt->
        //m_trk_ndof[i]           = trkIt->
        //m_trk_nster[i]          = trkIt->
        //m_trk_nlayer[i]         = trkIt->
        //m_trk_firstLayer[i]     = trkIt->
        //m_trk_lastLayer[i]      = trkIt->
        //m_trk_nCgemXClusters[i] = trkIt->
        //m_trk_nCgemVClusters[i] = trkIt->
 
        vector<HoughHit*> hotList = trkIt->getHotList();
        m_trk_nHot              = hotList.size();
        int i(0);
        for(vector<HoughHit*>::iterator hotIt = hotList.begin(); hotIt != hotList.end(); hotIt++){
            m_hot_hitID[i]      = (*hotIt)->getHitID();
            m_hot_hitType[i]    = (*hotIt)->getHitType();
            m_hot_layer[i]      = (*hotIt)->getLayer();
            m_hot_wire[i]       = (*hotIt)->getWire();
            m_hot_flag[i]       = (*hotIt)->getFlag();
            m_hot_halfCircle[i] = (*hotIt)->getHalfCircle();
            m_hot_x[i]          = (*hotIt)->getHitPosition().x();
            m_hot_y[i]          = (*hotIt)->getHitPosition().y();
            m_hot_z[i]          = (*hotIt)->getHitPosition().z();
            m_hot_drift[i]      = (*hotIt)->getDriftDist();
            m_hot_residual[i]   = (*hotIt)->residual(&(*trkIt),positionOntrack, positionOnDetector);
 
            if((*hotIt)->getPairHit()!=NULL){
                m_mcHot_hitID[i]      = (*hotIt)->getPairHit()->getHitID();
                m_mcHot_hitType[i]    = (*hotIt)->getPairHit()->getHitType();
                m_mcHot_layer[i]      = (*hotIt)->getPairHit()->getLayer();
                m_mcHot_wire[i]       = (*hotIt)->getPairHit()->getWire();
                m_mcHot_flag[i]       = (*hotIt)->getPairHit()->getFlag();
                m_mcHot_halfCircle[i] = (*hotIt)->getPairHit()->getHalfCircle();
                m_mcHot_x[i]          = (*hotIt)->getPairHit()->getHitPosition().x();
                m_mcHot_y[i]          = (*hotIt)->getPairHit()->getHitPosition().y();
                m_mcHot_z[i]          = (*hotIt)->getPairHit()->getHitPosition().z();
                m_mcHot_drift[i]      = (*hotIt)->getPairHit()->getDriftDist();
            }
            i++;
        }
 
        if(trkIt->getMcTrack()!=NULL){
            m_mcTrk_trackID        = trkIt->getMcTrack()->getTrkID();
            m_mcTrk_charge         = trkIt->getMcTrack()->getCharge();
            m_mcTrk_flag           = trkIt->getMcTrack()->getFlag();
            m_mcTrk_angle          = trkIt->getMcTrack()->getAngle();
            m_mcTrk_rho            = trkIt->getMcTrack()->getRho();
            m_mcTrk_dAngle         = trkIt->getMcTrack()->getDAngle();
            m_mcTrk_dRho           = trkIt->getMcTrack()->getDRho();
            m_mcTrk_dTanl          = trkIt->getMcTrack()->getDTanl();
            m_mcTrk_dDz            = trkIt->getMcTrack()->getDDz();
            m_mcTrk_Xc             = trkIt->getMcTrack()->center().x();
            m_mcTrk_Yc             = trkIt->getMcTrack()->center().y();
            m_mcTrk_R              = trkIt->getMcTrack()->radius();
            m_mcTrk_dr             = trkIt->getMcTrack()->dr();
            m_mcTrk_phi0           = trkIt->getMcTrack()->phi0();
            m_mcTrk_kappa          = trkIt->getMcTrack()->kappa();
            m_mcTrk_dz             = trkIt->getMcTrack()->dz();
            m_mcTrk_tanl           = trkIt->getMcTrack()->tanl();
            m_mcTrk_pxy            = trkIt->getMcTrack()->pt();
            m_mcTrk_px             = trkIt->getMcTrack()->momentum(0).x();
            m_mcTrk_py             = trkIt->getMcTrack()->momentum(0).y();
            m_mcTrk_pz             = trkIt->getMcTrack()->momentum(0).z();
            m_mcTrk_p              = trkIt->getMcTrack()->momentum(0).mag();
            m_mcTrk_phi            = trkIt->getMcTrack()->direction(0).phi();
            m_mcTrk_theta          = trkIt->getMcTrack()->direction(0).theta();
            m_mcTrk_cosTheta       = cos(trkIt->getMcTrack()->direction(0).theta());
            m_mcTrk_vx             = trkIt->getMcTrack()->x(0).x();
            m_mcTrk_vy             = trkIt->getMcTrack()->x(0).y();
            m_mcTrk_vz             = trkIt->getMcTrack()->x(0).z();
            m_mcTrk_vr             = trkIt->getMcTrack()->x(0).perp();
            //m_mcTrk_chi2           = trkIt->getMcTrack()->
            //m_mcTrk_fiTerm         = trkIt->getMcTrack()->
            //m_mcTrk_nhit           = trkIt->getMcTrack()->
            //m_mcTrk_ncluster       = trkIt->getMcTrack()->
            //m_mcTrk_stat           = trkIt->getMcTrack()->
            //m_mcTrk_ndof           = trkIt->getMcTrack()->
            //m_mcTrk_nster          = trkIt->getMcTrack()->
            //m_mcTrk_nlayer         = trkIt->getMcTrack()->
            //m_mcTrk_firstLayer     = trkIt->getMcTrack()->
            //m_mcTrk_lastLayer      = trkIt->getMcTrack()->
            //m_mcTrk_nCgemXClusters = trkIt->getMcTrack()->
            //m_mcTrk_nCgemVClusters = trkIt->getMcTrack()->
 
            vector<HoughHit*> mcHotList = trkIt->getMcTrack()->getHotList();
            m_mcTrk_nHot            = mcHotList.size();
            int j(0);
            for(vector<HoughHit*>::iterator mcHotIt = mcHotList.begin(); mcHotIt != mcHotList.end(); mcHotIt++){
                m_mcTrkHot_hitID[j]      = (*mcHotIt)->getHitID();
                m_mcTrkHot_hitType[j]    = (*mcHotIt)->getHitType();
                m_mcTrkHot_layer[j]      = (*mcHotIt)->getLayer();
                m_mcTrkHot_wire[j]       = (*mcHotIt)->getWire();
                m_mcTrkHot_flag[j]       = (*mcHotIt)->getFlag();
                m_mcTrkHot_halfCircle[j] = (*mcHotIt)->getHalfCircle();
                m_mcTrkHot_x[j]          = (*mcHotIt)->getHitPosition().x();
                m_mcTrkHot_y[j]          = (*mcHotIt)->getHitPosition().y();
                m_mcTrkHot_z[j]          = (*mcHotIt)->getHitPosition().z();
                m_mcTrkHot_drift[j]      = (*mcHotIt)->getDriftDist();
                j++;
            }
        }
        ntuple_track->write();
    }
    return m_houghTrackList.size();
}
 
int HoughFinder::dumpHoughEvent()
{
    if(m_houghTrackList.size()==0)return 0;
    m_evt_run                   = m_run;
    m_evt_event                 = m_event;
    m_evt_nTrack                = m_houghTrackList.size();
    int i(0);
    for(vector<HoughTrack>::iterator trkIt = m_houghTrackList.begin(); trkIt != m_houghTrackList.end(); trkIt++){
 
        m_evtTrk_trackID[i]        = trkIt->getTrkID();
        m_evtTrk_charge[i]         = trkIt->getCharge();
        m_evtTrk_flag[i]           = trkIt->getFlag();
        m_evtTrk_angle[i]          = trkIt->getAngle();
        m_evtTrk_rho[i]            = trkIt->getRho();
        m_evtTrk_dAngle[i]         = trkIt->getDAngle();
        m_evtTrk_dRho[i]           = trkIt->getDRho();
        m_evtTrk_dTanl[i]          = trkIt->getDTanl();
        m_evtTrk_dDz[i]            = trkIt->getDDz();
        m_evtTrk_Xc[i]             = trkIt->center().x();
        m_evtTrk_Yc[i]             = trkIt->center().y();
        m_evtTrk_R[i]              = trkIt->radius();
        m_evtTrk_dr[i]             = trkIt->dr();
        m_evtTrk_phi0[i]           = trkIt->phi0();
        m_evtTrk_kappa[i]          = trkIt->kappa();
        m_evtTrk_dz[i]             = trkIt->dz();
        m_evtTrk_tanl[i]           = trkIt->tanl();
        m_evtTrk_pxy[i]            = trkIt->pt();
        m_evtTrk_px[i]             = trkIt->momentum(0).x();
        m_evtTrk_py[i]             = trkIt->momentum(0).y();
        m_evtTrk_pz[i]             = trkIt->momentum(0).z();
        m_evtTrk_p[i]              = trkIt->momentum(0).mag();
        m_evtTrk_phi[i]            = trkIt->direction(0).phi();
        m_evtTrk_theta[i]          = trkIt->direction(0).theta();
        m_evtTrk_cosTheta[i]       = cos(trkIt->direction(0).theta());
        m_evtTrk_vx[i]             = trkIt->x(0).x();
        m_evtTrk_vy[i]             = trkIt->x(0).y();
        m_evtTrk_vz[i]             = trkIt->x(0).z();
        m_evtTrk_vr[i]             = trkIt->x(0).perp();
 
        //m_evtTrk_chi2[i]           = trkIt->
        //m_evtTrk_fiTerm[i]         = trkIt->
        //m_evtTrk_nhit[i]           = trkIt->
        //m_evtTrk_ncluster[i]       = trkIt->
        //m_evtTrk_stat[i]           = trkIt->
        //m_evtTrk_ndof[i]           = trkIt->
        //m_evtTrk_nster[i]          = trkIt->
        //m_evtTrk_nlayer[i]         = trkIt->
        //m_evtTrk_firstLayer[i]     = trkIt->
        //m_evtTrk_lastLayer[i]      = trkIt->
        //m_evtTrk_nCgemXClusters[i] = trkIt->
        //m_evtTrk_nCgemVClusters[i] = trkIt->
        if(trkIt->getMcTrack()!=NULL){
            m_mcEvtTrk_trackID[i]        = trkIt->getMcTrack()->getTrkID();
            m_mcEvtTrk_charge[i]         = trkIt->getMcTrack()->getCharge();
            m_mcEvtTrk_flag[i]           = trkIt->getMcTrack()->getFlag();
            m_mcEvtTrk_angle[i]          = trkIt->getMcTrack()->getAngle();
            m_mcEvtTrk_rho[i]            = trkIt->getMcTrack()->getRho();
            m_mcEvtTrk_dAngle[i]         = trkIt->getMcTrack()->getDAngle();
            m_mcEvtTrk_dRho[i]           = trkIt->getMcTrack()->getDRho();
            m_mcEvtTrk_dTanl[i]          = trkIt->getMcTrack()->getDTanl();
            m_mcEvtTrk_dDz[i]            = trkIt->getMcTrack()->getDDz();
            m_mcEvtTrk_Xc[i]             = trkIt->getMcTrack()->center().x();
            m_mcEvtTrk_Yc[i]             = trkIt->getMcTrack()->center().y();
            m_mcEvtTrk_R[i]              = trkIt->getMcTrack()->radius();
            m_mcEvtTrk_dr[i]             = trkIt->getMcTrack()->dr();
            m_mcEvtTrk_phi0[i]           = trkIt->getMcTrack()->phi0();
            m_mcEvtTrk_kappa[i]          = trkIt->getMcTrack()->kappa();
            m_mcEvtTrk_dz[i]             = trkIt->getMcTrack()->dz();
            m_mcEvtTrk_tanl[i]           = trkIt->getMcTrack()->tanl();
            m_mcEvtTrk_pxy[i]            = trkIt->getMcTrack()->pt();
            m_mcEvtTrk_px[i]             = trkIt->getMcTrack()->momentum(0).x();
            m_mcEvtTrk_py[i]             = trkIt->getMcTrack()->momentum(0).y();
            m_mcEvtTrk_pz[i]             = trkIt->getMcTrack()->momentum(0).z();
            m_mcEvtTrk_p[i]              = trkIt->getMcTrack()->momentum(0).mag();
            m_mcEvtTrk_phi[i]            = trkIt->getMcTrack()->direction(0).phi();
            m_mcEvtTrk_theta[i]          = trkIt->getMcTrack()->direction(0).theta();
            m_mcEvtTrk_cosTheta[i]       = cos(trkIt->getMcTrack()->direction(0).theta());
            m_mcEvtTrk_vx[i]             = trkIt->getMcTrack()->x(0).x();
            m_mcEvtTrk_vy[i]             = trkIt->getMcTrack()->x(0).y();
            m_mcEvtTrk_vz[i]             = trkIt->getMcTrack()->x(0).z();
            m_mcEvtTrk_vr[i]             = trkIt->getMcTrack()->x(0).perp();
            //m_mcEvtTrk_chi2[i]           = trkIt->getMcTrack()->
            //m_mcEvtTrk_fiTerm[i]         = trkIt->getMcTrack()->
            //m_mcEvtTrk_nhit[i]           = trkIt->getMcTrack()->
            //m_mcEvtTrk_ncluster[i]       = trkIt->getMcTrack()->
            //m_mcEvtTrk_stat[i]           = trkIt->getMcTrack()->
            //m_mcEvtTrk_ndof[i]           = trkIt->getMcTrack()->
            //m_mcEvtTrk_nster[i]          = trkIt->getMcTrack()->
            //m_mcEvtTrk_nlayer[i]         = trkIt->getMcTrack()->
            //m_mcEvtTrk_firstLayer[i]     = trkIt->getMcTrack()->
            //m_mcEvtTrk_lastLayer[i]      = trkIt->getMcTrack()->
            //m_mcEvtTrk_nCgemXClusters[i] = trkIt->getMcTrack()->
            //m_mcEvtTrk_nCgemVClusters[i] = trkIt->getMcTrack()->
        }
        i++;
    }
    ntuple_event->write();
    return m_houghTrackList.size();
}
 
int HoughFinder::dumpTdsTrack()
{
    SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
    if(!recMdcTrackCol)return 0;
    m_evt_run                   = m_run;
    m_evt_event                 = m_event;
    m_evt_nTrack                = recMdcTrackCol->size();
    for(RecMdcTrackCol::iterator iter = recMdcTrackCol->begin();iter!=recMdcTrackCol->end();iter++){
        m_trk_trackID             = (*iter)->trackId()      ;
        m_trk_charge              = (*iter)->charge()       ;
        //m_trk_flag                = (*iter)->              
        //m_trk_angle               = (*iter)->              
        //m_trk_rho                 = (*iter)->              
        //m_trk_dAngle              = (*iter)->              
        //m_trk_dRho                = (*iter)->              
        //m_trk_dTanl               = (*iter)->              
        //m_trk_dDz                 = (*iter)->              
        //m_trk_Xc                  = (*iter)->              
        //m_trk_Yc                  = (*iter)->              
        //m_trk_R                   = (*iter)->              
        m_trk_dr                  = (*iter)->helix(0)       ;
        m_trk_phi0                = (*iter)->helix(1)       ;
        m_trk_kappa               = (*iter)->helix(2)       ;
        m_trk_dz                  = (*iter)->helix(3)       ;
        m_trk_tanl                = (*iter)->helix(4)       ;
        m_trk_pxy                 = (*iter)->pxy()          ;
        m_trk_px                  = (*iter)->px()           ;
        m_trk_py                  = (*iter)->py()           ;
        m_trk_pz                  = (*iter)->pz()           ;
        m_trk_p                   = (*iter)->p()            ;
        m_trk_phi                 = (*iter)->theta()        ;
        m_trk_theta               = (*iter)->phi()          ;
        m_trk_cosTheta            = cos((*iter)->theta())   ;
        m_trk_vx                  = (*iter)->x()            ;
        m_trk_vy                  = (*iter)->y()            ;
        m_trk_vz                  = (*iter)->z()            ;
        m_trk_vr                  = (*iter)->r()            ;
        m_trk_chi2                = (*iter)->chi2()         ;
        m_trk_fiTerm              = (*iter)->getFiTerm()    ;
        m_trk_nhit                = (*iter)->getNhits()     ;
        m_trk_ncluster            = (*iter)->getNcluster()  ;
        m_trk_stat                = (*iter)->stat()         ;
        m_trk_ndof                = (*iter)->ndof()         ;
        m_trk_nster               = (*iter)->nster()        ;
        m_trk_nlayer              = (*iter)->nlayer()       ;
        m_trk_firstLayer          = (*iter)->firstLayer()   ;
        m_trk_lastLayer           = (*iter)->lastLayer()    ;
        m_trk_nCgemXClusters      = (*iter)->nCgemXCluster();
        m_trk_nCgemVClusters      = (*iter)->nCgemVCluster();
 
        HitRefVec hitRefVec         = (*iter)->getVecHits();
        HitRefVec::iterator hotIt   = hitRefVec.begin();
        ClusterRefVec clusterRefVec = (*iter)->getVecClusters();
        ClusterRefVec::iterator clusterIt = clusterRefVec.begin();
        m_trk_nHot = clusterRefVec.size() + hitRefVec.size();
        int i(0);
        for(clusterIt=clusterRefVec.begin();clusterIt!=clusterRefVec.end();clusterIt++){
            for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end(); hitIt++){
                if(hitIt->getHitType()!=HoughHit::CGEMHIT)continue;
                if((*clusterIt)->getclusterid()==(hitIt)->getCgemCluster()->getclusterid()){
                    m_hot_hitID[i]      = (hitIt)->getHitID();
                    m_hot_hitType[i]    = (hitIt)->getHitType();
                    m_hot_layer[i]      = (hitIt)->getLayer();
                    m_hot_wire[i]       = (hitIt)->getWire();
                    m_hot_flag[i]       = (hitIt)->getFlag();
                    m_hot_halfCircle[i] = (hitIt)->getHalfCircle();
                    m_hot_x[i]          = (hitIt)->getHitPosition().x();
                    m_hot_y[i]          = (hitIt)->getHitPosition().y();
                    m_hot_z[i]          = (hitIt)->getHitPosition().z();
                    m_hot_drift[i]      = (hitIt)->getDriftDist();
                    m_hot_residual[i]   = -999;
 
                    if((hitIt)->getPairHit()!=NULL){
                        m_mcHot_hitID[i]      = (hitIt)->getPairHit()->getHitID();
                        m_mcHot_hitType[i]    = (hitIt)->getPairHit()->getHitType();
                        m_mcHot_layer[i]      = (hitIt)->getPairHit()->getLayer();
                        m_mcHot_wire[i]       = (hitIt)->getPairHit()->getWire();
                        m_mcHot_flag[i]       = (hitIt)->getPairHit()->getFlag();
                        m_mcHot_halfCircle[i] = (hitIt)->getPairHit()->getHalfCircle();
                        m_mcHot_x[i]          = (hitIt)->getPairHit()->getHitPosition().x();
                        m_mcHot_y[i]          = (hitIt)->getPairHit()->getHitPosition().y();
                        m_mcHot_z[i]          = (hitIt)->getPairHit()->getHitPosition().z();
                        m_mcHot_drift[i]      = (hitIt)->getPairHit()->getDriftDist();
                    }
                    i++;
                }
            }
        }
        for(hotIt=hitRefVec.begin();hotIt!=hitRefVec.end();hotIt++){
            int layer = MdcID::layer((*hotIt)->getMdcId());
            int wire = MdcID::wire((*hotIt)->getMdcId());
            for(HitVector_Iterator hitIt = m_houghHitList.begin(); hitIt != m_houghHitList.end(); hitIt++){
                if(hitIt->getHitType()!=HoughHit::MDCHIT)continue;
                if(layer == (hitIt)->getLayer() && wire == (hitIt)->getWire()){
                    m_hot_hitID[i]      = (hitIt)->getHitID();
                    m_hot_hitType[i]    = (hitIt)->getHitType();
                    m_hot_layer[i]      = (hitIt)->getLayer();
                    m_hot_wire[i]       = (hitIt)->getWire();
                    m_hot_flag[i]       = (hitIt)->getFlag();
                    m_hot_halfCircle[i] = (hitIt)->getHalfCircle();
                    m_hot_x[i]          = (hitIt)->getHitPosition().x();
                    m_hot_y[i]          = (hitIt)->getHitPosition().y();
                    m_hot_z[i]          = (hitIt)->getHitPosition().z();
                    m_hot_drift[i]      = (hitIt)->getDriftDist();
                    m_hot_residual[i]   = 999;
 
                    if((hitIt)->getPairHit()!=NULL){
                        m_mcHot_hitID[i]      = (hitIt)->getPairHit()->getHitID();
                        m_mcHot_hitType[i]    = (hitIt)->getPairHit()->getHitType();
                        m_mcHot_layer[i]      = (hitIt)->getPairHit()->getLayer();
                        m_mcHot_wire[i]       = (hitIt)->getPairHit()->getWire();
                        m_mcHot_flag[i]       = (hitIt)->getPairHit()->getFlag();
                        m_mcHot_halfCircle[i] = (hitIt)->getPairHit()->getHalfCircle();
                        m_mcHot_x[i]          = (hitIt)->getPairHit()->getHitPosition().x();
                        m_mcHot_y[i]          = (hitIt)->getPairHit()->getHitPosition().y();
                        m_mcHot_z[i]          = (hitIt)->getPairHit()->getHitPosition().z();
                        m_mcHot_drift[i]      = (hitIt)->getPairHit()->getDriftDist();
                    }
                }
            }
            i++;
        }
 
        int maxSameHitNumber(0);
        vector<HoughTrack>::iterator sameTrkIt = m_mcTrackCol.end();
        vector<HoughTrack>::iterator trkIter = m_mcTrackCol.begin();
        for(;trkIter!=m_mcTrackCol.end();trkIter++){
            int sameHitNumber(0);
            vector<HoughHit*> hitList = trkIter->getHotList();
            for(vector<HoughHit*>::iterator mchotIt = hitList.begin(); mchotIt != hitList.end();mchotIt++){
                if((*mchotIt)->getHitType()==HoughHit::CGEMHIT){
                    for(clusterIt=clusterRefVec.begin();clusterIt!=clusterRefVec.end();clusterIt++){
                        if((*clusterIt)->getclusterid()==(*mchotIt)->getCgemCluster()->getclusterid()){
                            sameHitNumber++;
                        }
                    }
                }else if((*mchotIt)->getHitType()==HoughHit::MDCHIT){
                    for(hotIt=hitRefVec.begin();hotIt!=hitRefVec.end();hotIt++){
                        int layer = MdcID::layer((*hotIt)->getMdcId());
                        int wire = MdcID::wire((*hotIt)->getMdcId());
                        if(layer == (*mchotIt)->getLayer() && wire == (*mchotIt)->getWire()){
                            sameHitNumber++;
                        }
                    }
                }
            }
            if(sameHitNumber>maxSameHitNumber){
                maxSameHitNumber = sameHitNumber;
                sameTrkIt = trkIter;
            }
 
        }
 
        if(sameTrkIt!=m_mcTrackCol.end()){
            m_mcTrk_trackID           = sameTrkIt->getTrkID();
            m_mcTrk_charge            = sameTrkIt->getCharge();
            m_mcTrk_flag              = sameTrkIt->getFlag();
            m_mcTrk_angle             = sameTrkIt->getAngle();
            m_mcTrk_rho               = sameTrkIt->getRho();
            m_mcTrk_dAngle            = sameTrkIt->getDAngle();
            m_mcTrk_dRho              = sameTrkIt->getDRho();
            m_mcTrk_dTanl             = sameTrkIt->getDTanl();
            m_mcTrk_dDz               = sameTrkIt->getDDz();
            m_mcTrk_Xc                = sameTrkIt->center().x();
            m_mcTrk_Yc                = sameTrkIt->center().y();
            m_mcTrk_R                 = sameTrkIt->radius();
            m_mcTrk_dr                = sameTrkIt->dr();
            m_mcTrk_phi0              = sameTrkIt->phi0();
            m_mcTrk_kappa             = sameTrkIt->kappa();
            m_mcTrk_dz                = sameTrkIt->dz();
            m_mcTrk_tanl              = sameTrkIt->tanl();
            m_mcTrk_pxy               = sameTrkIt->pt();
            m_mcTrk_px                = sameTrkIt->momentum(0).x();
            m_mcTrk_py                = sameTrkIt->momentum(0).y();
            m_mcTrk_pz                = sameTrkIt->momentum(0).z();
            m_mcTrk_p                 = sameTrkIt->momentum(0).mag();
            m_mcTrk_phi               = sameTrkIt->direction(0).phi();
            m_mcTrk_theta             = sameTrkIt->direction(0).theta();
            m_mcTrk_cosTheta          = cos(sameTrkIt->direction(0).theta());
            m_mcTrk_vx                = sameTrkIt->x(0).x();
            m_mcTrk_vy                = sameTrkIt->x(0).y();
            m_mcTrk_vz                = sameTrkIt->x(0).z();
            m_mcTrk_vr                = sameTrkIt->x(0).perp();
            //m_mcTrk_chi2              = sameTrkIt->getChi2();
            //m_mcTrk_fiTerm            = sameTrkIt->
            //m_mcTrk_nhit              = (sameTrkIt->getHotList()).size();
            //m_mcTrk_ncluster          = sameTrkIt->
            //m_mcTrk_stat              = sameTrkIt->
            //m_mcTrk_ndof              = sameTrkIt->
            //m_mcTrk_nster             = sameTrkIt->
            //m_mcTrk_nlayer            = sameTrkIt->
            //m_mcTrk_firstLayer        = sameTrkIt->
            //m_mcTrk_lastLayer         = sameTrkIt->
            //m_mcTrk_nC emXClusters    = sameTrkIt-> 
            //m_mcTrk_nCgemVClusters    = sameTrkIt-> 
            vector<HoughHit*> mcHotList = sameTrkIt->getHotList();
            int j(0);
            for(vector<HoughHit*>::iterator mcHotIt = mcHotList.begin(); mcHotIt != mcHotList.end(); mcHotIt++){
                m_mcTrkHot_hitID[j]      = (*mcHotIt)->getHitID();
                m_mcTrkHot_hitType[j]    = (*mcHotIt)->getHitType();
                m_mcTrkHot_layer[j]      = (*mcHotIt)->getLayer();
                m_mcTrkHot_wire[j]       = (*mcHotIt)->getWire();
                m_mcTrkHot_flag[j]       = (*mcHotIt)->getFlag();
                m_mcTrkHot_halfCircle[j] = (*mcHotIt)->getHalfCircle();
                m_mcTrkHot_x[j]          = (*mcHotIt)->getHitPosition().x();
                m_mcTrkHot_y[j]          = (*mcHotIt)->getHitPosition().y();
                m_mcTrkHot_z[j]          = (*mcHotIt)->getHitPosition().z();
                m_mcTrkHot_drift[j]      = (*mcHotIt)->getDriftDist();
                j++;
            }
        }
        ntuple_track->write();
    }
    return recMdcTrackCol->size();
}
 
int HoughFinder::dumpTdsEvent() 
{
    SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
    if(!recMdcTrackCol)return 0;
    m_evt_run                   = m_run;
    m_evt_event                 = m_event;
    m_evt_nTrack                = recMdcTrackCol->size();
    int i(0);
    for(RecMdcTrackCol::iterator iter = recMdcTrackCol->begin();iter!=recMdcTrackCol->end();iter++){
        m_evtTrk_trackID[i]            = (*iter)->trackId()      ;
        m_evtTrk_charge[i]             = (*iter)->charge()       ;
        //m_evtTrk_flag[i]               = (*iter)->              
        //m_evtTrk_angle[i]              = (*iter)->              
        //m_evtTrk_rho[i]                = (*iter)->              
        //m_evtTrk_dAngle[i]             = (*iter)->              
        //m_evtTrk_dRho[i]               = (*iter)->              
        //m_evtTrk_dTanl[i]              = (*iter)->              
        //m_evtTrk_dDz[i]                = (*iter)->              
        //m_evtTrk_Xc[i]                 = (*iter)->              
        //m_evtTrk_Yc[i]                 = (*iter)->              
        //m_evtTrk_R[i]                  = (*iter)->              
        m_evtTrk_dr[i]                 = (*iter)->helix(0)       ;
        m_evtTrk_phi0[i]               = (*iter)->helix(1)       ;
        m_evtTrk_kappa[i]              = (*iter)->helix(2)       ;
        m_evtTrk_dz[i]                 = (*iter)->helix(3)       ;
        m_evtTrk_tanl[i]               = (*iter)->helix(4)       ;
        m_evtTrk_pxy[i]                = (*iter)->pxy()          ;
        m_evtTrk_px[i]                 = (*iter)->px()           ;
        m_evtTrk_py[i]                 = (*iter)->py()           ;
        m_evtTrk_pz[i]                 = (*iter)->pz()           ;
        m_evtTrk_p[i]                  = (*iter)->p()            ;
        m_evtTrk_phi[i]                = (*iter)->theta()        ;
        m_evtTrk_theta[i]              = (*iter)->phi()          ;
        m_evtTrk_cosTheta[i]           = cos((*iter)->theta())   ;
        m_evtTrk_vx[i]                 = (*iter)->x()            ;
        m_evtTrk_vy[i]                 = (*iter)->y()            ;
        m_evtTrk_vz[i]                 = (*iter)->z()            ;
        m_evtTrk_vr[i]                 = (*iter)->r()            ;
        m_evtTrk_chi2[i]               = (*iter)->chi2()         ;
        m_evtTrk_fiTerm[i]             = (*iter)->getFiTerm()    ;
        m_evtTrk_nhit[i]               = (*iter)->getNhits()     ;
        m_evtTrk_ncluster[i]           = (*iter)->getNcluster()  ;
        m_evtTrk_stat[i]               = (*iter)->stat()         ;
        m_evtTrk_ndof[i]               = (*iter)->ndof()         ;
        m_evtTrk_nster[i]              = (*iter)->nster()        ;
        m_evtTrk_nlayer[i]             = (*iter)->nlayer()       ;
        m_evtTrk_firstLayer[i]         = (*iter)->firstLayer()   ;
        m_evtTrk_lastLayer[i]          = (*iter)->lastLayer()    ;
        m_evtTrk_nCgemXClusters[i]     = (*iter)->nCgemXCluster();
        m_evtTrk_nCgemVClusters[i]     = (*iter)->nCgemVCluster();
 
        HitRefVec hitRefVec         = (*iter)->getVecHits();
        HitRefVec::iterator hitIt   = hitRefVec.begin();
        ClusterRefVec clusterRefVec = (*iter)->getVecClusters();
        ClusterRefVec::iterator clusterIt = clusterRefVec.begin();
 
        int maxSameHitNumber(0);
        vector<HoughTrack>::iterator sameTrkIt = m_mcTrackCol.end();
        vector<HoughTrack>::iterator trkIter = m_mcTrackCol.begin();
        for(;trkIter!=m_mcTrackCol.end();trkIter++){
            int sameHitNumber(0);
            vector<HoughHit*> hitList = trkIter->getHotList();
            for(vector<HoughHit*>::iterator hotIt = hitList.begin(); hotIt != hitList.end();hotIt++){
                if((*hotIt)->getHitType()==HoughHit::CGEMHIT){
                    for(;clusterIt!=clusterRefVec.end();clusterIt++){
                        if((*clusterIt)->getclusterid()==(*hotIt)->getCgemCluster()->getclusterid()){
                            sameHitNumber++;
                        }
                    }
                }else if((*hotIt)->getHitType()==HoughHit::MDCHIT){
                    for(;hitIt!=hitRefVec.end();hitIt++){
                        int layer = MdcID::layer((*hitIt)->getMdcId());
                        int wire = MdcID::wire((*hitIt)->getMdcId());
                        if(layer == (*hotIt)->getLayer() && wire == (*hotIt)->getWire()){
                            sameHitNumber++;
                        }
                    }
                }
            }
            if(sameHitNumber>maxSameHitNumber){
                maxSameHitNumber = sameHitNumber;
                sameTrkIt = trkIter;
            }
 
        }
 
        if(sameTrkIt!=m_mcTrackCol.end()){
            m_mcEvtTrk_trackID[i]          = sameTrkIt->getTrkID();
            m_mcEvtTrk_charge[i]           = sameTrkIt->getCharge();
            m_mcEvtTrk_flag[i]             = sameTrkIt->getFlag();
            m_mcEvtTrk_angle[i]            = sameTrkIt->getAngle();
            m_mcEvtTrk_rho[i]              = sameTrkIt->getRho();
            m_mcEvtTrk_dAngle[i]           = sameTrkIt->getDAngle();
            m_mcEvtTrk_dRho[i]             = sameTrkIt->getDRho();
            m_mcEvtTrk_dTanl[i]            = sameTrkIt->getDTanl();
            m_mcEvtTrk_dDz[i]              = sameTrkIt->getDDz();
            m_mcEvtTrk_Xc[i]               = sameTrkIt->center().x();
            m_mcEvtTrk_Yc[i]               = sameTrkIt->center().y();
            m_mcEvtTrk_R[i]                = sameTrkIt->radius();
            m_mcEvtTrk_dr[i]               = sameTrkIt->dr();
            m_mcEvtTrk_phi0[i]             = sameTrkIt->phi0();
            m_mcEvtTrk_kappa[i]            = sameTrkIt->kappa();
            m_mcEvtTrk_dz[i]               = sameTrkIt->dz();
            m_mcEvtTrk_tanl[i]             = sameTrkIt->tanl();
            m_mcEvtTrk_pxy[i]              = sameTrkIt->pt();
            m_mcEvtTrk_px[i]               = sameTrkIt->momentum(0).x();
            m_mcEvtTrk_py[i]               = sameTrkIt->momentum(0).y();
            m_mcEvtTrk_pz[i]               = sameTrkIt->momentum(0).z();
            m_mcEvtTrk_p[i]                = sameTrkIt->momentum(0).mag();
            m_mcEvtTrk_phi[i]              = sameTrkIt->direction(0).phi();
            m_mcEvtTrk_theta[i]            = sameTrkIt->direction(0).theta();
            m_mcEvtTrk_cosTheta[i]         = cos(sameTrkIt->direction(0).theta());
            m_mcEvtTrk_vx[i]                = sameTrkIt->x(0).x();
            m_mcEvtTrk_vy[i]                = sameTrkIt->x(0).y();
            m_mcEvtTrk_vz[i]                = sameTrkIt->x(0).z();
            m_mcEvtTrk_vr[i]                = sameTrkIt->x(0).perp();
            //m_mcEvtTrk_chi2[i]             = sameTrkIt->getChi2();
            //m_mcEvtTrk_fiTerm[i]           = sameTrkIt->
            //m_mcEvtTrk_nhit[i]             = (sameTrkIt->getHotList()).size();
            //m_mcEvtTrk_ncluster[i]         = sameTrkIt->
            //m_mcEvtTrk_stat[i]             = sameTrkIt->
            //m_mcEvtTrk_ndof[i]             = sameTrkIt->
            //m_mcEvtTrk_nster[i]            = sameTrkIt->
            //m_mcEvtTrk_nlayer[i]           = sameTrkIt->
            //m_mcEvtTrk_firstLayer[i]       = sameTrkIt->
            //m_mcEvtTrk_lastLayer[i]        = sameTrkIt->
            //m_mcEvtTrk_nCgemXClusters[i]   = sameTrkIt-> 
            //m_mcEvtTrk_nCgemVClusters[i]   = sameTrkIt-> 
        }
        i++;
    }
    ntuple_event->write();
    return recMdcTrackCol->size();
}
 
int HoughFinder::printTdsTrack()
{
    SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
    if(!recMdcTrackCol)return 0;
    for(RecMdcTrackCol::iterator iter = recMdcTrackCol->begin();iter!=recMdcTrackCol->end();iter++){
        double dr                  = (*iter)->helix(0)       ;
        double phi0                = (*iter)->helix(1)       ;
        double kappa               = (*iter)->helix(2)       ;
        double dz                  = (*iter)->helix(3)       ;
        double tanl                = (*iter)->helix(4)       ;
        cout<<setw(12)<<"dr:"    <<setw(15)<<dr
            <<setw(12)<<"phi0:"  <<setw(15)<<phi0
            <<setw(12)<<"kappa:" <<setw(15)<<kappa
            <<setw(12)<<"dz:"    <<setw(15)<<dz
            <<setw(12)<<"tanl:"  <<setw(15)<<tanl
            <<endl;
 
        ClusterRefVec clusterRefVec = (*iter)->getVecClusters();
        ClusterRefVec::iterator clusterIt = clusterRefVec.begin();
        for(clusterIt=clusterRefVec.begin();clusterIt!=clusterRefVec.end();clusterIt++){
            //cout<<(*clusterIt)->;
        }
 
        HitRefVec hitRefVec         = (*iter)->getVecHits();
        HitRefVec::iterator hotIt   = hitRefVec.begin();
        //cout<<"("<<"layer"<<","<<"wire"<<") "
        cout<<"("<<"l "<<","<<" w "<<") "
            //<<"Stat  "
            //<<"FlagLR  "
            <<"S  "
            <<"F  "
            <<"DriftLeft  "
            <<" DriftRight  "
            <<"ErrDrift_L  "
            <<"ErrDrift_R  "
            <<"ChisqAdd  "
            //<<"  Tdc   "
            //<<"  Adc       "
            <<"  DriftT  "
            <<"  Doca      "
            <<"  Entra     "
            <<"  Zhit    "
            <<"  FltLen  "
            <<endl;
        for(hotIt=hitRefVec.begin();hotIt!=hitRefVec.end();hotIt++){
            int layer = MdcID::layer((*hotIt)->getMdcId());
            int wire = MdcID::wire((*hotIt)->getMdcId());
            cout<<"("<<setw( 2)<<layer<<","<<setw( 3)<<wire<<") "
                <<setw( 3)<<(*hotIt)->getStat()
                <<setw( 3)<<(*hotIt)->getFlagLR()
                <<setw(12)<<(*hotIt)->getDriftDistLeft()
                <<setw(12)<<(*hotIt)->getDriftDistRight()
                <<setw(12)<<(*hotIt)->getErrDriftDistLeft()
                <<setw(12)<<(*hotIt)->getErrDriftDistRight()
                <<setw(12)<<(*hotIt)->getChisqAdd()
                //<<setw( 6)<<(*hotIt)->getTdc()
                //<<setw( 6)<<(*hotIt)->getAdc()
                <<setw(10)<<(*hotIt)->getDriftT()
                <<setw(12)<<(*hotIt)->getDoca()
                <<setw(12)<<(*hotIt)->getEntra()
                <<setw(10)<<(*hotIt)->getZhit()
                <<setw(8 )<<(*hotIt)->getFltLen() 
                <<endl;
        }
        cout<<endl;
    }
    return recMdcTrackCol->size();
}
 
//cout<<__FILE__<<"   "<<__LINE__<<endl;