MdcHoughFinder/MdcHoughFinder-00-00-12/src/Hough.cxx

Go to the documentation of this file.

#include "MdcHoughFinder/Hough.h" 
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/IPartPropSvc.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/IService.h"
#include "GaudiKernel/NTuple.h"                                        
#include "GaudiKernel/INTupleSvc.h"
#include "EventModel/EventHeader.h"
#include "MdcRawEvent/MdcDigi.h"
#include "Identifier/Identifier.h"
#include "Identifier/MdcID.h"
 
#include <iostream>
#include <math.h>
 
#include "EvTimeEvent/RecEsTime.h"
#include "MdcGeom/EntranceAngle.h"
#include "RawEvent/RawDataUtil.h"
#include "MdcData/MdcHit.h"
 
#include "McTruth/DecayMode.h"
#include "McTruth/McParticle.h"
#include "TrackUtil/Helix.h"
#include "MdcRecoUtil/Pdt.h"
 
#include "TrkBase/TrkFit.h"
#include "TrkBase/TrkHitList.h"
#include "TrkBase/TrkExchangePar.h"
#include "TrkFitter/TrkHelixMaker.h"
#include "TrkFitter/TrkCircleMaker.h"
#include "TrkFitter/TrkLineMaker.h"
#include "TrkFitter/TrkHelixFitter.h"
#include "MdcxReco/Mdcxprobab.h"
#include "MdcData/MdcRecoHitOnTrack.h"
#include "MdcPrintSvc/IMdcPrintSvc.h"
 
#include "MdcHoughFinder/HoughTrackList.h"
#include "MdcHoughFinder/HoughTrack.h"
#include <math.h>
 
#include "EvTimeEvent/RecEsTime.h"
//#include "MdcCalibFunSvc/IMdcCalibFunSvc.h"
#include "McTruth/MdcMcHit.h" 
#include "TrkBase/TrkFit.h"
#include "TrkBase/TrkHitList.h"
#include "TrkBase/TrkExchangePar.h"
#include "TrkFitter/TrkHelixMaker.h"
#include "TrkFitter/TrkCircleMaker.h"
#include "TrkFitter/TrkLineMaker.h"
#include "TrkFitter/TrkHelixFitter.h"
#include "MdcPrintSvc/IMdcPrintSvc.h"
#include "MdcGeom/EntranceAngle.h"
#include "TrackUtil/Helix.h"
#include "GaudiKernel/IPartPropSvc.h"
#include "MdcRecoUtil/Pdt.h"
#include "RawEvent/RawDataUtil.h"
#include "MdcData/MdcHit.h"
#include "MdcTrkRecon/MdcMap.h"
 
#include "TTree.h"
#include "TH2D.h"
 
#include "BesTimerSvc/IBesTimerSvc.h"                                                         
#include "BesTimerSvc/BesTimerSvc.h"                                                          
#include "BesTimerSvc/BesTimer.h"
 
 
 
using namespace std;
using namespace Event;
 
MdcHoughFinder::MdcHoughFinder(const std::string& name, ISvcLocator* pSvcLocator) :
    Algorithm(name, pSvcLocator)
{
    // Declare the properties  
    declareProperty("debug",          m_debug = 0);
    declareProperty("debugMap",       m_debugMap = 0);
    declareProperty("debug2D",          m_debug2D = 0);
    declareProperty("debugTrack",          m_debugTrack = 0);
    declareProperty("debugPeak",          m_debugPeak = 0);
    declareProperty("debugStereo",          m_debugStereo= 0);
    declareProperty("debugZs",          m_debugZs= 0);
    declareProperty("debug3D",          m_debug3D= 0);
    declareProperty("debugArbHit",          m_debugArbHit= 0);
    declareProperty("hist",           m_hist  = 0);
    declareProperty("filter",           m_filter= 0);
    //read raw data setup
    declareProperty("keepBadTdc",     m_keepBadTdc = 0);
    declareProperty("dropHot",        m_dropHot= 0);
    declareProperty("keepUnmatch",    m_keepUnmatch = 0);
    // combine pattsf
    declareProperty("combineTracking",m_combineTracking = false);
    declareProperty("removeBadTrack", m_removeBadTrack = 0);
    declareProperty("dropTrkDrCut",   m_dropTrkDrCut= 1.);
    declareProperty("dropTrkDzCut",   m_dropTrkDzCut= 10.);
    declareProperty("dropTrkPtCut",   m_dropTrkPtCut= 0.12);
    declareProperty("dropTrkChi2Cut", m_dropTrkChi2Cut = 10000.);
    //input setup
    declareProperty("inputType",      m_inputType = 0);
    //set MdcHoughFinder map 
    declareProperty("mapCharge",      m_mapCharge= -1);    //0 use whole ; 1 only half 
    declareProperty("useHalf",        m_useHalf= 0);    //0 use whole ; 1 only half 
    declareProperty("mapHitStyle",    m_mapHit= 0);    //0 : all ; 1 :axial
    declareProperty("nbinTheta",      m_nBinTheta = 100);
    declareProperty("nbinRho",        m_nBinRho = 100);
    declareProperty("rhoRange",       m_rhoRange = 0.1);
    declareProperty("peakWidth",      m_peakWidth= 3);
    declareProperty("peakHigh",       m_peakHigh= 1);
    declareProperty("hitPro",         m_hitPro= 0.4);
 
    declareProperty("recpos",         m_recpos= 1);
    declareProperty("recneg",         m_recneg= 1);
    declareProperty("combineSelf",    m_combine= 1);
    declareProperty("z0CutCompareHough",    m_z0Cut_compareHough= 10);
 
    //split drift circle
    declareProperty("n1",          m_npart= 100);
    declareProperty("n2",          m_n2= 30);
 
    declareProperty("d1",          m_d1= 0.2);
    declareProperty("d2",          m_d2= 0.2);
 
    declareProperty("pdtFile",        m_pdtFile = "pdt.table");
    declareProperty("eventFile",      m_evtFile= "EventList");
 
    declareProperty("cgem",           m_cgem=true);
    declareProperty("skipMDCIT",      m_skipMDCIT=true);
    declareProperty("globalfit",      m_globalfit=true);
    declareProperty("zsRecMethod",    m_recMethod=1);
    declareProperty("factor2D",       m_factor2D=5.0);
    declareProperty("factor3D",       m_factor3D=5.0);
    declareProperty("cut2D",          m_cut_2D);
    declareProperty("cut3D",          m_cut_3D);
    declareProperty("maxGapLength",   m_maxGapLength=99);
    declareProperty("nHitDeleted",    m_nHitDeleted=99);
    declareProperty("nPoint3D",     m_nPoint3D=5);
    declareProperty("nPoint2D",     m_nPoint2D=5);
    declareProperty("nRMS",         m_nRMS=3);
    declareProperty("fillMap",      m_fillMap=1);
    declareProperty("printflag",    m_printflag=0);
    //declareProperty("cgemfit",        m_cgemfit=1); //0:not fit, 1:fit both 2D and 3D  , 2:fit 2D only  , 3:fit 3D only
 
}
//
StatusCode MdcHoughFinder::beginRun(){
    //Initailize MdcDetector
    Global::m_gm = MdcDetector::instance(0);
    if(NULL == Global::m_gm) return StatusCode::FAILURE;
 
    return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode MdcHoughFinder::initialize(){
 
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "in initialize()" << endreq;
 
    StatusCode sc;
 
    //particle
    IPartPropSvc* p_PartPropSvc;
    static const bool CREATEIFNOTTHERE(true);
    sc = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
    if (!sc.isSuccess() || 0 == p_PartPropSvc) {
        log << MSG::ERROR << " Could not initialize PartPropSvc" << endreq;
        return sc;
    }
    m_particleTable = p_PartPropSvc->PDT();
 
    // RawData
    IRawDataProviderSvc* irawDataProviderSvc;
    sc = service ("RawDataProviderSvc", irawDataProviderSvc);
    m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
    if ( sc.isFailure() ){
        log << MSG::FATAL << name()<<" Could not load RawDataProviderSvc!" << endreq;
        return StatusCode::FAILURE;
    }
 
    //  Geometry
    IMdcGeomSvc*     imdcGeomSvc;
    sc = service ("MdcGeomSvc", imdcGeomSvc);
    m_mdcGeomSvc = dynamic_cast<MdcGeomSvc*> (imdcGeomSvc);
    if ( sc.isFailure() ){
        log << MSG::FATAL << "Could not load MdcGeoSvc!" << endreq;
        return StatusCode::FAILURE;
    }
 
    //---CGEM geometry
    if(m_cgem){
        ICgemGeomSvc*    icgemGeomSvc;
        sc = service ("CgemGeomSvc", icgemGeomSvc);
        m_cgemGeomSvc = dynamic_cast<CgemGeomSvc*> (icgemGeomSvc);
        if ( sc.isFailure() ){
            log << MSG::FATAL << "Could not load CgemGeomSvc!" << endreq;
            return StatusCode::FAILURE;
        }
        HoughHit::setCgemGeomSvc(m_cgemGeomSvc);
        //CgemHitOnTrack::cgemGeomSvc = m_cgemGeomSvc;
 
        ICgemCalibFunSvc* icgemCalibSvc;
        sc = service ("CgemCalibFunSvc", icgemCalibSvc);
        m_cgemCalibFunSvc = dynamic_cast<CgemCalibFunSvc*>(icgemCalibSvc);
        if ( sc.isFailure() ){
            log << MSG::FATAL << "Could not load CgemCalibFunSvc!" << endreq;
            return StatusCode::FAILURE;
        }
 
    }
 
    //  magneticfield
    sc = service ("MagneticFieldSvc",m_pIMF);
    if(sc!=StatusCode::SUCCESS) {
        log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
    }
    m_bfield = new BField(m_pIMF);
    log << MSG::INFO << "field z = "<<m_bfield->bFieldNominal()<< endreq;
    m_context = new TrkContextEv(m_bfield);
 
    //read pdt
    Pdt::readMCppTable(m_pdtFile);
 
    //Get MdcCalibFunSvc
    IMdcCalibFunSvc* imdcCalibSvc; 
    sc = service ("MdcCalibFunSvc", imdcCalibSvc);
    m_mdcCalibFunSvc = dynamic_cast<MdcCalibFunSvc*>(imdcCalibSvc);
    if ( sc.isFailure() ){
        log << MSG::FATAL << "Could not load MdcCalibFunSvc!" << endreq;
        return StatusCode::FAILURE;
    }
 
    //initialize MdcPrintSvc
    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;
    }
 
    //time 
    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");
 
    //initialize static 
    HoughHit::setMdcCalibFunSvc(m_mdcCalibFunSvc);
    HoughHit::setMdcGeomSvc(m_mdcGeomSvc);
    HoughHit::setBunchTime(m_bunchT0);
    Hough2D::setContext(m_context);
    Hough2D::setCalib(m_mdcCalibFunSvc);
    Hough2D::setCalib(m_cgemCalibFunSvc);
    Hough2D::setGeomSvc(m_cgemGeomSvc);
    //Hough2D::setcgemfitflag(m_cgemfit);
    Hough3D::setContext(m_context);
    Hough3D::setCalib(m_mdcCalibFunSvc);
    Hough3D::setCalib(m_cgemCalibFunSvc);
    Hough3D::setGeomSvc(m_cgemGeomSvc);
    //Hough3D::setcgemfitflag(m_cgemfit);
    HoughZsFit::setmethod(m_recMethod);
 
    HoughHit::_npart=m_npart;
    HoughMap::m_useHalfCir=m_useHalf;
    HoughMap::m_N1=m_npart;
    HoughMap::m_N2=m_n2;
    HoughMap::m_nPoint2D=m_nPoint2D;
    HoughMap::m_nRMS=m_nRMS;
    HoughMap::m_method=m_fillMap;
 
    HoughPeak::m_dcut1=m_d1;
    HoughPeak::m_dcut2=m_d2;
    HoughPeak::m_factor=m_factor2D;
    HoughPeak::m_cut_deltaD=m_cut_2D;
    HoughZsFit::m_factor=m_factor3D;
    HoughZsFit::m_cut_deltaZ=m_cut_3D;
    HoughZsFit::m_nPoint3D=m_nPoint3D;
 
    if(m_debugMap> 0)   HoughMap::m_debug = 1;
    if(m_debug2D > 0)   Hough2D::m_debug = 1;
    if(m_debug3D > 0)   Hough3D::m_debug = 1;
    if(m_debugTrack> 0) HoughTrack::m_debug = 1;
    if(m_globalfit> 0)HoughTrack::m_globalfit = 1;
    if(m_debugPeak> 0)  HoughPeak::m_debug = 1;
    if(m_debugPeak> 0)  HoughPeak::m_debug = 1;
    if(m_debugStereo> 0)HoughStereo::m_debug = 1;
    if(m_debugZs> 0)    HoughZsFit::m_debug = 1;
    //if(m_recMethod> 0)  HoughZsFit::m_recMethod = 1;
    if(m_debug3D > 4)   TrkHelixFitter::m_debug = 1;
 
    if(m_hist) sc = bookTuple();
    if ( sc.isFailure() ){
        log << MSG::FATAL << "Could not book Tuple !" << endreq;
        return StatusCode::FAILURE;
    }
    return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode MdcHoughFinder::execute() {
 
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "in execute()" << endreq;
    cout.precision(6);
 
 
    m_timer_all->start();
 
    //event start time
    SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
    if (!evTimeCol) {
        //log << MSG::WARNING<< "Could not retrieve RecEsTimeCol , use t0=0" << endreq;
        //m_bunchT0=0.;
        log << MSG::WARNING<< "Could not retrieve RecEsTimeCol" << endreq;
        return StatusCode::SUCCESS;
    }else{
        RecEsTimeCol::iterator iter_evt = evTimeCol->begin();
        if (iter_evt != evTimeCol->end()){
            m_bunchT0 = (*iter_evt)->getTest()*1.e-9;//m_bunchT0-s, getTest-ns
        }
    }
    HoughHit::setBunchTime(m_bunchT0);
 
    //-- Get the event header, print out event and run number
    SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
    if (!eventHeader) {
        log << MSG::FATAL << "Could not find Event Header" << endreq;
        return StatusCode::FAILURE;
    }
 
    //-- Get CGEM cluster
    SmartDataPtr<RecCgemClusterCol> recCgemClusterCol(eventSvc(),  "/Event/Recon/RecCgemClusterCol");
    if(m_cgem){
        if (!recCgemClusterCol){
            log << MSG::WARNING << "Could not retrieve Cgem cluster list" << endreq;
            return StatusCode::FAILURE;
        }
    }
 
    m_event=eventHeader->eventNumber();
    m_run=eventHeader->runNumber();
    if(m_debug>0) cout<<"begin evt "<<eventHeader->eventNumber()<<endl;
    if(m_event%1000==0) cout  <<  "  run No: "<< m_run << "  event No: " << m_event << endl;
 
    //prepare tds 
    RecMdcTrackCol* trackList_tds;
    RecMdcHitCol* hitList_tds;
    vector<RecMdcTrack*> vec_trackPatTds;
    int nTdsTk = storeTracks(trackList_tds,hitList_tds,vec_trackPatTds);
    //print track in pattsf with bad vertex
    if(m_debug>0){
        RecMdcTrackCol::iterator iteritrk_pattsf = vec_trackPatTds.begin();
        for(;iteritrk_pattsf!=vec_trackPatTds.end();iteritrk_pattsf++){
            cout<<"in PATTSF LOST: "<<(*iteritrk_pattsf)->helix(0)<<" "<<(*iteritrk_pattsf)->helix(1)<<" "<<(*iteritrk_pattsf)->helix(2)<<" "<<(*iteritrk_pattsf)->helix(3)<<" "<<(*iteritrk_pattsf)->helix(4)<<" chi2 "<<(*iteritrk_pattsf)->chi2()<<endl;
        }
    }
 
    //for arbi hits
    MdcTrackParams  m_par;
    if(m_debugArbHit>0 ) m_par.lPrint=8;
    m_par.lRemoveInActive=1;
    //m_par.lUseQualCuts=0;
    m_par.maxGapLength=m_maxGapLength;
    m_par.nHitDeleted=m_nHitDeleted;
    //m_par.maxChisq=99;
    //m_par.minHits=99;
 
    // if filter read eventNum in file
    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);
    }
    //m_printflag = 0;
    if(m_printflag==1)cout<<"Begin Event ============================================================================================================== "<<eventHeader->eventNumber()<<endl;
 
    if(m_inputType == -1) GetMcInfo();
 
    //-- Get MDC digi vector
    if(m_debug>0) cout<<"step1 :  prepare digi "<<endl;
    MdcDigiVec mdcDigiVec = prepareDigi();
 
    //-- Create MdcHoughFinder hit list
    bool debugTruth = false;
    if(m_inputType == -1) debugTruth= true;
    if(m_debug>0) cout<<"step2 : hits-> hough hit list "<<endl;
    HoughHitList houghHitList(mdcDigiVec);
    if(m_cgem) houghHitList.addCgemClusterList(recCgemClusterCol); // add CGEM cluster
    //houghHitList.print();
    //if(m_cgem) houghHitList.addCgemClusterList(); // add CGEM cluster
    //  if( mdcDigiVec.size() < 10 ) return StatusCode::SUCCESS;
    if( houghHitList.nHit() < 10 || houghHitList.nHit()>500 ) return StatusCode::SUCCESS;
    if(debugTruth) houghHitList.addTruthInfo(g_tkParTruth);
    if(m_hist) dumpHoughHitList(houghHitList);
    if(m_debug>0) houghHitList.printAll();
 
    vector<MdcHit*>  mdcHitCol_neg;
    vector<MdcHit*>  mdcHitCol_pos;
    MdcDigiVec::iterator iter = mdcDigiVec.begin();
    for (;iter != mdcDigiVec.end(); iter++) {
        const MdcDigi* digi = (*iter);
        if( HoughHit(digi).driftTime()>1000 || HoughHit(digi).driftTime()<=0 ) continue;
        MdcHit *mdcHit= new MdcHit(digi, Global::m_gm);
        MdcHit *mdcHit_pos= new MdcHit(digi, Global::m_gm);
        mdcHitCol_neg.push_back(mdcHit);
        mdcHitCol_pos.push_back(mdcHit_pos);
    }
 
    HoughMap *m_map = new HoughMap(-1,houghHitList,m_mapHit,m_nBinTheta,m_nBinRho,-1.*m_rhoRange,m_rhoRange,m_peakWidth,m_peakHigh,m_hitPro);//, 2dOr3d);
 
    //m_map->printPeak();
    //m_map->printTrack();
    //if(m_hist) m_nHit = houghHitList.nHit();
    //if(m_hist) dumpHoughMap(*m_map);
 
    //track
    if(m_debug>0) cout<<"step3 : neg track list "<<endl;
    vector<HoughTrack*> vec_track_neg; 
    vector<HoughTrack*> vec_track2D_neg; 
    MdcTrackList  mdcTrackList_neg(m_par) ;
    if(m_recneg){
        HoughTrackList trackList_neg(*m_map);
        int trackList_size = trackList_neg.getTrackNum();
        vector< vector<int> > stat_2d(2,vector<int>(trackList_size,0) );
        vector< vector<int> > stat_3d(2,vector<int>(trackList_size,0) );
        //int ifit=0;
        //int ifit3D=0;
        if(m_printflag==1)cout<<"----------------------------------------------------------------------------------------------------"<<endl;
        if(m_printflag==1)cout<<" suppose charge negetive,ntracks: "<<trackList_size<<endl;
        for(int itrack=0;itrack<trackList_size;itrack++){
            if(m_debug>0) cout<<"begin track: "<<itrack<<endl;
            if(m_printflag==1)cout<<endl;
            if(m_printflag==1)cout<<"begin track: "<<itrack<<endl;
            //suppose charge -
            if(m_debug>0) cout<<" suppose charge -1 "<<endl;
            HoughTrack &track_neg = trackList_neg.getTrack(itrack);
            track_neg.setMdcHit( &mdcHitCol_neg);
            track_neg.setHoughHitList(houghHitList.getList());
            track_neg.setCharge(-1);
            track_neg.setMcPar(g_tkParTruth);
            track_neg.setbunchTime(m_bunchT0);
            track_neg.fitLeast();
            if(m_printflag==1)cout<<"---before 2D fitting---"<<endl;
            if(m_printflag==1)track_neg.printHoughTrack();
            stat_2d[0][itrack] = track_neg.fit2D(m_bunchT0);
            int track_charge_2d = track_neg.trackCharge2D();
            if(m_debug>0) cout<<" charge -1 stat2d "<<stat_2d[0][itrack]<<" "<<track_charge_2d<<endl; 
            //cout<<" charge  stat2d "<<stat_2d[0][itrack]<<" "<<track_charge_2d<<endl; 
            if( stat_2d[0][itrack] == 0 || track_charge_2d ==0  ){
                if( stat_2d[0][itrack] == 0 && m_printflag==1)cout<<"==============================> 2D fit fail"<<endl;
                if(track_charge_2d ==0 && m_printflag==1)cout<<"==============================> 2D charge wrong"<<endl;
                continue;
            }
 
            if(m_hist==2)dumpHoughTrack(track_neg,itrack,trackList_size);
            if(m_hist==2)dumpHitOnTrack(track_neg);
            //vec_track2D_neg.push_back( &track_neg );
            //fill 2D inf
            //ifit++;
            //3D fit
            int nHit3d = track_neg.find_stereo_hit();
            int npair= track_neg.find_pair_hit();
            //cout<<"npair "<<npair<<endl;
 
            int track_charge_3d = track_neg.trackCharge3D();
            if(m_debug>0) cout<<" nhitstereo -1 "<<nHit3d<<" "<<track_charge_3d<<endl; 
            //cout<<" nhitstereo    "<<nHit3d<<" "<<track_charge_3d<<endl; 
            if( nHit3d <3 || track_charge_3d==0 ){
                if(nHit3d <3 && m_printflag==1)cout<<"==============================> stereo hit <3"<<endl;
                if(track_charge_3d==0 && m_printflag==1)cout<<"==============================> 3D charge wrong"<<endl;
                continue;
            }
 
            //choose fit method
            if( npair==0 ) stat_3d[0][itrack] = track_neg.fit3D();  
            else stat_3d[0][itrack] = track_neg.fit3D_inner();  
            //else stat_3d[0][itrack] = track_neg.fit3D();  
 
            if(m_debug>0) cout<<" charge -1 stat3d "<<stat_3d[0][itrack]<<" "<<endl;
            //cout<<" charge  stat3d "<<stat_3d[0][itrack]<<" "<<endl;
            if( stat_3d[0][itrack]==0 ){
                if(m_printflag==1)cout<<"==============================> 3D fit fail"<<endl;
                continue;
            }
            vec_track_neg.push_back( &track_neg );
            //fill 3D inf
            if(m_hist==3)dumpHoughTrack(track_neg,itrack,trackList_size);
            if(m_hist==3)dumpHitOnTrack(track_neg);
            //track_neg.printRecHit();
 
        }
 
        std::sort ( vec_track_neg.begin(),vec_track_neg.end(),more_pt); 
        if(m_globalfit){
            if(m_printflag==1)cout<<"tracks before fitting: "<<trackList_neg.getTrackNum()<<endl<<"tracks after  fittng: "<<vec_track_neg.size()<<endl;
            //track for clear
            vector<MdcTrack*> vec_MdcTrack_neg;
            for( unsigned int i =0;i<vec_track_neg.size();i++){
                MdcTrack *mdcTrack = new MdcTrack(vec_track_neg[i]->getTrk());
                vec_MdcTrack_neg.push_back(mdcTrack);
                if(m_debug>0) cout<<"trackneg: "<<i<<" pt "<<vec_track_neg[i]->getPt3D()<<endl;
                if(m_debug>0) vec_track_neg[i]->print();
                if(m_printflag==1)vec_track_neg[i]->getTrk()->printAll(std::cout);
            }
            if(m_debug>0) cout<<"step4 : judge neg track list "<<endl;
            int nDeleted = judgeHit(mdcTrackList_neg,vec_MdcTrack_neg);
            if(m_printflag==1)cout<<"tracks before judge hits: "<<vec_MdcTrack_neg.size()<<endl;
            //printTrack(vec_MdcTrack_neg);
            if(m_printflag==1)cout<<"tracks after  judge hits: "<<mdcTrackList_neg.length()<<endl;
            if(m_printflag==1)printTrack(mdcTrackList_neg);
            //if(nDeleted>0)cout<<nDeleted<<"  tracks have been deleted by MdcTrackList.arbitrateHits()"<<endl;
            if(m_debug>0 && nDeleted>0)cout<<nDeleted<<"  tracks have been deleted by MdcTrackList.arbitrateHits()"<<endl;
            if(m_debug>0) cout<<"finish - charge "<<endl;
        }else{
            storeTracks(trackList_tds, hitList_tds, trackList_neg);
        }
    }
 
    //do rec pos charge
 
    HoughMap *m_map2 = new HoughMap(1,houghHitList,m_mapHit,m_nBinTheta,m_nBinRho,-1.*m_rhoRange,m_rhoRange,m_peakWidth,m_peakHigh,m_hitPro);//, 2dOr3d);
    if(m_debug>0) cout<<"step5 : pos track list "<<endl;
    vector<HoughTrack*> vec_track_pos; 
    MdcTrackList  mdcTrackList_pos(m_par) ;
    if(m_recpos){
        HoughTrackList tracklist_pos(*m_map2);
        int tracklist_size2 = tracklist_pos.getTrackNum();
        vector< vector<int> > stat_2d2(2,vector<int>(tracklist_size2,0) );
        vector< vector<int> > stat_3d2(2,vector<int>(tracklist_size2,0) );
        //int ifit=0;
        //int ifit3D=0;
        if(m_printflag==1)cout<<endl<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
        if(m_printflag==1)cout<<" suppose charge positive,ntracks: "<<tracklist_size2<<endl;
        for(int itrack=0;itrack<tracklist_size2;itrack++){
            if(m_printflag==1)cout<<endl;
            if(m_printflag==1)cout<<"begin track: "<<itrack<<endl;
            //suppose charge +
            if(m_debug>0) cout<<" suppose charge +1 "<<endl;
            HoughTrack &track_pos = tracklist_pos.getTrack(itrack);
            //track_pos.setMdcHit( &mdcHitCol_pos);
            track_pos.setMdcHit( &mdcHitCol_neg);
            track_pos.setHoughHitList(houghHitList.getList());
            track_pos.setCharge(1);
            track_pos.setMcPar(g_tkParTruth);
            track_pos.setbunchTime(m_bunchT0);
            track_pos.fitLeast();
            if(m_printflag==1)cout<<"---before 2D fitting---"<<endl;
            if(m_printflag==1)track_pos.printHoughTrack();
            stat_2d2[0][itrack] = track_pos.fit2D(m_bunchT0);
            int track_charge_2d = track_pos.trackCharge2D();
            if(m_debug>0) cout<<" charge +1 stat2d "<<stat_2d2[0][itrack]<<" "<<track_charge_2d<<endl; 
            if( stat_2d2[0][itrack] == 0 || track_charge_2d==0 ){
                if( stat_2d2[0][itrack] == 0 && m_printflag==1)cout<<"==============================> 2D fit fail"<<endl;
                if(track_charge_2d ==0 && m_printflag==1)cout<<"==============================> 2D charge wrong"<<endl;
                continue;
            }
            //fill 2d inf
            //ifit++;
            //3d fit
            int nHit3d = track_pos.find_stereo_hit();
            int npair= track_pos.find_pair_hit();
            //cout<<"npair "<<npair<<endl;
 
            int track_charge_3d = track_pos.trackCharge3D();
            if(m_debug>0) cout<<" nhitstereo +1 "<<nHit3d<<" "<<track_charge_3d<<endl; 
            if( nHit3d <3 || track_pos.trackCharge3D()==0 ){
                if(nHit3d <3 && m_printflag==1)cout<<"==============================> stereo hit <3"<<endl;
                if(track_charge_3d==0 && m_printflag==1)cout<<"==============================> 3D charge wrong"<<endl;
                continue;
            }
            //choose fit method
            if( npair==0 ) stat_3d2[0][itrack] = track_pos.fit3D();  
            else stat_3d2[0][itrack] = track_pos.fit3D_inner();  
            //else stat_3d2[0][itrack] = track_neg.fit3D();  
            //dumpHitOnTrack(track_pos);
 
            if(m_debug>0) cout<<" charge +1 stat3d "<<stat_3d2[0][itrack]<<" "<<endl;
            if( stat_3d2[0][itrack]==0 ){
                if(m_printflag==1)cout<<"==============================> 3D fit fail"<<endl;
                continue;
            }
            vec_track_pos.push_back( &track_pos );
            //track_pos.printRecHit();
            //fill 3d inf
            //ifit3D++;
            //dumpHoughTrack(track_pos,itrack,tracklist_size2);
        }
 
        //sort pos tracklist
        std::sort ( vec_track_pos.begin(),vec_track_pos.end(),more_pt); 
 
        // clear pos track
        if(m_globalfit){
            if(m_printflag==1)cout<<"tracks before fitting: "<<tracklist_pos.getTrackNum()<<endl<<"tracks after  fittng: "<<vec_track_pos.size()<<endl;
            vector<MdcTrack*> vec_MdcTrack_pos;
            for( unsigned int i =0;i<vec_track_pos.size();i++){
                MdcTrack *mdcTrack = new MdcTrack(vec_track_pos[i]->getTrk());
                vec_MdcTrack_pos.push_back(mdcTrack);
                if(m_debug>0) cout<<"trackpos : "<<i<<" pt "<<vec_track_pos[i]->getPt3D()<<endl;
                if(m_debug>0) vec_track_pos[i]->print();
                if(m_printflag==1)vec_track_pos[i]->getTrk()->printAll(std::cout);
            }
            if(m_debug>0) cout<<"step6 : judge pos track list "<<endl;
            int nDeleted = judgeHit(mdcTrackList_pos,vec_MdcTrack_pos);
            if(m_printflag==1)cout<<"tracks before judge hits: "<<vec_MdcTrack_pos.size()<<endl;
            //printTrack(vec_MdcTrack_pos);
            if(m_printflag==1)cout<<"tracks after  judge hits: "<<mdcTrackList_pos.length()<<endl;
            if(m_printflag==1)printTrack(mdcTrackList_pos);
            //if(nDeleted>0)cout<<nDeleted<<"  tracks have been deleted by MdcTrackList.arbitrateHits()"<<endl;
            if(m_debug>0 && nDeleted>0)cout<<nDeleted<<"  tracks have been deleted by MdcTrackList.arbitrateHits()"<<endl;
 
        }else{
            storeTracks(trackList_tds, hitList_tds, tracklist_pos);
        }
    }
 
    if(m_globalfit){
        //combine hough itself
        //  if(m_debug>0) cout<<"step7 : combine neg&pos track list "<<endl;
        if(m_printflag==1)cout<<endl<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
        if(m_printflag==1)cout<<endl<<"-/+tracks before merge: "<<mdcTrackList_neg.length()<<"  "<<mdcTrackList_pos.length()<<endl;
        if(m_printflag==1)printTrack(mdcTrackList_neg);
        if(m_printflag==1)printTrack(mdcTrackList_pos);
        if(m_combine){
            compareHough(mdcTrackList_neg);
            compareHough(mdcTrackList_pos);
        }
        if(m_printflag==1)cout<<endl<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
        if(m_printflag==1)cout<<"-/+tracks after  merge: "<<mdcTrackList_neg.length()<<"  "<<mdcTrackList_pos.length()<<endl;
        if(m_printflag==1)printTrack(mdcTrackList_neg);
        if(m_printflag==1)printTrack(mdcTrackList_pos);
        if( mdcTrackList_neg.length()!=0 && mdcTrackList_pos.length()!=0 ) judgeChargeTrack(mdcTrackList_neg,mdcTrackList_pos);
        if(m_printflag==1)cout<<endl<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
        if(m_printflag==1)cout<<"tracks after -/+ judge: "<<mdcTrackList_neg.length()<<"  "<<mdcTrackList_pos.length()<<endl;
        if(m_printflag==1)printTrack(mdcTrackList_neg);
        if(m_printflag==1)printTrack(mdcTrackList_pos);
        //cout<<"pos  "<<mdcTrackList_pos.length()<<endl;
        //cout<<"neg  "<<mdcTrackList_neg.length()<<endl;
        //add tracklist
        mdcTrackList_neg+=mdcTrackList_pos;   //neg -> all charge
        //if(m_cgem || m_skipMDCIT)mdcTrackList_neg.setNoInner(true);
        //MdcTrackParams  new_par;
        //new_par.lPrint=8;
        //mdcTrackList_neg.newParams(new_par);
        int nDeleted = mdcTrackList_neg.arbitrateHits();
        if(m_printflag==1)cout<<endl<<"~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~"<<endl;
        if(m_printflag==1 && nDeleted>0)cout<<nDeleted<<"  tracks have been deleted by MdcTrackList.arbitrateHits()"<<endl;
        if(m_printflag==1)cout<<"event:"<<eventHeader->eventNumber()<<"  "<<"total tracks: "<<mdcTrackList_neg.length()<<endl<<endl;
        if(m_printflag==1)printTrack(mdcTrackList_neg);
 
        //compare pattsf&hough self
        if(m_combineTracking) nTdsTk = comparePatTsf(mdcTrackList_neg,trackList_tds);
 
        // store tds
        if(m_debug>0) cout<<"step8 : store tds  "<<endl;
        if(m_debug>0) cout<<"store tds "<<endl;
        int tkId = nTdsTk ;   //trkid
        for( unsigned int i =0;i<mdcTrackList_neg.length();i++){
            if(m_debug>0) cout<<"- charge size i : "<<i<<" "<<mdcTrackList_neg.length()<<endl;
            int tkStat = 4;//track find by houghspace set stat=4
            mdcTrackList_neg[i]->storeTrack(tkId, trackList_tds, hitList_tds, tkStat);
            tkId++;
            delete  mdcTrackList_neg[i];
        }
    }
    if(m_debug>0)  m_mdcPrintSvc->printRecMdcTrackCol();
 
    m_timer_all->stop();
    double t_teventTime = m_timer_all->elapsed();
    //RecMdcTrackCol::iterator iteritrk = trackList_tds->begin();
    //for(;iteritrk!=trackList_tds->end();iteritrk++){
        //cout<<"nlayer:"<<(*iteritrk)->nlayer()<<endl;
    //}
    if(m_hist==1)dumpHoughTrack(trackList_tds);
    if(m_hist==1)dumpHitOnTrack(trackList_tds);
 
    if(m_hist)ntuple_hit->write();
    if(m_hist)ntuple_trk->write();
    delete m_map;
    delete m_map2;
    if(m_debug>0) cout<<"after delete map "<<endl;
    for(int ihit=0;ihit<mdcHitCol_neg.size();ihit++){
        delete mdcHitCol_neg[ihit];
        delete mdcHitCol_pos[ihit];
    }
 
    if(m_debug>0) cout<<"after delete hit"<<endl;
    //clearMem(mdcTrackList_neg,mdcTrackList_pos);
    if(m_debug>0) cout<<"end event "<<endl;
 
    return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode MdcHoughFinder::finalize() {
    MsgStream log(msgSvc(), name());
    delete m_bfield ;
    log << MSG::INFO<< "in finalize()" << endreq;
 
    return StatusCode::SUCCESS;
}
 
int MdcHoughFinder::GetMcInfo(){
    MsgStream log(msgSvc(), name());
    StatusCode sc;
    SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
    if (!mcParticleCol) {
        log << MSG::ERROR << "Could not find McParticle" << endreq;
        return StatusCode::FAILURE;
    }
 
 
    g_tkParTruth.clear();//clear track param.
 
    //t_t0Truth=-1;
    int t_qTruth = 0;
    int t_pidTruth = -999;
    int nMcTk=0;
    McParticleCol::iterator iter_mc = mcParticleCol->begin();
    for (;iter_mc != mcParticleCol->end(); iter_mc++ ) {
        //cout<<"trackid: "<<setw(4)<<(*iter_mc)->trackIndex()<<"  "
            //<<setw(8)<<(*iter_mc)->particleProperty()<<"  "
            //<<(*iter_mc)->primaryParticle()<<"  "
            //<<(*iter_mc)->leafParticle()<<"  "
            //<<(*iter_mc)->decayFromGenerator()<<"  "
            //<<(*iter_mc)->decayInFlight()<<"  "
            //<<endl;
        t_pidTruth = (*iter_mc)->particleProperty();
        //if(!(*iter_mc)->primaryParticle()) continue;
        if((*iter_mc)->primaryParticle()) continue;
        if(!(*iter_mc)->decayFromGenerator()) continue;
        if(!((*iter_mc)->particleProperty() == 211 || (*iter_mc)->particleProperty() == -211 || (*iter_mc)->particleProperty() ==  11 ||(*iter_mc)->particleProperty() == -11))continue;
        //if((m_pid!=0) && (t_pidTruth != m_pid)){ continue; }
        //t_t0Truth=(*iter_mc)->initialPosition().t();
        std::string name;
        int pid = t_pidTruth; 
        if( pid == 0 ) { 
            log << MSG::WARNING << "Wrong particle id" <<endreq;
            continue;
        }else{
            IPartPropSvc* p_PartPropSvc;
            static const bool CREATEIFNOTTHERE(true);
            StatusCode PartPropStatus = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
            if (!PartPropStatus.isSuccess() || 0 == p_PartPropSvc) {
                std::cout<< " Could not initialize Particle Properties Service" << std::endl;
            }
            HepPDT::ParticleDataTable* p_particleTable = p_PartPropSvc->PDT(); 
            if( p_particleTable->particle(pid) ){
                name = p_particleTable->particle(pid)->name();
                t_qTruth = p_particleTable->particle(pid)->charge();
            }else if( p_particleTable->particle(-pid) ){
                name = "anti " + p_particleTable->particle(-pid)->name();
                t_qTruth = (-1)*p_particleTable->particle(-pid)->charge();
            }
        }
        t_q = t_qTruth;
        //cout<<"decay: "<<(*iter_mc)->decayInFlight()<<endl; 
 
        //helix
        if(m_debug>1) std::cout<<__FILE__<<"   "<<__LINE__<<" new helix with pos "<<(*iter_mc)->initialPosition().v()<<" mom "<<(*iter_mc)->initialFourMomentum().v()<<" q "<<t_qTruth<<std::endl;
        Helix mchelix = Helix((*iter_mc)->initialPosition().v(), (*iter_mc)->initialFourMomentum().v(), t_qTruth);//cm
        mchelix.pivot( HepPoint3D(0.,0.,0.) );
        //m_helix = &mchelix;
        m_helix = new Helix((*iter_mc)->initialPosition().v(), (*iter_mc)->initialFourMomentum().v(), t_qTruth);
        m_helix->pivot( HepPoint3D(0.,0.,0.) );
        //m_helix = &(Helix((*iter_mc)->initialPosition().v(), (*iter_mc)->initialFourMomentum().v(), t_qTruth));//cm;
 
        int mcTkId = (*iter_mc)->trackIndex();
 
        pair<int, const HepVector> p(mcTkId,mchelix.a());
        g_tkParTruth.insert(p);
        // exchange to rho theta   test negtive charge
        double rho , theta;
        if( mchelix.phi0()>M_PI) { rho=t_qTruth/fabs(mchelix.radius()); theta = mchelix.phi0()-M_PI;} 
        else {rho =-1.*t_qTruth/fabs(mchelix.radius()) ; theta = mchelix.phi0();} 
        t_d0=-mchelix.dr();
        t_z0=mchelix.dz();
        t_pt=mchelix.pt();
        t_p=mchelix.momentum(0.).mag();
        t_tanl=mchelix.tanl();
        t_cos= mchelix.cosTheta();
 
        //m_pzTruth[nMcTk]=mchelix.momentum(0.).z();
        //m_pidTruth[nMcTk] = t_pidTruth;
        //m_costaTruth[nMcTk] = mchelix.cosTheta();
        //m_ptTruth[nMcTk] = mchelix.pt();
        //m_pTruth[nMcTk] = mchelix.momentum(0.).mag();
        //m_qTruth[nMcTk] = t_qTruth;
        //m_drTruth[nMcTk] = mchelix.dr();
        //m_phi0Truth[nMcTk] = mchelix.phi0();
        //m_omegaTruth[nMcTk]  =1./(mchelix.radius());         //Q
        //m_dzTruth[nMcTk] = mchelix.dz();
        //m_tanl_mc[nMcTk]  =mchelix.tanl();
        //m_rho_mc[nMcTk]  = rho;
        //m_theta_mc[nMcTk]  = theta;
 
 
        double phi0 = mchelix.phi0()+M_PI/2;
        while(phi0> M_PI)phi0-=2*M_PI;
        while(phi0<-M_PI)phi0+=2*M_PI;
 
        if(m_printflag==1){
            cout<<" charge "<<t_qTruth<<endl;
            double rho = fabs(mchelix.kappa()/333.567);
            double theta = atan2((mchelix.dr()-(333.567/mchelix.kappa()))*sin(mchelix.phi0()),(mchelix.dr()-(333.567/mchelix.kappa()))*cos(mchelix.phi0()));
            if(theta<0){theta+=M_PI;rho=-rho;}
            cout<<"(Xc, Yc, R):"<<"("<<(mchelix.dr()-(333.567/mchelix.kappa()))*cos(mchelix.phi0())<<", "<<(mchelix.dr()-(333.567/mchelix.kappa()))*sin(mchelix.phi0())<<", "<<(-333.567/mchelix.kappa())<<")   (theta,rho):"<<"("<<theta<<", "<<rho<<")  ,pt="<<1./mchelix.kappa()<<endl;
            cout<<"truth barbar: "<<setw(15)<<-mchelix.dr()<<setw(15)<<phi0<<setw(15)<<mchelix.kappa()/333.567<<setw(15)<<mchelix.dz()<<setw(15)<<mchelix.tanl()<<endl;
            cout<<"truth besIII: "<<setw(15)<<mchelix.dr()<<setw(15)<<mchelix.phi0()<<setw(15)<<mchelix.kappa()<<setw(15)<<mchelix.dz()<<setw(15)<<mchelix.tanl()<<endl;
            cout<<endl;
        }
        t_phi = phi0;
        t_omega = mchelix.kappa()/333.567;
        //  if(m_debug>0){
        //    std::cout<<"Truth tk "<<nMcTk<<" " <<name<<" pid "<<pid<<" charge "<<t_qTruth << " helix "<< mchelix.a()<<" p "<<mchelix.momentum(0.)<<" p "<<mchelix.momentum(0.).mag()<<" pt "<<mchelix.momentum(0.).perp()<<" pz "<<mchelix.momentum(0.).z()<<std::endl;
        //
        //    cout<<"ptTruth "<<mchelix.pt()<<endl;
        //    cout<<"tanlTruth"<<mchelix.tanl()<<endl;
        //    cout<<"d0Truth"<<mchelix.dr()<<endl;
        //  }
        //  nMcTk++;
    }
    //  m_nTrkMC  = nMcTk;
    //if(m_debug>2) m_mdcPrintSvc->printMdcMcHitCol();
 
    g_firstHit.set(0,0,0);
    SmartDataPtr<MdcMcHitCol> mdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
    if(mdcMcHitCol){
        //cout<<"size:"<<mdcMcHitCol->size()<<endl;
        MdcMcHitCol::iterator iter_mchit = mdcMcHitCol->begin();
        for (;iter_mchit != mdcMcHitCol->end(); iter_mchit++) {
            const Identifier id= (*iter_mchit)->identify();
            int l = MdcID::layer(id);
            int w = MdcID::wire(id);
            //cout<<"MDC trackID: "<<(*iter_mchit)->getTrackIndex()<<"   layer: "<<l<<"   wire: "<<w<<endl;
            if((*iter_mchit)->getTrackIndex()==0) {
                g_firstHit.setX((*iter_mchit)->getPositionX()/10.);//mm to cm
                g_firstHit.setY((*iter_mchit)->getPositionY()/10.);
                g_firstHit.setZ((*iter_mchit)->getPositionZ()/10.);
                //break;
            }
        }
    }else{
        std::cout<<__FILE__<<" Could not get MdcMcHitCol  "<<std::endl;
        return StatusCode::FAILURE;
    }
 
    //add truth info. to HoughHit in list
    return StatusCode::SUCCESS;
}
 
MdcDigiVec MdcHoughFinder::prepareDigi(){
    //  retrieve Mdc digi vector form RawDataProviderSvc 
    uint32_t getDigiFlag = 0;
    if(m_dropHot || m_combineTracking )getDigiFlag |= MdcRawDataProvider::b_dropHot;
    if(m_keepBadTdc)  getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
    if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
 
    MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
    return mdcDigiVec;
}
 
void MdcHoughFinder::dumpHoughHitList(const HoughHitList& houghhitlist){
 
    m_hit_run = m_run;
    m_hit_evt = m_event;
    m_hit_nhit = houghhitlist.nHit();
    std::vector<HoughHit>::const_iterator iter = houghhitlist.getList().begin();
    for(int iHit=0;iter!= houghhitlist.getList().end(); iter++,iHit++){
        const HoughHit h = (*iter);
        if( m_debug>0 ) h.printTruth();
        m_hit_hitid[iHit]       = h.getHitId();
        m_hit_layer[iHit]       = h.getLayerId();
        m_hit_wire[iHit]        = h.getWireId();
        m_hit_x[iHit]           = h.getMidX();
        m_hit_y[iHit]           = h.getMidY();
        m_hit_z[iHit]           = h.getMidPoint().z();
        m_hit_driftdist[iHit]   = h.getDriftDist();
        m_hit_drifttime[iHit]   = h.getDriftTime();
        m_hit_flag[iHit]        = h.getSlayerType();
        m_hit_truth_x[iHit]     = h.getXTruth();
        m_hit_truth_y[iHit]     = h.getYTruth();
        m_hit_truth_z[iHit]     = h.getZTruth();
        m_hit_truth_drift[iHit] = h.getDriftDistTruth();
        m_hit_truth_ambig[iHit] = h.getLrTruth();
    }
    //  m_nHit = houghhitlist.nHit();
}
void MdcHoughFinder::dumpHitOnTrack( HoughTrack& houghTrack){
    std::vector<HoughRecHit> recHitVec = houghTrack.getHoughHitList();
    m_hot_run = m_run;
    m_hot_evt = m_event;
    m_hot_trk = houghTrack.getTrkid();
    m_hot_nhot = recHitVec.size();
    std::vector<HoughRecHit>::const_iterator iter = recHitVec.begin();
    for(int iHit=0;iter!= recHitVec.end(); iter++,iHit++){
        const HoughRecHit h = (*iter);
        if( m_debug>0 ) h.printTruth();
        m_hot_hitid[iHit] = h.getHitId();
        m_hot_layer[iHit] = h.getLayerId();
        m_hot_wire[iHit]  = h.getWireId();
        m_hot_x[iHit]     = h.getMidX();
        m_hot_y[iHit]     = h.getMidY();
        m_hot_z[iHit]     = h.getMidPoint().z();
        m_hot_x0[iHit]    = h.getLeftPoint().x();
        m_hot_y0[iHit]    = h.getLeftPoint().y();
        m_hot_z0[iHit]    = h.getLeftPoint().z();
        m_hot_s0[iHit]    = h.getsAmb(0);
        m_hot_x1[iHit]    = h.getRightPoint().x();
        m_hot_y1[iHit]    = h.getRightPoint().y();
        m_hot_z1[iHit]    = h.getRightPoint().z();
        m_hot_s1[iHit]    = h.getsAmb(1);
        m_hot_drift[iHit] = h.getDriftDist();
        m_hot_flag[iHit]  = h.getSlayerType();
        m_hot_deltaD[iHit]= h.getDeltaD();
        m_hot_truth_x[iHit]  = h.getXTruth();
        m_hot_truth_y[iHit]  = h.getYTruth();
        m_hot_truth_z[iHit]  = h.getZTruth();
        m_hot_truth_drift[iHit] = h.getDriftDistTruth();
        m_hot_truth_ambig[iHit] = h.getLrTruth();
    }
    if(m_hist)ntuple_hot->write();
}
void MdcHoughFinder::dumpHitOnTrack( RecMdcTrackCol* trackList_tds){
    m_hot_run = m_run;
    m_hot_evt = m_event;
    RecMdcTrackCol::iterator iter = trackList_tds->begin();
    for(;iter!=trackList_tds->end();iter++){
        m_hot_trk = (*iter)->trackId();
        m_hot_nhot = ((*iter)->getVecHits()).size()+(*iter)->getNcluster();
        
        int iHit=0;
        ClusterRefVec clusterRefVec = (*iter)->getVecClusters();
        ClusterRefVec::iterator itCluster = clusterRefVec.begin();
        for( ; itCluster != clusterRefVec.end(); itCluster++,iHit++){
            int layer = (*itCluster)->getlayerid();
            int wire  = (*itCluster)->getsheetid();
            m_hot_hitid[iHit] = (*itCluster)->getclusterid();
            m_hot_layer[iHit] = layer;
            m_hot_wire[iHit]  = wire;
            m_hot_x[iHit]     = 0;
            m_hot_y[iHit]     = 0;
            m_hot_z[iHit]     = (*itCluster)->getRecZ();
            m_hot_x0[iHit]    = 0;
            m_hot_y0[iHit]    = 0;
            m_hot_z0[iHit]    = 0;
            m_hot_s0[iHit]    = 0;
            m_hot_x1[iHit]    = 0;
            m_hot_y1[iHit]    = 0;
            m_hot_z1[iHit]    = 0;
            m_hot_s1[iHit]    = 0;
            m_hot_drift[iHit] = 0;
            m_hot_flag[iHit]  = (*itCluster)->getflag();
            m_hot_deltaD[iHit]= 0;
            m_hot_truth_x[iHit]  = 0; 
            m_hot_truth_y[iHit]  = 0;
            m_hot_truth_z[iHit]  = 0;
            m_hot_truth_drift[iHit] = 0;
            m_hot_truth_ambig[iHit] = 0;
        }
 
        HitRefVec vechits = (*iter)->getVecHits();
        HitRefVec::iterator itHit = vechits.begin();
        for(; itHit != vechits.end(); itHit++,iHit++){
            const Identifier mdcid = (*itHit)->getMdcId();
            int layer = MdcID::layer(mdcid);
            int wire  = MdcID::wire(mdcid);
            m_hot_hitid[iHit] = (*itHit)->getId();
            m_hot_layer[iHit] = layer;
            m_hot_wire[iHit]  = wire;
            m_hot_x[iHit]     = 0;
            m_hot_y[iHit]     = 0;
            m_hot_z[iHit]     = (*itHit)->getZhit();
            m_hot_x0[iHit]    = 0;
            m_hot_y0[iHit]    = 0;
            m_hot_z0[iHit]    = 0;
            m_hot_s0[iHit]    = 0;
            m_hot_x1[iHit]    = 0;
            m_hot_y1[iHit]    = 0;
            m_hot_z1[iHit]    = 0;
            m_hot_s1[iHit]    = 0;
            m_hot_drift[iHit] = 0;
            m_hot_flag[iHit]  = (*itHit)->getStat();
            m_hot_deltaD[iHit]= 0;
            m_hot_truth_x[iHit]  = 0; 
            m_hot_truth_y[iHit]  = 0;
            m_hot_truth_z[iHit]  = 0;
            m_hot_truth_drift[iHit] = 0;
            m_hot_truth_ambig[iHit] = 0;
        }
        if(m_hist)ntuple_hot->write();
    }
}
void MdcHoughFinder::dumpHoughTrack(RecMdcTrackCol* trackList_tds){
    m_trk_run  = m_run;
    m_trk_evt  = m_event;
    m_trk_ntrk = trackList_tds->size();
    m_trk_size = trackList_tds->size();
    double alpha = -333.567;
    RecMdcTrackCol::iterator iter = trackList_tds->begin();
    for(int itrack=0;iter!=trackList_tds->end();iter++,itrack++){
        m_trk_trackId[itrack]           = (*iter)->trackId()        ;
        m_trk_charge[itrack]            = (*iter)->charge()         ;
        m_trk_dr[itrack]                = (*iter)->helix(0)         ;
        m_trk_phi0[itrack]              = (*iter)->helix(1)         ;
        m_trk_kappa[itrack]             = (*iter)->helix(2)         ;
        m_trk_dz[itrack]                = (*iter)->helix(3)         ;
        m_trk_tanl[itrack]              = (*iter)->helix(4)         ;
        m_trk_pxy[itrack]               = (*iter)->pxy()            ;
        m_trk_px[itrack]                = (*iter)->px()             ;
        m_trk_py[itrack]                = (*iter)->py()             ;
        m_trk_pz[itrack]                = (*iter)->pz()             ;
        m_trk_p[itrack]                 = (*iter)->p()              ;
        m_trk_theta[itrack]             = (*iter)->theta()          ;
        m_trk_phi[itrack]               = (*iter)->phi()            ;
        m_trk_x[itrack]                 = (*iter)->x()              ;
        m_trk_y[itrack]                 = (*iter)->y()              ;
        m_trk_z[itrack]                 = (*iter)->z()              ;
        m_trk_r[itrack]                 = (*iter)->r()              ;
        m_trk_chi2[itrack]              = (*iter)->chi2()           ;
        m_trk_fiTerm[itrack]            = (*iter)->getFiTerm()      ;
        //m_trk_matchChi2[itrack]         = (*iter)->getMatchChi2()   ;
        m_trk_nhit[itrack]              = (*iter)->getNhits()       ;
        m_trk_ncluster[itrack]          = (*iter)->getNcluster()    ;
        m_trk_stat[itrack]              = (*iter)->stat()           ;
        m_trk_ndof[itrack]              = (*iter)->ndof()           ;
        m_trk_nster[itrack]             = (*iter)-> nster()         ;
        m_trk_nlayer[itrack]            = (*iter)->nlayer()         ;
        m_trk_firstLayer[itrack]        = (*iter)->firstLayer()     ;
        m_trk_lastLayer[itrack]         = (*iter)->lastLayer()      ;
        //m_trk_nCgemXClusters[itrack]    = (*iter)->nCgemXCluster()  ;
        //m_trk_nCgemVClusters[itrack]    = (*iter)->nCgemVCluster()  ;
        m_trk_nhop[itrack]              = 0                         ;
        m_trk_nhot[itrack]              = 0                         ;
        m_trk_Xc[itrack]                = (*iter)->x() + ( (*iter)->helix(0) + alpha/(*iter)->helix(2) )*cos((*iter)->helix(1));
        m_trk_Yc[itrack]                = (*iter)->y() + ( (*iter)->helix(0) + alpha/(*iter)->helix(2) )*sin((*iter)->helix(1));
        m_trk_R[itrack]                 = fabs(alpha/(*iter)->helix(2));
 
        m_trk_truth_charge[itrack]   = t_q;
        m_trk_truth_dr[itrack]       = m_helix->dr();
        m_trk_truth_phi0[itrack]     = m_helix->phi0();
        m_trk_truth_kappa[itrack]    = m_helix->kappa();
        m_trk_truth_dz[itrack]       = m_helix->dz();
        m_trk_truth_tanl[itrack]     = m_helix->tanl();
        m_trk_truth_pxy[itrack]      = m_helix->pt();
        m_trk_truth_px[itrack]       = m_helix->momentum(0.).x();
        m_trk_truth_py[itrack]       = m_helix->momentum(0.).y();
        m_trk_truth_pz[itrack]       = m_helix->momentum(0.).z();
        m_trk_truth_p[itrack]        = m_helix->momentum(0.).mag();
        m_trk_truth_theta[itrack]    = acos(m_helix->momentum(0.).z()/m_helix->momentum(0.).mag());
        m_trk_truth_phi[itrack]      = atan2(m_helix->momentum(0.).y(),m_helix->momentum(0.).x());
        m_trk_truth_x[itrack]        = m_helix->pivot().x();
        m_trk_truth_y[itrack]        = m_helix->pivot().y();
        m_trk_truth_z[itrack]        = m_helix->pivot().z();
        m_trk_truth_r[itrack]        = sqrt(m_helix->pivot().x()*m_helix->pivot().x()+m_helix->pivot().y()*m_helix->pivot().y());
        m_trk_truth_cosTheta[itrack] = m_helix->cosTheta();
        //m_trk_truth_Xc[itrack]       = m_helix->center().x();
        //m_trk_truth_Yc[itrack]       = m_helix->center().y() ;
        //m_trk_truth_R[itrack]        = m_helix->radius();
        m_trk_truth_Xc[itrack]       = m_helix->pivot().x() + ( m_helix->dr() + alpha/m_helix->kappa() )*cos(m_helix->phi0());
        m_trk_truth_Yc[itrack]       = m_helix->pivot().y() + ( m_helix->dr() + alpha/m_helix->kappa() )*sin(m_helix->phi0());
        m_trk_truth_R[itrack]        = fabs(alpha/m_helix->kappa());
        //cout<<(*iter)->helix(0)<<endl;
        //cout<<(*iter)->helix(1)<<endl;
        //cout<<(*iter)->helix(2)<<endl;
        //cout<<(*iter)->helix(3)<<endl;
        //cout<<(*iter)->helix(4)<<endl;
        //cout<<endl;
    }
    if(m_trk_ntrk==0){
        int itrack = 0;
        m_trk_size = 1;
        m_trk_truth_charge[itrack]   = t_q;
        m_trk_truth_dr[itrack]       = m_helix->dr();
        m_trk_truth_phi0[itrack]     = m_helix->phi0();
        m_trk_truth_kappa[itrack]    = m_helix->kappa();
        m_trk_truth_dz[itrack]       = m_helix->dz();
        m_trk_truth_tanl[itrack]     = m_helix->tanl();
        m_trk_truth_pxy[itrack]      = m_helix->pt();
        m_trk_truth_px[itrack]       = m_helix->momentum(0.).x();
        m_trk_truth_py[itrack]       = m_helix->momentum(0.).y();
        m_trk_truth_pz[itrack]       = m_helix->momentum(0.).z();
        m_trk_truth_p[itrack]        = m_helix->momentum(0.).mag();
        m_trk_truth_theta[itrack]    = acos(m_helix->momentum(0.).z()/m_helix->momentum(0.).mag());
        m_trk_truth_phi[itrack]      = atan2(m_helix->momentum(0.).y(),m_helix->momentum(0.).x());
        m_trk_truth_x[itrack]        = m_helix->pivot().x();
        m_trk_truth_y[itrack]        = m_helix->pivot().y();
        m_trk_truth_z[itrack]        = m_helix->pivot().z();
        m_trk_truth_r[itrack]        = sqrt(m_helix->pivot().x()*m_helix->pivot().x()+m_helix->pivot().y()*m_helix->pivot().y());
        m_trk_truth_cosTheta[itrack] = m_helix->cosTheta();
        //m_trk_truth_Xc[itrack]       = m_helix->center().x();
        //m_trk_truth_Yc[itrack]       = m_helix->center().y() ;
        //m_trk_truth_R[itrack]        = m_helix->radius();
        m_trk_truth_Xc[itrack]       = m_helix->pivot().x() + ( m_helix->dr() + alpha/m_helix->kappa() )*cos(m_helix->phi0());
        m_trk_truth_Yc[itrack]       = m_helix->pivot().y() + ( m_helix->dr() + alpha/m_helix->kappa() )*sin(m_helix->phi0());
        m_trk_truth_R[itrack]        = fabs(alpha/m_helix->kappa());
    }
}
void MdcHoughFinder::dumpHoughTrack( HoughTrack& houghTrack, int itrack, int ntrack){
    double alpha = -333.567;
    HepVector barbar(5,0);
    double  d0             = houghTrack.getD0();
    double  phi0           = houghTrack.getPhi0();
    double  omega          = houghTrack.getOmega();
    double  z0             = houghTrack.getZ0();
    double  tanl           = houghTrack.getTanl();
    barbar(1)              = d0;
    barbar(2)              = phi0;
    barbar(3)              = omega;
    barbar(4)              = z0;
    barbar(5)              = tanl;
    HepVector bes          = barbar2bes(barbar);
 
    m_trk_run  = m_run;
    m_trk_evt  = m_event;
    m_trk_ntrk = ntrack;
    m_trk_size = ntrack;
    m_trk_trackId[itrack]  = houghTrack.getTrkid();
    m_trk_charge[itrack]   = houghTrack.getCharge();
    m_trk_dr[itrack]       = bes(1);
    m_trk_phi0[itrack]     = bes(2);
    m_trk_kappa[itrack]    = bes(3);
    m_trk_dz[itrack]       = bes(4);
    m_trk_tanl[itrack]     = bes(5);
    m_trk_pxy[itrack]      = houghTrack.getPt_least();
    m_trk_nhop[itrack]     = houghTrack.getCenterPeak().peakHeight();
    m_trk_nhot[itrack]     = houghTrack.getHitNumA(0);
    m_trk_Xc[itrack]       = houghTrack.getCirX();
    m_trk_Yc[itrack]       = houghTrack.getCirY();
    m_trk_R[itrack]        = houghTrack.getCirR();
 
    m_trk_truth_charge[itrack]   = t_q;
    m_trk_truth_dr[itrack]       = m_helix->dr();
    m_trk_truth_phi0[itrack]     = m_helix->phi0();
    m_trk_truth_kappa[itrack]    = m_helix->kappa();
    m_trk_truth_dz[itrack]       = m_helix->dz();
    m_trk_truth_tanl[itrack]     = m_helix->tanl();
    m_trk_truth_pxy[itrack]      = m_helix->pt();
    m_trk_truth_px[itrack]       = m_helix->momentum(0.).x();
    m_trk_truth_py[itrack]       = m_helix->momentum(0.).y();
    m_trk_truth_pz[itrack]       = m_helix->momentum(0.).z();
    m_trk_truth_p[itrack]        = m_helix->momentum(0.).mag();
    m_trk_truth_theta[itrack]    = acos(m_helix->momentum(0.).z()/m_helix->momentum(0.).mag());
    m_trk_truth_phi[itrack]      = atan2(m_helix->momentum(0.).y(),m_helix->momentum(0.).x());
    m_trk_truth_x[itrack]        = m_helix->pivot().x();
    m_trk_truth_y[itrack]        = m_helix->pivot().y();
    m_trk_truth_z[itrack]        = m_helix->pivot().z();
    m_trk_truth_r[itrack]        = sqrt(m_helix->pivot().x()*m_helix->pivot().x()+m_helix->pivot().y()*m_helix->pivot().y());
    m_trk_truth_cosTheta[itrack] = m_helix->cosTheta();
    //m_trk_truth_Xc[itrack]       = m_helix->center().x();
    //m_trk_truth_Yc[itrack]       = m_helix->center().y() ;
    //m_trk_truth_R[itrack]        = m_helix->radius();
    m_trk_truth_Xc[itrack]       = m_helix->pivot().x() + ( m_helix->dr() + alpha/m_helix->kappa() )*cos(m_helix->phi0());
    m_trk_truth_Yc[itrack]       = m_helix->pivot().y() + ( m_helix->dr() + alpha/m_helix->kappa() )*sin(m_helix->phi0());
    m_trk_truth_R[itrack]        = fabs(alpha/m_helix->kappa());
    //}
}
HepVector barbar2bes(HepVector a)
{
    double  dr    = -a(1);
    double  phi0  = a(2)-M_PI/2;
    while(phi0>2*M_PI)phi0-=2*M_PI;
    while(phi0<0)phi0+=2*M_PI;
    double  kappa = a(3)*333.567;
    double  dz    = a(4);
    double  tanl  = a(5);
    HepVector bes(5,0);
    bes(1) = dr;
    bes(2) = phi0;
    bes(3) = kappa;
    bes(4) = dz;
    bes(5) = tanl;
    return bes;
}
HepVector bes2barbar(HepVector a)
{
    double  d0    = -a(1);
    double  phi0  = a(2)+M_PI/2;
    while(phi0> M_PI)phi0-=2*M_PI;
    while(phi0<-M_PI)phi0+=2*M_PI;
    double  omega = a(3)/333.567;
    double  z0    = a(4);
    double  tanl  = a(5);
    HepVector barbar(5,0);
    barbar(1) = d0;
    barbar(2) = phi0;
    barbar(3) = omega;
    barbar(4) = z0;
    barbar(5) = tanl;
    return barbar;
}
/*
void MdcHoughFinder::diffMap(const HoughMap& houghmap,const HoughMap& houghmap2){
    //map 1
    m_nMap1Pk=houghmap.getPeakNumber();
    for(int iPk=0;iPk< houghmap.getPeakNumber(); iPk++){
        m_PkRho1[iPk]  =  houghmap.getPeak(iPk).getRho();
        m_PkTheta1[iPk]=  houghmap.getPeak(iPk).getTheta();
    }
    m_nMap1Tk=houghmap.getTrackNumber();
    for(int iTk=0;iTk< houghmap.getTrackNumber(); iTk++){
        m_TkRho1[iTk]  =  houghmap.getTrack(iTk).getRho();
        m_TkTheta1[iTk]=  houghmap.getTrack(iTk).getTheta();
    }
    //map 2
    m_nMap2Pk=houghmap2.getPeakNumber();
    for(int iPk=0;iPk< houghmap2.getPeakNumber(); iPk++){
        m_PkRho2[iPk]  =  houghmap2.getPeak(iPk).getRho();
        m_PkTheta2[iPk]=  houghmap2.getPeak(iPk).getTheta();
    }
    m_nMap2Tk=houghmap2.getTrackNumber();
    for(int iTk=0;iTk< houghmap2.getTrackNumber(); iTk++){
        m_TkRho2[iTk]  =  houghmap2.getTrack(iTk).getRho();
        m_TkTheta2[iTk]=  houghmap2.getTrack(iTk).getTheta();
    }
}
 
void MdcHoughFinder::Leastfit(vector<double> u,vector<double> v,double& k ,double& b){
    int N = u.size();
    if( N <= 2 )  return;
    double *x=new double[N];
    double *y=new double[N];
    for(int i=0;i<N;i++){
        x[i]=u[i];
        y[i]=v[i];
    }
 
    TF1 *fpol1=new TF1("fpol1","pol1",-50,50);
    TGraph *tg=new TGraph(N,x,y);
    tg->Fit("fpol1","QN");
    double ktemp  =fpol1->GetParameter(1);
    double btemp  =fpol1->GetParameter(0);
    k = ktemp;
    b = btemp;
    delete []x;
    delete []y;
    delete fpol1;
    delete tg;
}
*/
 
bool more_pt(const HoughTrack* tracka,const HoughTrack* trackb){
    return tracka->getPt3D() > trackb->getPt3D();
}
int MdcHoughFinder::storeTracks(RecMdcTrackCol*& trackList_tds ,RecMdcHitCol*& hitList_tds,vector<RecMdcTrack*>& vec_trackPatTds){
    MsgStream log(msgSvc(), name());
    StatusCode sc;
    // tds
    DataObject *aTrackCol;
    DataObject *aRecHitCol;
    SmartIF<IDataManagerSvc> dataManSvc(eventSvc()); 
    if(!m_combineTracking){
        eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
        if(aTrackCol != NULL) {
            dataManSvc->clearSubTree("/Event/Recon/RecMdcTrackCol");
            eventSvc()->unregisterObject("/Event/Recon/RecMdcTrackCol");
        }
        eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
        if(aRecHitCol != NULL) {
            dataManSvc->clearSubTree("/Event/Recon/RecMdcHitCol");
            eventSvc()->unregisterObject("/Event/Recon/RecMdcHitCol");
        }
    }
 
    //RecMdcTrackCol* trackList_tds;
    eventSvc()->findObject("/Event/Recon/RecMdcTrackCol",aTrackCol);
    if (aTrackCol) {
        trackList_tds = dynamic_cast<RecMdcTrackCol*> (aTrackCol);
    }else{
        trackList_tds = new RecMdcTrackCol;
        StatusCode sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcTrackCol, trackList_tds);
        if(!sc.isSuccess()) {
            log << MSG::FATAL << " Could not register RecMdcTrack collection" <<endreq;
            return StatusCode::FAILURE;
        }
    }
    //RecMdcHitCol* hitList_tds;
    eventSvc()->findObject("/Event/Recon/RecMdcHitCol",aRecHitCol);
    if (aRecHitCol) {
        hitList_tds = dynamic_cast<RecMdcHitCol*> (aRecHitCol);
    }else{
        hitList_tds = new RecMdcHitCol;
        StatusCode sc =  eventSvc()->registerObject(EventModel::Recon::RecMdcHitCol, hitList_tds);
        if(!sc.isSuccess()) {
            log << MSG::FATAL << " Could not register RecMdcHit collection" <<endreq;
            return StatusCode::FAILURE;
        }
    }
    //remove bad quality or low pt track(s)
    if(m_removeBadTrack && m_combineTracking) {
        if( m_debug>0 ) cout<<"PATTSF collect "<<trackList_tds->size()<<" track. "<<endl;
        if(trackList_tds->size()!=0){
            RecMdcTrackCol::iterator iter_pat = trackList_tds->begin();
            for(;iter_pat!=trackList_tds->end();iter_pat++){
                double d0=(*iter_pat)->helix(0);
                double kap = (*iter_pat)->helix(2);
                double pt = 0.00001;
                if(fabs(kap)>0.00001) pt = fabs(1./kap);
                double dz=(*iter_pat)->helix(3);
                double chi2=(*iter_pat)->chi2();
                if( m_debug>0) cout<<"d0: "<<d0<<" z0: "<<dz<<" pt "<<pt<<" chi2 "<<chi2<<endl;
                //if( (fabs(d0)>m_dropTrkDrCut || fabs(dz)>m_dropTrkDzCut || chi2>m_dropTrkChi2Cut) && pt<=m_dropTrkPtCut) 
                if(  pt<=m_dropTrkPtCut ){ 
                    vec_trackPatTds.push_back(*iter_pat);
                    //remove hits on track
                    HitRefVec rechits = (*iter_pat)->getVecHits();
                    HitRefVec::iterator hotIter = rechits.begin();
                    while (hotIter!=rechits.end()) {
                        Identifier id = (*hotIter)->getMdcId();
                        int layer = MdcID::layer(id);
                        int wire  = MdcID::wire(id);
                        if(m_debug>0) cout <<"remove hit " << (*hotIter)->getId()<<"("<<layer<<","<<wire<<")"<<endl;
                        hitList_tds->remove(*hotIter);
                        hotIter++;
                    }
                    //if( (fabs(d0)>3*m_dropTrkDrCut || fabs(dz)>3*m_dropTrkDzCut) ){  // drop pattsf fate track
                    //  trackList_tds->remove(*iter_pat);
                    //  iter_pat--;
                    //}
                }
                if( m_debug>0 ) cout<<"after delete trackcol size : "<<trackList_tds->size()<<endl;
            }
        } else {
            if(m_debug>0) cout<<" PATTSF find 0 track. "<<endl;
        }
    }//end of remove bad quality high pt track
 
    int nTdsTk=trackList_tds->size();
    return nTdsTk;
}//end of stroeTracks
 
void MdcHoughFinder::clearMem(MdcTrackList& list1,MdcTrackList& list2){
    cout<<"length in clearMem "<<list1.length()<<" "<<list2.length()<<endl;
    cout<<"in list "<<vectrk_for_clean.size()<<endl;
 
    for(unsigned int i=0;i<vectrk_for_clean.size();i++){
        int sametrk=0;
        for(unsigned int j=0;j<list1.length();j++){
            cout<<"-i j "<<i<<" "<<j<<" "<<vectrk_for_clean[i]<<" "<<&(list1[j]->track())<<endl;
            if(vectrk_for_clean[i] == &(list1[j]->track()) )  {sametrk=1;cout<<"not delete new trk "<<endl;}
            //delete list1[j];
        }
        for(unsigned int k=0;k<list2.length();k++){
            cout<<"+i k "<<i<<" "<<k<<" "<<vectrk_for_clean[i]<<" "<<&(list2[k]->track())<<endl;
            if(vectrk_for_clean[i] == &(list2[k]->track()) )  {sametrk=1;cout<<"not delete new trk "<<endl;}
            //delete list2[k];
        }
        if( sametrk ==0 ) {
            cout<<"delete i "<<i<<endl;
            delete vectrk_for_clean[i];
            vectrk_for_clean[i] = NULL;
        }
    }
    vectrk_for_clean.clear();
}
 
int MdcHoughFinder::judgeChargeTrack(MdcTrackList& list1 ,MdcTrackList& list2){
    if(m_debug>0) cout<<"in judgeChargeTrack"<<endl;
    if(m_debug>0) cout<<"ntrack length neg  : "<<list1.length()<<endl;
    if(m_debug>0) cout<<"ntrack length pos  : "<<list2.length()<<endl;
    double R_cut = 81.0;//radius of MDC
    for(int itrack1=0; itrack1<list1.nTrack(); itrack1++){
        MdcTrack* track_1 = list1[itrack1];
        TrkExchangePar par1 = track_1->track().fitResult()->helix(0.);
        double d0_1 = par1.d0();
        double phi0_1 = par1.phi0();
        double omega_1 = par1.omega();
        double cr=fabs(1./par1.omega());
        double cx= sin(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1.; //z axis verse,x_babar = -x_belle
        double cy= -1.*cos(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1;//???
        double z0_1 = par1.z0();
        if(d0_1+2.0*cr>R_cut)continue;
        for(int itrack2=0; itrack2<list2.nTrack(); itrack2++){
            MdcTrack* track_2 = list2[itrack2];
            TrkExchangePar par2 = track_2->track().fitResult()->helix(0.);
            double d0_2 = par2.d0();
            double phi0_2 = par2.phi0();
            double omega_2 = par2.omega();
            double cR=fabs(1./par2.omega());
            double cX= sin(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1.; //z axis verse,x_babar = -x_belle
            double cY= -1.*cos(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1;//???
            double z0_2= par2.z0();
            if(d0_2+2.0*cR>R_cut)continue;
            double   bestDiff = 1.0e+20;
            if(fabs(cr)*(1.-0.25) <= fabs(cR) && fabs(cR) <= fabs(cr)*(1.+0.25)){
                if(fabs(cx-cX) <= fabs(cr)*(0.25)&& fabs(cy-cY) <= fabs(cr)*(0.25) ){
                    double diff = fabs((fabs(cr)-fabs(cR))*(cx-cX)*(cy-cY));
                    if(diff < bestDiff){
                        bestDiff = diff;
                    }
                }
            }
 
            if(bestDiff != 1.0e20) { 
                if(m_debug>0) {
                    cout<<"There is ambiguous +- charge in the track list "<<endl;
                    cout<<"z0_1 : "<<z0_1<<" z0_2 : "<<z0_2<<endl;
                }
                //iter_hough2 = vec_track.erase(iter_hough2); iter_hough2--;
                if( fabs(z0_1) >= fabs(z0_2) ){
                    //cout<<"remove track:"<<track_1->track().id()<<"  "<<(track_1->track().fitResult()->pt())*(track_1->track().fitResult()->charge())<<endl;
                    list1.remove( track_1 ); 
                    itrack1--; 
                    break;
                }
                if( fabs(z0_1) < fabs(z0_2) ){
                    //cout<<"remove track:"<<track_2->track().id()<<"  "<<(track_2->track().fitResult()->pt())*(track_2->track().fitResult()->charge())<<endl;
                    list2.remove( track_2 ); 
                    itrack2--; 
                    break;
                }
            }
            if(bestDiff == 1.0e20) {  if(m_debug>0) cout<<" no (+-) track in hough"<<endl; }
        }
    }
}
 
void MdcHoughFinder::compareHough(MdcTrackList& mdcTrackList){
    if(m_debug>0)  cout<<"in compareHough "<<endl;
    if(m_debug>0) cout<<"ntrack  : "<<mdcTrackList.length()<<endl;
    for(int itrack=0; itrack<mdcTrackList.length(); itrack++){
        MdcTrack* track = mdcTrackList[itrack];
        TrkExchangePar par = track->track().fitResult()->helix(0.);
        int nhit = track->track().hots()->nHit();
        double pt = (1./par.omega())/333.567;  
        if(m_debug>0) cout<<"i Track : "<<itrack<<" nHit: "<<nhit<<" pt: "<<pt<<endl;
    }
    // vector<HoughTrack>::iterator iter_hough = vec_track.begin();
    //  vector<HoughTrack>::iterator iter_hough2 = vec_track.begin()+1;
    for(int itrack1=0; itrack1<mdcTrackList.length(); itrack1++){
        MdcTrack* track_1 = mdcTrackList[itrack1];
        TrkExchangePar par1 = track_1->track().fitResult()->helix(0.);
        int nhit1 = track_1->track().hots()->nHit();
        double d0_1 = par1.d0();
        double phi0_1 = par1.phi0();
        double omega_1 = par1.omega();
        double z0_1= par1.z0();
        double cr=fabs(1./par1.omega());
        double cx= sin(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1.; //z axis verse,x_babar = -x_belle
        double cy= -1.*cos(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1;//???
        if(m_debug>0) cout<<"circle 1 : "<<cr<<","<<cx<<","<<cy<<endl;
 
        for(int itrack2=itrack1+1; itrack2<mdcTrackList.length(); itrack2++){
            MdcTrack* track_2 = mdcTrackList[itrack2];
            TrkExchangePar par2 = track_2->track().fitResult()->helix(0.);
            int nhit2 = track_2->track().hots()->nHit();
            double d0_2 = par2.d0();
            double phi0_2 = par2.phi0();
            double omega_2 = par2.omega();
            double z0_2= par2.z0();
            double cR=fabs(1./par2.omega());
            double cX= sin(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1.; //z axis verse,x_babar = -x_belle
            double cY= -1.*cos(par2.phi0()) *(par2.d0()+1/par2.omega()) * -1;//???
            if(m_debug>0)     cout<<"circle 2 : "<<cR<<","<<cX<<","<<cY<<endl;
 
            double   bestDiff = 1.0e+20;
            if(fabs(cr)*(1.-0.25) <= fabs(cR) && fabs(cR) <= fabs(cr)*(1.+0.25)){
                if(fabs(cx-cX) <= fabs(cr)*(0.25)&& fabs(cy-cY) <= fabs(cr)*(0.25) ){
                    double diff = fabs((fabs(cr)-fabs(cR))*(cx-cX)*(cy-cY));
                    if(diff < bestDiff){
                        bestDiff = diff;
                    }
                }
            }
 
            double zdiff = z0_1-z0_2;
            if(bestDiff != 1.0e20 &&  fabs(zdiff)<25.){ 
                if(m_debug>0) cout<<"z0_1 : "<<z0_1<<" z0_2 : "<<z0_2<<endl;
                //if( fabs(z0_1)>fabs(z0_2) && fabs(z0_1)>m_z0Cut_compareHough)
                if( nhit1<nhit2 && (fabs(z0_1)>fabs(z0_2) && fabs(z0_1)>m_z0Cut_compareHough) ){   //FIXME
                    if(m_debug>0) cout<<"remove track1 "<<1./par1.omega()/333.567<<endl;
                    //remove 1 
                    mdcTrackList.remove( track_1 );
                    itrack1--;
                    break;
                }
                else{
                    //remove 2 
                    if(m_debug>0) cout<<"remove track2 "<<1./par2.omega()/333.567<<endl;
                    //remove 1 
                    mdcTrackList.remove( track_2 );
                    itrack2--;
                }
            }
            if(bestDiff == 1.0e20) {  if(m_debug>0) cout<<" no same track in hough"<<endl; }
        }
    }
}
int MdcHoughFinder::comparePatTsf(MdcTrackList& mdcTrackList, RecMdcTrackCol* trackList_tds){
    //  vector<HoughTrack>::iterator iter_hough = vec_track.begin();
    //  int itrk_hough=0;
    //  for(;iter_hough != vec_track.end();iter_hough++)
    //  if(m_debug>0) cout<<"being hough trk "<<itrk_hough<<endl;
    //  itrk_hough++;
    //  double cr= (*iter_hough).getCirR();
    //  double cx= (*iter_hough).getCirX(); 
    //  double cy= (*iter_hough).getCirY(); 
    for(int itrack1=0; itrack1<mdcTrackList.length(); itrack1++){
        MdcTrack* track_1 = mdcTrackList[itrack1];
        TrkExchangePar par1 = track_1->track().fitResult()->helix(0.);
        int nhit1 = track_1->track().hots()->nHit();
        double d0_1 = par1.d0();
        double phi0_1 = par1.phi0();
        double omega_1 = par1.omega();
        double z0_1= par1.z0();
        double cr=fabs(1./par1.omega());
        double cx= sin(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1.; //z axis verse,x_babar = -x_belle
        double cy= -1.*cos(par1.phi0()) *(par1.d0()+1/par1.omega()) * -1;//???
 
        //track2
        double   bestDiff = 1.0e+20;
        double   ptDiff = 0.0;
        double   cR, cX, cY;
        vector<double> a0,a1,a2,a3,a4;
        RecMdcTrackCol::iterator iterBest; 
        RecMdcTrackCol::iterator iteritrk = trackList_tds->begin();
        int itrk=0;
        for(;iteritrk!=trackList_tds->end();iteritrk++){
            if(m_debug>0) cout<<"being pattsf trk "<<itrk<<endl;
            double pt=(*iteritrk)->pxy();
            a0.push_back( (*iteritrk)->helix(0) );
            a1.push_back( (*iteritrk)->helix(1) );
            a2.push_back( (*iteritrk)->helix(2) );
            a3.push_back( (*iteritrk)->helix(3) );
            a4.push_back( (*iteritrk)->helix(4) );
            cR=(-333.567)/a2[itrk];
            cX=(cR+a0[itrk])*cos(a1[itrk]);                   
            cY=(cR+a0[itrk])*sin(a1[itrk]);
            if(m_debug>0) 
            {
                cout<<" compare PATTSF and HOUGH "<<endl;
                cout<<"  fabs(cX) "<<  fabs(cX)<< "fabs(cx) "<<fabs(cx)<<endl;
                cout<<"  fabs(cY) "<<  fabs(cY)<< "fabs(cy) "<<fabs(cy)<<endl;
                cout<<"  fabs(cR) "<<  fabs(cR)<< "fabs(cr) "<<fabs(cr)<<endl;
                cout<<"  fabs(cx-cX) "<<  fabs(cx-cX)<< "fabs(cy-cY) "<<fabs(cy-cY)<<" fabs(cr-cR) "<< fabs(cr-cR)<<endl;
            }
 
            if(fabs(cr)*(1.-0.25) <= fabs(cR) && fabs(cR) <= fabs(cr)*(1.+0.25)){
                if(fabs(cx-cX) <= fabs(cr)*(0.25)&& fabs(cy-cY) <= fabs(cr)*(0.25) ){
                    double diff = fabs((fabs(cr)-fabs(cR))*(cx-cX)*(cy-cY));
                    if(m_debug>0) cout<<" diff "<<diff<<" pt "<<pt<<endl;
                    if( fabs(pt) >  ptDiff){
                        ptDiff= pt;
                        iterBest=iteritrk;
                        bestDiff = diff;
                    }
                    //if(diff < bestDiff){
                    //  bestDiff = diff;
                    //  //          bestIndex = itrk;
                    //  iterBest=iteritrk;
                    //}
                }
            }
            itrk++;
        }
        if(bestDiff != 1.0e20) { 
            if(m_debug>0) cout<<" same track pattsf & hough , then compare nhit. "<<endl;
            double d0_pattsf =(*iterBest)->helix(0);
            double z0_pattsf =(*iterBest)->helix(3);
            double d0_hough = d0_1;
            double z0_hough = z0_1;
            int use_pattsf(1),use_hough(1);
            if( fabs(d0_pattsf)>m_dropTrkDrCut || fabs(z0_pattsf)>m_dropTrkDzCut ) use_pattsf=0;
            if( fabs(d0_hough)>m_dropTrkDrCut || fabs(z0_hough)>m_dropTrkDzCut ) use_hough=0;
            if( use_pattsf==0 && use_hough==1 ) {
                trackList_tds->remove(*iterBest); 
                if(m_debug>0) cout<<" use houghTrack, vertex pattsf wrong"<<endl;
            }
            if( use_pattsf==1 && use_hough==0 ) {
                mdcTrackList.remove( track_1 );
                itrack1--;
                if(m_debug>0) cout<<" use pattsfTrack, vertex hough wrong"<<endl;
            }
            if( (use_pattsf==0 && use_hough==0) || (use_pattsf==1 && use_hough==1) ){
                int nhit_pattsf=( (*iterBest)->ndof()+5 );
                int nhit_hough =  nhit1;
                if(m_debug>0) cout<<" nhit "<<nhit_pattsf<<" "<<nhit_hough<<endl;
                if( nhit_pattsf>nhit_hough ) {
                    mdcTrackList.remove( track_1 );
                    itrack1--;
                    if(m_debug>0) cout<<" use pattsfTrack "<<endl;
                }
                else{
                    trackList_tds->remove(*iterBest); 
                    if(m_debug>0) cout<<" use houghTrack "<<endl;
                }
            }
        }
        //  if(m_debug>0) cout<<" d0"<<d0_pattsf<<" "<<d0_hough<<endl;
        //  if(m_debug>0) cout<<" z0"<<z0_pattsf<<" "<<z0_hough<<endl;
        //  if( fabs(d0_pattsf) >= fabs(d0_hough) && fabs(z0_pattsf) >= fabs(z0_hough) ) {
        //    trackList_tds->remove(*iterBest); 
        //    if(m_debug>0) cout<<" use houghTrack "<<endl;
        //  }
        //  else{
        //    iter_hough = vec_track.erase(iter_hough); iter_hough--;
        //    if(m_debug>0) cout<<" use pattsfTrack "<<endl;
        //  }
 
        if(bestDiff == 1.0e20) { if(m_debug>0) cout<<" no same track in pattsf"<<endl; }
        //if(bestDiff != 1.0e20) { iter_hough = vec_track.erase(iter_hough); iter_hough--;}
        //if(bestDiff == 1.0e20) { if(m_debug>0) cout<<" no merge "<<endl; }
    }
    //cout<<"size after combine "<<trackList_tds->size()<<endl;
    if(trackList_tds->size()!=0){
        int digi=0;
        RecMdcTrackCol::iterator iter_pat = trackList_tds->begin();
        for(;iter_pat!=trackList_tds->end();iter_pat++,digi++){
            (*iter_pat)->setTrackId(digi);
            //  cout<<"digi "<<(*iter_pat)->trackId()<<endl;
        }
    }
    return trackList_tds->size();
}
//int MdcHoughFinder::judgeHit(MdcTrackList& list,const vector<HoughTrack>& trackCol)
int MdcHoughFinder::judgeHit(MdcTrackList& list,vector<MdcTrack*>& trackCol){
    if(m_debug>0) cout<<"in judgHit: "<<endl;
    for( unsigned int i =0;i<trackCol.size();i++){
        //MdcTrack *mdcTrack = new  MdcTrack(trackCol[i].getTrk());
        //list.append(mdcTrack);                                                            
        list.append(trackCol[i]);
    } 
    if(m_cgem || m_skipMDCIT)list.setNoInner(true);
    int nDeleted(0);
    //nDeleted = list.
    return nDeleted;
}
bool sortCluster(const RecCgemCluster* clusterA , const RecCgemCluster* clusterB){
    return clusterA->getlayerid()<clusterB->getlayerid();
}
int MdcHoughFinder::storeTracks(RecMdcTrackCol* trackList, RecMdcHitCol* hitList, HoughTrackList& houghTrackList){
 
    int tracklist_size = houghTrackList.getTrackNum();
    for(int itrack=0;itrack<tracklist_size;itrack++){
        HoughTrack &houghTrack = houghTrackList.getTrack(itrack);
        RecMdcTrack* recMdcTrack = new RecMdcTrack();
        recHitCol recHitVec = houghTrack.getHoughHitList();
        int    hitId = 0;
        int    nSt=0;
        int    nHits = recHitVec.size();
        int    tkStat = 4;
        int    q = houghTrack.getCharge();
        double phi0 = houghTrack.getPhi0();
        double tanDip = houghTrack.getTanl();
        double d0 = houghTrack.getD0();
        //momenta and position
        double helixPar[5];
        helixPar[0] = -d0;
        double tphi0 = phi0 - Constants::pi/2.;
        if (tphi0 < 0. ) tphi0 += Constants::pi * 2.;
        helixPar[1] = tphi0;
        helixPar[2] = houghTrack.getOmega()*333.567;
        helixPar[3] = houghTrack.getZ0();
        helixPar[4] = tanDip;
        recMdcTrack->setTrackId(itrack);
        recMdcTrack->setHelix(helixPar); 
        recMdcTrack->setCharge(q);
        recMdcTrack->setStat(tkStat);
        recMdcTrack->setNhits(nHits);
        //-----------hit list-------------
        HitRefVec  hitRefVec;
        ClusterRefVec clusterRefVec;
        for (int i=0;i<nHits;i++){
            if(recHitVec[i].getflag()!=0 ) continue;
            if(recHitVec[i].getDetectorType()==MDC){
                RecMdcHit* recMdcHit = new RecMdcHit;
                Identifier mdcId = recHitVec[i].digi()->identify();
                recMdcHit->setMdcId(mdcId);
                double fltLen = recHitVec[i].getRet().second/(1./sqrt(1+tanDip*tanDip)); 
                //fltLen = recHitVec[i].getRet().second; 
                recMdcHit->setFltLen(fltLen);
                recMdcHit->setTdc(recHitVec[i].digi()->getTimeChannel());
                recMdcHit->setAdc(recHitVec[i].digi()->getChargeChannel()); 
                if(recHitVec[i].getDetectorType()!=0)nSt++;
                hitList->push_back(recMdcHit);
                SmartRef<RecMdcHit> refHit(recMdcHit);
                hitRefVec.push_back(refHit);
            }else if(recHitVec[i].getDetectorType()==CGEM && recHitVec[i].getSlayerType()==2){
                clusterRefVec.push_back(recHitVec[i].getCluster());
            }
        }
        std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
        recMdcTrack->setNster(nSt);  
        recMdcTrack->setVecHits(hitRefVec);
        recMdcTrack->setVecClusters(clusterRefVec);
        recMdcTrack->setNcluster(clusterRefVec.size());
        trackList->push_back(recMdcTrack);
        //cout<<helixPar[2]<<endl;
    }
}//end of storeTrack
 
void MdcHoughFinder::printTrack(vector<MdcTrack*> vec_MdcTrack)
{
    //cout<<"track number: "<<vec_MdcTrack.length()<<endl;
    for( unsigned int i =0;i<vec_MdcTrack.size();i++){
        cout<<"trackID:"<<vec_MdcTrack[i]->track().id()<<"  ";
        vec_MdcTrack[i]->track().printAll(cout);
            //<<(vec_MdcTrack[i]->track().fitResult()->pt())*(vec_MdcTrack[i]->track().fitResult()->charge())<<endl;
    }
    //cout<<endl;
}
void MdcHoughFinder::printTrack(MdcTrackList& mdcTrackList)
{
    //cout<<"track number: "<<mdcTrackList.length()<<endl;
    for( unsigned int i =0;i<mdcTrackList.length();i++){
        cout<<"trackID:"<<mdcTrackList[i]->track().id()<<"  "
            <<(mdcTrackList[i]->track().fitResult()->pt())*(mdcTrackList[i]->track().fitResult()->charge())<<endl;
        mdcTrackList[i]->track().printAll(cout);
    }
    //cout<<endl;
}
 
StatusCode MdcHoughFinder::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_evt",          m_hit_evt);
            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_layer",        m_hit_nhit, m_hit_layer);
            ntuple_hit->addItem("hit_wire",         m_hit_nhit, m_hit_wire);
            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_driftdist",    m_hit_nhit, m_hit_driftdist);
            ntuple_hit->addItem("hit_drifttime",    m_hit_nhit, m_hit_drifttime);
            ntuple_hit->addItem("hit_flag",         m_hit_nhit, m_hit_flag);
            ntuple_hit->addItem("hit_truth_x",      m_hit_nhit, m_hit_truth_x);
            ntuple_hit->addItem("hit_truth_y",      m_hit_nhit, m_hit_truth_y);
            ntuple_hit->addItem("hit_truth_z",      m_hit_nhit, m_hit_truth_z);
            ntuple_hit->addItem("hit_truth_drift",  m_hit_nhit, m_hit_truth_drift);
            ntuple_hit->addItem("hit_truth_ambig",  m_hit_nhit, m_hit_truth_ambig);
        } else { log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/hit" <<endmsg;
            return StatusCode::FAILURE;
        }
    }
 
    NTuplePtr nt2(ntupleSvc(), "mdcHoughFinder/hot");
    if ( nt2 ){
        ntuple_hot = nt2;
    } else {
        ntuple_hot = ntupleSvc()->book("mdcHoughFinder/hot", CLID_ColumnWiseTuple, "hot");
        if(ntuple_hot){
            ntuple_hot->addItem("hot_run",          m_hot_run);
            ntuple_hot->addItem("hot_evt",          m_hot_evt);
            ntuple_hot->addItem("hot_trk",          m_hot_trk);
            ntuple_hot->addItem("hot_nhot",         m_hot_nhot, 0, 1000 );
            ntuple_hot->addItem("hot_hitid",        m_hot_nhot, m_hot_hitid);
            ntuple_hot->addItem("hot_layer",        m_hot_nhot, m_hot_layer);
            ntuple_hot->addItem("hot_wire",         m_hot_nhot, m_hot_wire);
            ntuple_hot->addItem("hot_x",            m_hot_nhot, m_hot_x);
            ntuple_hot->addItem("hot_y",            m_hot_nhot, m_hot_y);
            ntuple_hot->addItem("hot_z",            m_hot_nhot, m_hot_z);
            ntuple_hot->addItem("hot_x0",           m_hot_nhot, m_hot_x0);
            ntuple_hot->addItem("hot_y0",           m_hot_nhot, m_hot_y0);
            ntuple_hot->addItem("hot_z0",           m_hot_nhot, m_hot_z0);
            ntuple_hot->addItem("hot_s0",           m_hot_nhot, m_hot_s0);
            ntuple_hot->addItem("hot_x1",           m_hot_nhot, m_hot_x1);
            ntuple_hot->addItem("hot_y1",           m_hot_nhot, m_hot_y1);
            ntuple_hot->addItem("hot_z1",           m_hot_nhot, m_hot_z1);
            ntuple_hot->addItem("hot_s1",           m_hot_nhot, m_hot_s1);
            ntuple_hot->addItem("hot_drift",        m_hot_nhot, m_hot_drift);
            ntuple_hot->addItem("hot_flag",         m_hot_nhot, m_hot_flag);
            ntuple_hot->addItem("hot_deltaD",       m_hot_nhot, m_hot_deltaD);
            ntuple_hot->addItem("hot_truth_x",      m_hot_nhot, m_hot_truth_x);
            ntuple_hot->addItem("hot_truth_y",      m_hot_nhot, m_hot_truth_y);
            ntuple_hot->addItem("hot_truth_z",      m_hot_nhot, m_hot_truth_z);
            ntuple_hot->addItem("hot_truth_drift",  m_hot_nhot, m_hot_truth_drift);
            ntuple_hot->addItem("hot_truth_ambig",  m_hot_nhot, m_hot_truth_ambig);
 
        } else { log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/hot" <<endmsg;
            return StatusCode::FAILURE;
        }
    }
 
    NTuplePtr nt3(ntupleSvc(), "mdcHoughFinder/trk");
    if ( nt3 ){
        ntuple_trk = nt3;
    } else {
        ntuple_trk = ntupleSvc()->book("mdcHoughFinder/trk", CLID_ColumnWiseTuple, "trk");
        if(ntuple_trk){
            ntuple_trk->addItem("trk_run",                  m_trk_run);
            ntuple_trk->addItem("trk_evt",                  m_trk_evt);
            ntuple_trk->addItem("trk_ntrk",                 m_trk_ntrk);
            ntuple_trk->addItem("trk_size",                 m_trk_size, 0, 100);
            ntuple_trk->addItem("trk_trackId",              m_trk_size, m_trk_trackId);
            ntuple_trk->addItem("trk_charge",               m_trk_size, m_trk_charge);
            ntuple_trk->addItem("trk_dr",                   m_trk_size, m_trk_dr);
            ntuple_trk->addItem("trk_phi0",                 m_trk_size, m_trk_phi0);
            ntuple_trk->addItem("trk_kappa",                m_trk_size, m_trk_kappa);
            ntuple_trk->addItem("trk_dz",                   m_trk_size, m_trk_dz);
            ntuple_trk->addItem("trk_tanl",                 m_trk_size, m_trk_tanl);
            ntuple_trk->addItem("trk_pxy",                  m_trk_size, m_trk_pxy);
            ntuple_trk->addItem("trk_px",                   m_trk_size, m_trk_px);
            ntuple_trk->addItem("trk_py",                   m_trk_size, m_trk_py);
            ntuple_trk->addItem("trk_pz",                   m_trk_size, m_trk_pz);
            ntuple_trk->addItem("trk_p",                    m_trk_size, m_trk_p);
            ntuple_trk->addItem("trk_theta",                m_trk_size, m_trk_theta);
            ntuple_trk->addItem("trk_phi",                  m_trk_size, m_trk_phi);
            ntuple_trk->addItem("trk_x",                    m_trk_size, m_trk_x);
            ntuple_trk->addItem("trk_y",                    m_trk_size, m_trk_y);
            ntuple_trk->addItem("trk_z",                    m_trk_size, m_trk_z);
            ntuple_trk->addItem("trk_r",                    m_trk_size, m_trk_r);
            ntuple_trk->addItem("trk_chi2",                 m_trk_size, m_trk_chi2);
            ntuple_trk->addItem("trk_fiTerm",               m_trk_size, m_trk_fiTerm);
            ntuple_trk->addItem("trk_matchChi2",            m_trk_size, m_trk_matchChi2);
            ntuple_trk->addItem("trk_nhit",                 m_trk_size, m_trk_nhit);
            ntuple_trk->addItem("trk_ncluster",             m_trk_size, m_trk_ncluster);
            ntuple_trk->addItem("trk_stat",                 m_trk_size, m_trk_stat);
            ntuple_trk->addItem("trk_ndof",                 m_trk_size, m_trk_ndof);
            ntuple_trk->addItem("trk_nster",                m_trk_size, m_trk_nster);
            ntuple_trk->addItem("trk_nlayer",               m_trk_size, m_trk_nlayer);
            ntuple_trk->addItem("trk_firstLayer",           m_trk_size, m_trk_firstLayer);
            ntuple_trk->addItem("trk_lastLayer",            m_trk_size, m_trk_lastLayer);
            ntuple_trk->addItem("trk_nCgemXClusters",       m_trk_size, m_trk_nCgemXClusters);
            ntuple_trk->addItem("trk_nCgemVClusters",       m_trk_size, m_trk_nCgemVClusters);
            ntuple_trk->addItem("trk_nhop",                 m_trk_size, m_trk_nhop);
            ntuple_trk->addItem("trk_nhot",                 m_trk_size, m_trk_nhot);
            ntuple_trk->addItem("trk_Xc",                   m_trk_size, m_trk_Xc);
            ntuple_trk->addItem("trk_Yc",                   m_trk_size, m_trk_Yc);
            ntuple_trk->addItem("trk_R",                    m_trk_size, m_trk_R);
 
            ntuple_trk->addItem("trk_truth_charge",         m_trk_size, m_trk_truth_charge);
            ntuple_trk->addItem("trk_truth_dr",             m_trk_size, m_trk_truth_dr);
            ntuple_trk->addItem("trk_truth_phi0",           m_trk_size, m_trk_truth_phi0);
            ntuple_trk->addItem("trk_truth_kappa",          m_trk_size, m_trk_truth_kappa);
            ntuple_trk->addItem("trk_truth_dz",             m_trk_size, m_trk_truth_dz);
            ntuple_trk->addItem("trk_truth_tanl",           m_trk_size, m_trk_truth_tanl);
            ntuple_trk->addItem("trk_truth_pxy",            m_trk_size, m_trk_truth_pxy);
            ntuple_trk->addItem("trk_truth_px",             m_trk_size, m_trk_truth_px);
            ntuple_trk->addItem("trk_truth_py",             m_trk_size, m_trk_truth_py);
            ntuple_trk->addItem("trk_truth_pz",             m_trk_size, m_trk_truth_pz);
            ntuple_trk->addItem("trk_truth_p",              m_trk_size, m_trk_truth_p);
            ntuple_trk->addItem("trk_truth_theta",          m_trk_size, m_trk_truth_theta);
            ntuple_trk->addItem("trk_truth_phi",            m_trk_size, m_trk_truth_phi);
            ntuple_trk->addItem("trk_truth_x",              m_trk_size, m_trk_truth_x);
            ntuple_trk->addItem("trk_truth_y",              m_trk_size, m_trk_truth_y);
            ntuple_trk->addItem("trk_truth_z",              m_trk_size, m_trk_truth_z);
            ntuple_trk->addItem("trk_truth_r",              m_trk_size, m_trk_truth_r);
            ntuple_trk->addItem("trk_truth_cosTheta",       m_trk_size, m_trk_truth_cosTheta);
            ntuple_trk->addItem("trk_truth_Xc",             m_trk_size, m_trk_truth_Xc);
            ntuple_trk->addItem("trk_truth_Yc",             m_trk_size, m_trk_truth_Yc);
            ntuple_trk->addItem("trk_truth_R",              m_trk_size, m_trk_truth_R);
 
        } else { log << MSG::ERROR << "Cannot book tuple mdcHoughFinder/trk" <<endmsg;
            return StatusCode::FAILURE;
        }
    }
    return StatusCode::SUCCESS;
}