MdcHoughFinder Class Reference

#include <Hough.h>

Inheritance diagram for MdcHoughFinder:

Public Member Functions
	MdcHoughFinder (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
StatusCode	beginRun ()
StatusCode	bookTuple ()
	MdcHoughFinder (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
StatusCode	beginRun ()
StatusCode	bookTuple ()

Detailed Description

Definition at line 49 of file InstallArea/include/MdcHoughFinder/MdcHoughFinder/Hough.h.

Constructor & Destructor Documentation

MdcHoughFinder() [1/2]

MdcHoughFinder::MdcHoughFinder	(	const std::string &	name,
		ISvcLocator *	pSvcLocator
	)

Definition at line 76 of file MdcHoughFinder/MdcHoughFinder-00-00-12/src/Hough.cxx.

                                                                              :
    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
 
}

MdcHoughFinder() [2/2]

MdcHoughFinder::MdcHoughFinder	(	const std::string &	name,
		ISvcLocator *	pSvcLocator
	)

Member Function Documentation

beginRun() [1/2]

StatusCode MdcHoughFinder::beginRun

(

)

Definition at line 149 of file MdcHoughFinder/MdcHoughFinder-00-00-12/src/Hough.cxx.

                                   {
    //Initailize MdcDetector
    Global::m_gm = MdcDetector::instance(0);
    if(NULL == Global::m_gm) return StatusCode::FAILURE;
 
    return StatusCode::SUCCESS;
}

beginRun() [2/2]

StatusCode MdcHoughFinder::beginRun

(

)

bookTuple() [1/2]

StatusCode MdcHoughFinder::bookTuple

(

)

Definition at line 1727 of file MdcHoughFinder/MdcHoughFinder-00-00-12/src/Hough.cxx.

                                    {
    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;
}

Referenced by initialize().

bookTuple() [2/2]

StatusCode MdcHoughFinder::bookTuple

(

)

execute() [1/2]

StatusCode MdcHoughFinder::execute

(

)

Definition at line 307 of file MdcHoughFinder/MdcHoughFinder-00-00-12/src/Hough.cxx.

                                   {
 
    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;
}

execute() [2/2]

StatusCode MdcHoughFinder::execute

(

)

finalize() [1/2]

StatusCode MdcHoughFinder::finalize

(

)

Definition at line 709 of file MdcHoughFinder/MdcHoughFinder-00-00-12/src/Hough.cxx.

                                    {
    MsgStream log(msgSvc(), name());
    delete m_bfield ;
    log << MSG::INFO<< "in finalize()" << endreq;
 
    return StatusCode::SUCCESS;
}

finalize() [2/2]

StatusCode MdcHoughFinder::finalize

(

)

initialize() [1/2]

StatusCode MdcHoughFinder::initialize

(

)

Definition at line 158 of file MdcHoughFinder/MdcHoughFinder-00-00-12/src/Hough.cxx.

                                     {
 
    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;
}

initialize() [2/2]

StatusCode MdcHoughFinder::initialize

(

)

---

The documentation for this class was generated from the following files:

• InstallArea/include/MdcHoughFinder/MdcHoughFinder/Hough.h
• Reconstruction/MdcHoughFinder/MdcHoughFinder-00-00-12/MdcHoughFinder/Hough.h
• MdcHoughFinder/MdcHoughFinder-00-00-12/src/Hough.cxx