CgemMdcFitAlg Class Reference

#include <CgemMdcFitAlg.h>

Inheritance diagram for CgemMdcFitAlg:

Public Types
enum	TrackType { lineType = 2 , circleType = 3 , helixType = 5 }

Public Member Functions
	CgemMdcFitAlg (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
void	updateTracks (int trackId, TrkRecoTrk *trkRecoTrk, RecMdcTrack *recMdcTrack, int trkStat)
void	compareTracks (MdcTrackList &mdcTrackList, double dzCut)
TrkErrCode	fit (RecMdcTrack *recMdcTrack, TrkRecoTrk *trkRecoTrk, TrackType trackType)
TrkErrCode	check (TrkRecoTrk *trkRecoTrk, TrackType trackType)

Detailed Description

Definition at line 14 of file CgemMdcFitAlg.h.

Member Enumeration Documentation

TrackType

enum CgemMdcFitAlg::TrackType

Enumerator
lineType
circleType
helixType

Definition at line 17 of file CgemMdcFitAlg.h.

{ lineType = 2, circleType = 3, helixType = 5};

Constructor & Destructor Documentation

CgemMdcFitAlg()

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

Definition at line 60 of file CgemMdcFitAlg.cxx.

                                                                           :
    Algorithm(name, pSvcLocator)
{
    declareProperty("debug"         ,m_debug = 0            );
    declareProperty("tuple"         ,m_tuple = 0            );
    declareProperty("change"        ,m_changeTDS = 2        );// 0:not change, 1:update, 2:replace
    declareProperty("debugArbHit"   ,m_debugArbHit= 0       );
    declareProperty("fit2D"         ,m_fit2D= 1             );
    declareProperty("drCut"         ,m_drCut= 1.            );
    declareProperty("dzCut"         ,m_dzCut= 10.           );
    declareProperty("filter"        ,m_filter= 0            );
    declareProperty("eventFile"     ,m_evtFile= "EventList" );
    declareProperty("compareTrack"  ,m_compareTrack= 0      );
    declareProperty("removeBadTrack",m_removeBadTrack= 0    );
}

Member Function Documentation

check()

TrkErrCode CgemMdcFitAlg::check	(	TrkRecoTrk *	trkRecoTrk,
		TrackType	trackType
	)

Definition at line 563 of file CgemMdcFitAlg.cxx.

{
    m_qualityFactor=3;
    m_dropTrkDrCut=1;
    m_dropTrkDzCut=10;
    if(trackType==circleType)m_dropTrkNhitCut=5;
    if(trackType==helixType)m_dropTrkNhitCut=9;
    m_dropTrkChi2Cut=99999;
    m_dropTrkChi2NdfCut=99999;
 
    TrkErrCode err;
    int nActiveHit = trkRecoTrk->hots()->nActive();
    TrkHitList* qhits = trkRecoTrk->hits();
    const TrkFit* trkFit = trkRecoTrk->fitResult();
    //if (!trkFit){cout<<"fit failed"<<endl;}
    //else{cout<<"fit success"<<endl;}
    //cout<<trkFit<<endl;
    //TrkExchangePar par = trkFit->helix(0.);
    int fit_stat=false;
    double chi2=-999.;
    if (!trkFit){    
        err.setFailure(12,"fit failed");
        return err;
    }else{
        TrkExchangePar par = trkFit->helix(0.);
        if( abs( 1/(par.omega()) )>0.03) {
            fit_stat = 1;
            chi2=trkFit->chisq();
        }
        else {
            fit_stat = 0;
            chi2=-999;
            if(m_debug>0){
                std::cout<<__FILE__<<"   "<<__LINE__<<" drop track radius cut  "<<abs(1/(par.omega()))<<" < "<<0.03 <<std::endl;
            }
        }
 
        bool badQuality = false;
        if(fabs(chi2)>m_qualityFactor*m_dropTrkChi2Cut ){
            if(m_debug>0){
                std::cout<<__FILE__<<"   "<<__LINE__<<" drop track by chi2  "<<chi2<<" > "<<m_qualityFactor<<" * "<<m_dropTrkChi2Cut <<std::endl;
            }
            badQuality=1;
        }
        if( fabs(par.d0())>m_qualityFactor*m_dropTrkDrCut) {
            if(m_debug>0){
                std::cout<<__FILE__<<"   "<<__LINE__<<" drop track by d0 "<<par.d0()<<" > "<<m_qualityFactor<<" * "<<m_dropTrkDrCut <<std::endl;
            }
            badQuality=1;
        }
        if( fabs(par.z0())>m_qualityFactor*m_dropTrkDzCut && trackType == helixType) {
            if(m_debug>0){
                std::cout<<__FILE__<<"   "<<__LINE__<<" drop track by z0 "<<par.z0()<<" > "<<m_qualityFactor<<" * "<<m_dropTrkDzCut <<std::endl;
            }
            badQuality=1;
        }
        if( (fabs(chi2)/qhits->nHit()) > m_qualityFactor*m_dropTrkChi2NdfCut) {
            if(m_debug>0){
                std::cout<<__FILE__<<"   "<<__LINE__<<" drop track by chi2/ndf "<<(fabs(chi2)/qhits->nHit()) <<" > "<<m_qualityFactor<<" * "<<m_dropTrkChi2NdfCut<<std::endl;
            }
            badQuality=1;
        }
        if( nActiveHit <= m_dropTrkNhitCut) {
            if(m_debug>0){
                std::cout<<__FILE__<<"   "<<__LINE__<<" drop track by nhit "<<nActiveHit <<" <= "<<m_dropTrkNhitCut<<std::endl;
            }
            badQuality=1;
        }
        //badQuality = false;   //not delete track in this stage
        if( badQuality) fit_stat=0;
    }
    //if( fit_stat!=1 ){
        //delete trkRecoTrk;
        //vectrk_for_clean.pop_back();
    //}
    if( fit_stat==1 ){
 
        if(m_debug>1){
            TrkExchangePar par = trkFit->helix(0.);
            double _d0=par.d0();
            double _phi0=par.phi0();
            double _omega=par.omega();
            double _tanl=0;
            double _z0=0;
            if(trackType == helixType){
                _tanl=par.tanDip();
                _z0=par.z0();
            }
            //double _tanl=par.tanDip();
            //double _z0=par.z0();
            int _charge=trkFit->charge();
            cout << " global 3d fit success"<<endl;
            //cout<<__FILE__<<__LINE__<<"AFTER fit 3d "<<endl;
            trkRecoTrk->print(std::cout);
            double _circleR=_charge/par.omega();
            double _circleX= sin(par.phi0()) *(par.d0()+1/par.omega()) * -1.; //z axis verse,x_babar = -x_belle
            double _circleY= -1.*cos(par.phi0()) *(par.d0()+1/par.omega()) * -1;//???
            double pxy=fabs(_circleR/333.567);
            double pz=pxy*_tanl;
            double pxyz=_charge*sqrt(pxy*pxy+pz*pz);
            double _pt=pxy; 
            double _pz=pz;
            double _p=pxyz;
            cout<<" circle after globle 3d:  "<<"("<<_circleX<<","<<_circleY<<","<<_circleR<<")"<<endl;
            cout<<"pt_rec:   "<<_pt<<endl;
            cout<<"pz_rec:   "<<_pz<<endl;
            cout<<"p_rec:   "<<_p<<endl;
 
            int nfit3d=0;
            cout<<" hot list:"<<endl;
            TrkHotList::hot_iterator hotIter= qhits->hotList().begin();
            int lay=((MdcHit*)(hotIter->hit()))->layernumber();
            int wir=((MdcHit*)(hotIter->hit()))->wirenumber();
            while (hotIter!=qhits->hotList().end()) {
                if( typeid(*hotIter)==typeid(MdcRecoHitOnTrack) || typeid(*hotIter)==typeid(MdcHitOnTrack)){
                    cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber()
                        <<","<<((MdcHit*)(hotIter->hit()))->wirenumber() 
                        <<":"<<hotIter->isActive()<<")  ";
                }else if(typeid(*hotIter)==typeid(CgemHitOnTrack)){
                    const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hotIter));
                    const RecCgemCluster* recCgemCluster = recoHot->baseHit();
                    cout<<"(" <<recCgemCluster->getlayerid()
                        <<","<<recCgemCluster->getsheetid()
                        <<":"<<hotIter->isActive()<<")  ";
                    //<<"clusterID:"<<recCgemCluster->getclusterid()
                    //<<"  flag:"<<recCgemCluster->getflag()<<endl;
                    //hotIter++;
                    //continue;
                }
                //    cout << "nuse "<<hotIter->hit()->nUsedHits()<<endl;
                if (hotIter->isActive()==1) nfit3d++;
                hotIter++;
            }
            double _nfit=nfit3d;
            cout<<" in 3D fit number of Active hits : "<<_nfit<<endl;
        }
        err.setSuccess(11,"good track");
        return err;
    }else {
        err.setFailure(12,"drop by cut");
        return err;
    }
    //_chi2_aver=chi2/_nfit;
    //  for(int i=0;i<t_hitCol.size();i++){
    //  delete t_hitCol[i];
    //  }
    //return fit_stat;
 
}

Referenced by execute().

compareTracks()

void CgemMdcFitAlg::compareTracks	(	MdcTrackList &	mdcTrackList,
		double	dzCut
	)

Definition at line 1185 of file CgemMdcFitAlg.cxx.

                                                                         {
    if(m_debug>1) 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>1) cout<<"i Track : "<<itrack<<" nHit: "<<nhit<<" pt: "<<pt<<endl;
    }
 
    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>1) 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>1)     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>1) cout<<"z0_1 : "<<z0_1<<" z0_2 : "<<z0_2<<endl;
                if( nhit1<nhit2 && (fabs(z0_1)>fabs(z0_2) && fabs(z0_1)>dzCut) ){   //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>1) cout<<" no same track in hough"<<endl; }
        }
    }
}

Referenced by execute().

execute()

StatusCode CgemMdcFitAlg::execute

(

)

Definition at line 123 of file CgemMdcFitAlg.cxx.

{
//cout<<"++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"<<endl;
 
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "In CgemMdcFitAlg execute()" << endreq;
 
    SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
    if (!eventHeader) {
        log << MSG::FATAL << "Could not find Event Header" << endreq;
        return StatusCode::FAILURE;
    }
    int run = eventHeader->runNumber();
    int evt = eventHeader->eventNumber();
    SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),  "/Event/Recon/RecMdcTrackCol");
    if (!recMdcTrackCol){
        log << MSG::WARNING << "Could not find RecMdcTrackCol" << endreq;
        return StatusCode::SUCCESS;
    }
 
    SmartDataPtr<RecEsTimeCol> evTimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
    if (!evTimeCol) {
        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
        }
    }
 
    SmartDataPtr<MdcDigiCol> mdcDigiCol(evtSvc(),"/Event/Digi/MdcDigiCol");
    if (!mdcDigiCol) {
        log << MSG::FATAL << "Could not find MdcDigiCol!" << endreq;
        return StatusCode::SUCCESS;
    }
 
    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(),evt);
        if( iter_evt == evtlist.end() )  { setFilterPassed(false); return StatusCode::SUCCESS; }
        else setFilterPassed(true);
    }
 
    if(m_debug)cout<<endl;
    if(m_debug)cout<<"run:"<<run<<"    event:"<<evt<<endl;
    if(m_tuple)getMcTruth();
    if(m_debug)cout<<"tracks before fitting: "<<recMdcTrackCol->size()<<endl;
    //if(recMdcTrackCol->size()==0)return StatusCode::SUCCESS;
    if(m_tuple){
        m_run     = run;
        m_event   = evt;
        m_nTrkRec = recMdcTrackCol->size();
    }
    
    for(RecMdcTrackCol::iterator itTrk = recMdcTrackCol->begin();itTrk != recMdcTrackCol->end();itTrk++)
    {
        HitRefVec vechits = (*itTrk)->getVecHits();
        HitRefVec::iterator itHit = vechits.begin();
        for( ; itHit != vechits.end(); itHit++)
        {
            double tdc = (*itHit)->getTdc();
            double adc = (*itHit)->getAdc();
            const Identifier mdcid = (*itHit)->getMdcId();
            int layer = MdcID::layer(mdcid);
            int wire  = MdcID::wire(mdcid);
            //cout<<"layer: "<<layer<<"    wire: "<<wire<<"  tdc: "<<tdc<<"  adc: "<<adc<<endl;
            MdcDigiCol::iterator iterDigi = mdcDigiCol->begin();
            for(; iterDigi != mdcDigiCol->end(); iterDigi++ ) {
                if((*iterDigi)->identify() == mdcid) break;
            }
            if(tdc==0)tdc = (*iterDigi)->getTimeChannel();
            //cout<<"  tdc= "<<(*iterDigi)->getTimeChannel()<<"  adc= "<<(*iterDigi)->getChargeChannel()<<endl;
 
            bool used = false;
            vector<MdcHit*>::iterator imdcHit = mdcHitCol.begin(); 
            for(;imdcHit !=  mdcHitCol.end();imdcHit++){
                if((*imdcHit)->mdcId() == mdcid){
                    used = true;
                    //cout<<"layer: "<<layer<<"    view: "<<mdcHit->layer()->view()<<"    wire: "<<wire<<endl;
                    break;
                }
            }
            if(used)continue;
            //cout<<"layer: "<<layer<<"    wire: "<<wire<<"  tdc: "<<tdc<<"  adc: "<<adc<<endl;
            const MdcDigi *mdcDigi = new MdcDigi(mdcid,tdc,adc);
            MdcHit *mdcHit = new MdcHit(mdcDigi,m_mdcDetector);
            mdcHitCol.push_back(mdcHit);
        }
    }
    //cout<<mdcHitCol.size()<<endl;
    MdcTrackParams  m_par;
    if(m_debugArbHit>0 ) m_par.lPrint=8;
    m_par.lRemoveInActive=1;
    //m_par.lUseQualCuts=0;
    //m_par.maxGapLength=99;
    //m_par.maxChisq=99;
    //m_par.minHits=99;
    MdcTrackList  mdcTrackList(m_par);
 
    RecMdcTrackCol* trackList_tds;
    RecMdcHitCol* hitList_tds;
    //vector<MdcTrack*> vec_MdcTrack;
    RecMdcTrackCol::iterator itTrk = recMdcTrackCol->begin();
    for(int itrk=0;itTrk != recMdcTrackCol->end(); itrk++,itTrk++)
    {
/*
        int     trkId = (*itTrk)->trackId();
        double  dr    = (*itTrk)->helix(0);
        double  phi0  = (*itTrk)->helix(1);
        double  kappa = (*itTrk)->helix(2);
        double  dz    = (*itTrk)->helix(3);
        double  tanl  = (*itTrk)->helix(4);
        int     chargeRec   = (*itTrk)->charge()        ;
        int     statRec     = (*itTrk)->stat()          ;
        int     nhitsRec    = (*itTrk)->getNhits()      ;
        int     nclusterRec = (*itTrk)->getNcluster()   ;
        int     nsterRec    = (*itTrk)->nster()         ;
        int     ndofRec     = (*itTrk)->ndof()          ;
        double  chi2Rec     = (*itTrk)->chi2()          ;
        double  pxRec       = (*itTrk)->px()            ;
        double  pyRec       = (*itTrk)->py()            ;
        double  pzRec       = (*itTrk)->pz()            ;
        double  pxyRec      = (*itTrk)->pxy()           ;
        double  pRec        = (*itTrk)->p()             ;
        double  xRec        = (*itTrk)->x()             ;
        double  yRec        = (*itTrk)->y()             ;
        double  zRec        = (*itTrk)->z()             ;
        double  rRec        = (*itTrk)->r()             ;
        double  thetaRec    = (*itTrk)->theta()         ;
        double  phiRec      = (*itTrk)->phi()           ;
        double  fiTermRec   = (*itTrk)->getFiTerm()     ;
        if(m_tuple){
            m_trkidRec[itrk] = trkId;
            m_drRec[itrk]    = dr   ;
            m_phi0Rec[itrk]  = phi0 ;
            m_kappaRec[itrk] = kappa;
            m_dzRec[itrk]    = dz   ;
            m_tanlRec[itrk]  = tanl ;
            m_chargeRec[itrk]  = chargeRec  ;
            m_statRec[itrk]    = statRec    ;
            m_nhitsRec[itrk]   = nhitsRec   ;
            m_nclusterRec[itrk]= nclusterRec;
            m_nsterRec[itrk]   = nsterRec   ;
            m_ndofRec[itrk]    = ndofRec    ;
            m_chi2Rec[itrk]    = chi2Rec    ;
            m_pxRec[itrk]      = pxRec      ;
            m_pyRec[itrk]      = pyRec      ;
            m_pzRec[itrk]      = pzRec      ;
            m_pxyRec[itrk]     = pxyRec     ;
            m_pRec[itrk]       = pRec       ;
            m_xRec[itrk]       = xRec       ;
            m_yRec[itrk]       = yRec       ;
            m_zRec[itrk]       = zRec       ;
            m_rRec[itrk]       = rRec       ;
            m_thetaRec[itrk]   = thetaRec   ;
            m_phiRec[itrk]     = phiRec     ;
            m_fiTermRec[itrk]  = fiTermRec  ;
        }
*/
        if(m_tuple){
            m_trkidRec[itrk] = (*itTrk)->trackId();
            m_drRec[itrk]    = (*itTrk)->helix(0);
            m_phi0Rec[itrk]  = (*itTrk)->helix(1);
            m_kappaRec[itrk] = (*itTrk)->helix(2);
            m_dzRec[itrk]    = (*itTrk)->helix(3);
            m_tanlRec[itrk]  = (*itTrk)->helix(4);
            m_chargeRec[itrk]  = (*itTrk)->charge()        ;
            m_statRec[itrk]    = (*itTrk)->stat()          ;
            m_nhitsRec[itrk]   = (*itTrk)->getNhits()      ;
            m_nclusterRec[itrk]= (*itTrk)->getNcluster()   ;
            m_nsterRec[itrk]   = (*itTrk)->nster()         ;
            m_ndofRec[itrk]    = (*itTrk)->ndof()          ;
            m_chi2Rec[itrk]    = (*itTrk)->chi2()          ;
            m_pxRec[itrk]      = (*itTrk)->px()            ;
            m_pyRec[itrk]      = (*itTrk)->py()            ;
            m_pzRec[itrk]      = (*itTrk)->pz()            ;
            m_pxyRec[itrk]     = (*itTrk)->pxy()           ;
            m_pRec[itrk]       = (*itTrk)->p()             ;
            m_xRec[itrk]       = (*itTrk)->x()             ;
            m_yRec[itrk]       = (*itTrk)->y()             ;
            m_zRec[itrk]       = (*itTrk)->z()             ;
            m_rRec[itrk]       = (*itTrk)->r()             ;
            m_thetaRec[itrk]   = (*itTrk)->theta()         ;
            m_phiRec[itrk]     = (*itTrk)->phi()           ;
            m_fiTermRec[itrk]  = (*itTrk)->getFiTerm()     ;
        }
 
        if(m_debug)cout<<"fitting track: "<<itrk<<endl;
 
        //---make tracks
        float chisum =0.;
        bool permitFlips = true;
        bool lPickHits = true;
        int     trkId = (*itTrk)->trackId();
        double  dr    = (*itTrk)->helix(0);
        double  phi0  = (*itTrk)->helix(1);
        double  kappa = (*itTrk)->helix(2);
        double  dz    = (*itTrk)->helix(3);
        double  tanl  = (*itTrk)->helix(4);
        if(m_debug)cout<<"helix parameters before fitting: "<<setw(12)<<dr<<"  "<<setw(12)<<phi0<<"  "<<setw(12)<<kappa<<"  "<<setw(12)<<dz<<"  "<<setw(12)<<tanl<<endl;
 
        //TrkLineMaker   linefactory;
        //---fitting circle
        if(m_fit2D){
            TrkExchangePar circle(-dr,phi0+M_PI/2,kappa/333.567,0,0);
 
            trackType = circleType;
            TrkCircleMaker circlefactory;
            TrkRecoTrk* circleTrk = circlefactory.makeTrack(circle, chisum, *m_context, m_bunchT0);
            circlefactory.setFlipAndDrop(*circleTrk, permitFlips, lPickHits);
            TrkErrCode fit_circle = fit((*itTrk), circleTrk, trackType);
            TrkErrCode check_circle = check(circleTrk,trackType);
            if (m_debug && fit_circle.failure()){cout<<"fit: "<<fit_circle<<endl;}//FIXME
            if (m_debug && check_circle.failure()){cout<<"check: "<<check_circle<<endl;}//FIXME
            if (fit_circle.failure() || check_circle.failure()){continue;}//FIXME
 
            //---update initial fitting parameters of helix dr, phi0, kappa
            TrkExchangePar par = circleTrk->fitResult()->helix(0.);
            dr    = -par.d0();
            phi0  = par.phi0()-M_PI/2;
            if(phi0<0)phi0+=2*M_PI;
            kappa = par.omega()*333.567;
            if(m_debug)cout<<"helix parameters after 2D fitting: "<<setw(12)<<dr<<"  "<<setw(12)<<phi0<<"  "<<setw(12)<<kappa<<"  "<<setw(12)<<dz<<"  "<<setw(12)<<tanl<<endl;
        }
 
        //--- fitting helix
        TrkExchangePar helix(-dr,phi0+M_PI/2,kappa/333.567,dz,tanl);
 
        trackType = helixType;
        TrkHelixMaker  helixfactory;
        TrkRecoTrk* helixTrk = helixfactory.makeTrack(helix, chisum, *m_context, m_bunchT0);
        //const TrkId id = helixTrk->id();
        //cout<<"======================================================"<<long(id)<<endl;
        helixfactory.setFlipAndDrop(*helixTrk, permitFlips, lPickHits);
        TrkErrCode fit_helix = fit((*itTrk), helixTrk, trackType);
        TrkErrCode check_helix = check(helixTrk,trackType);
        if (m_debug && fit_helix.failure()){cout<<"fit: "<<fit_helix<<endl;}//FIXME
        if (m_debug && check_helix.failure()){cout<<"check: "<<check_helix<<endl;}//FIXME
        if (fit_helix.failure() || check_helix.failure()){continue;}//FIXME
 
        TrkExchangePar par = helixTrk->fitResult()->helix(0.);
        dr    = -par.d0();
        phi0  = par.phi0()-M_PI/2;
        if(phi0<0)phi0+=2*M_PI;
        kappa = par.omega()*333.567;
        dz=par.z0();
        tanl=par.tanDip();
        if(m_debug)cout<<"helix parameters after  fitting: "<<setw(12)<<dr<<"  "<<setw(12)<<phi0<<"  "<<setw(12)<<kappa<<"  "<<setw(12)<<dz<<"  "<<setw(12)<<tanl<<endl<<endl;
 
        //RecMdcTrack* recMdcTrack = (*itTrk);
        int trkStat = 5;//track fi by CgemMdcFitAlg set stat=5
        if(m_changeTDS ==1)updateTracks(trkId,helixTrk,(*itTrk),trkStat);
        //mdcTrack->storeTrack(trkId, trackList_tds, hitList_tds, trkStat);
        //vec_MdcTrack.push_back(mdcTrack);
        //trackList_tds->push_back(*itTrk);
        //MdcTrack *mdcTrack ;
        if(m_changeTDS ==2){
            MdcTrack *mdcTrack = new MdcTrack(helixTrk);
            //mdcTrack = new MdcTrack(helixTrk);
            mdcTrackList.append(mdcTrack);
        }
 
    }
 
    if(m_changeTDS ==2){
        mdcTrackList.setNoInner(true);
        if(m_removeBadTrack){
            int nDeleted=mdcTrackList.arbitrateHits();//FIXME
            if(m_debug>0 && nDeleted>0)cout<<nDeleted<<"  tracks have been deleted by MdcTrackList..arbitrateHits()"<<endl;//FIXME
        }
        if(m_compareTrack>0)compareTracks(mdcTrackList, m_dzCut);
        storeTracks(trackList_tds, hitList_tds,mdcTrackList);
    }
 
    if(m_tuple){
        SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(),  "/Event/Recon/RecMdcTrackCol");
        if (!recMdcTrackCol){
            log << MSG::WARNING << "Could not find RecMdcTrackCol" << endreq;
            return StatusCode::SUCCESS;
        }
        m_nTrkFit = recMdcTrackCol->size();
        RecMdcTrackCol::iterator itTrk = recMdcTrackCol->begin();
        for(int itrk=0;itTrk != recMdcTrackCol->end(); itrk++,itTrk++)
        {
            //RecMdcTrack* recMdcTrack = (*itTrk);
            m_trkidFit[itrk] = (*itTrk)->trackId();
            m_drFit[itrk]    = (*itTrk)->helix(0);
            m_phi0Fit[itrk]  = (*itTrk)->helix(1);
            m_kappaFit[itrk] = (*itTrk)->helix(2);
            m_dzFit[itrk]    = (*itTrk)->helix(3);
            m_tanlFit[itrk]  = (*itTrk)->helix(4);
            m_chargeFit[itrk]  = (*itTrk)->charge()      ;
            m_statFit[itrk]    = (*itTrk)->stat()        ;
            m_nhitsFit[itrk]   = (*itTrk)->getNhits()    ;
            m_nclusterFit[itrk]= (*itTrk)->getNcluster() ;
            m_nsterFit[itrk]   = (*itTrk)->nster()       ;
            m_ndofFit[itrk]    = (*itTrk)->ndof()        ;
            m_chi2Fit[itrk]    = (*itTrk)->chi2()        ;
            m_pxFit[itrk]      = (*itTrk)->px()          ;
            m_pyFit[itrk]      = (*itTrk)->py()          ;
            m_pzFit[itrk]      = (*itTrk)->pz()          ;
            m_pxyFit[itrk]     = (*itTrk)->pxy()         ;
            m_pFit[itrk]       = (*itTrk)->p()           ;
            m_xFit[itrk]       = (*itTrk)->x()           ;
            m_yFit[itrk]       = (*itTrk)->y()           ;
            m_zFit[itrk]       = (*itTrk)->z()           ;
            m_rFit[itrk]       = (*itTrk)->r()           ;
            m_thetaFit[itrk]   = (*itTrk)->theta()       ;
            m_phiFit[itrk]     = (*itTrk)->phi()         ;
            m_fiTermFit[itrk]  = (*itTrk)->getFiTerm()   ;
        }
    }
 
//cout<<"--------------------------------------------------------------------------------"<<endl;
    if(m_tuple)ntuple_evt->write();   
    //if(m_debug)cout<<endl;
    vector<MdcHit*>::iterator imdcHit = mdcHitCol.begin();
    for(;imdcHit !=  mdcHitCol.end();imdcHit++){
        delete (*imdcHit);
    }
    mdcHitCol.clear();
    return StatusCode::SUCCESS;
}

finalize()

StatusCode CgemMdcFitAlg::finalize

(

)

Definition at line 455 of file CgemMdcFitAlg.cxx.

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

fit()

TrkErrCode CgemMdcFitAlg::fit	(	RecMdcTrack *	recMdcTrack,
		TrkRecoTrk *	trkRecoTrk,
		TrackType	trackType
	)

Definition at line 463 of file CgemMdcFitAlg.cxx.

{
    //---make tracks
    //double  dr    = recMdcTrack->helix(0);
    //double  phi0  = recMdcTrack->helix(1);
    //double  kappa = recMdcTrack->helix(2);
    //double  dz    = recMdcTrack->helix(3);
    //double  tanl  = recMdcTrack->helix(4);
    //TrkHelixMaker trackfactory;
    //TrkLineMaker   trackfactory;
    //TrkCircleMaker trackfactory;
    //if(typeid(trackfactory)==typeid(TrkHelixMaker)){
        //TrkExchangePar tt(-dr,phi0+M_PI/2,kappa/333.567,dz,tanl);
    //}else if(typeid(trackfactory)==typeid(TrkCircleMaker)){
        //TrkExchangePar tt(-dr,phi0+M_PI/2,kappa/333.567,0,0);
    //}else if(typeid(trackfactory)==typeid(TrkLineMaker)){
        //FIXME
    //}else return TrkErrCode(TrkErrCode::fail);
    //float chisum =0.;
    //bool permitFlips = true;
    //bool lPickHits = true;
    //TrkExchangePar tt(-dr,phi0+M_PI/2,kappa/333.567,dz,tanl);
    //trkRecoTrk = trackfactory.makeTrack(tt, chisum, *m_context, m_bunchT0);
    //trackfactory.setFlipAndDrop(*trkRecoTrk, permitFlips, lPickHits);
    TrkHitList* trkHitList = trkRecoTrk->hits();
 
    //--- for ODC hits
    double ddCut=1;
    HitRefVec vechits = recMdcTrack->getVecHits();
    HitRefVec::iterator itHit = vechits.begin();
    for( ; itHit != vechits.end(); itHit++)
    {
        double flight = (*itHit)->getFltLen();
        double tdc = (*itHit)->getTdc();
        double adc = (*itHit)->getAdc();
        const Identifier mdcid = (*itHit)->getMdcId();
        int layer = MdcID::layer(mdcid);
        int wire  = MdcID::wire(mdcid);
        MdcHit* mdcHit;
        if(trackType == circleType){
            vector<MdcHit*>::iterator imdcHit = mdcHitCol.begin(); 
            for(;imdcHit !=  mdcHitCol.end();imdcHit++){
                if((*imdcHit)->mdcId() == mdcid){
                    //mdcHit = *imdcHit;
                    tdc = (*imdcHit)->tdcIndex();
                    break;
                }
            }
            const MdcDigi *mdcDigi = new MdcDigi(mdcid,tdc,adc);
            mdcHit = new MdcHit(mdcDigi,m_mdcDetector);
        }else if(trackType == helixType){
            vector<MdcHit*>::iterator imdcHit = mdcHitCol.begin(); 
            for(;imdcHit !=  mdcHitCol.end();imdcHit++){
                if((*imdcHit)->mdcId() == mdcid){
                    mdcHit = *imdcHit;
                    break;
                }
            }
        }else{
           TrkErrCode err; 
            err.setFailure(10,"fit failed,track type error!");
            return err;
        }
        //cout<<"layer: "<<layer<<"    wire: "<<wire<<"  tdc= "<<tdc<<"  adc: "<<adc<<endl;
 
        mdcHit->setCalibSvc(m_mdcCalibFunSvc);
        mdcHit->setCountPropTime(true);
 
        int newAmbig = 0;
        MdcRecoHitOnTrack temp( *mdcHit, newAmbig, m_bunchT0*1.e9);
        MdcHitOnTrack* newhot = &temp;
        newhot->setFltLen( flight );
        if(m_debug>1)cout<<"layer: "<<layer<<"    view: "<<mdcHit->layer()->view()<<"    wire: "<<wire<<endl;
        if(mdcHit->driftTime(m_bunchT0,0)>1000)continue;
        if(trackType==circleType && mdcHit->layer()->view())continue;
        if(trackType==circleType && mdcHit->driftDist(m_bunchT0,0)>ddCut)continue;
        trkHitList->appendHot(newhot);
        //hitList_tds->push_back((*itHit));   
    }
    //---for CGEM clusters
    ClusterRefVec clusterRefVec = recMdcTrack->getVecClusters();
    ClusterRefVec::iterator itCluster = clusterRefVec.begin();
    for( ; itCluster != clusterRefVec.end(); itCluster++)
    {
        RecCgemCluster *recCgemCluster = *itCluster;
        CgemHitOnTrack* newhot = new CgemHitOnTrack(*recCgemCluster,m_bunchT0*1.e9);
        newhot->setCgemGeomSvc(m_cgemGeomSvc);
        newhot->setCgemCalibFunSvc(m_cgemCalibFunSvc);
        trkHitList->appendHot(newhot);
        if(m_debug>1)cout<<"layer: "<<recCgemCluster->getlayerid()<<"    flag: "<<recCgemCluster->getflag()<<"  clusterId: "<<recCgemCluster->getclusterid()<<endl;
    }
 
    //---fitting
    TrkHitList* qhits = trkRecoTrk->hits();   
    TrkErrCode err=qhits->fit();
    return err;
}

Referenced by execute().

initialize()

StatusCode CgemMdcFitAlg::initialize

(

)

Definition at line 76 of file CgemMdcFitAlg.cxx.

{
    MsgStream log(msgSvc(), name());
    log << MSG::INFO << "In CgemMdcFitAlg initialize()" << endreq;
 
    StatusCode sc;
    
    IMagneticFieldSvc* pIMF;
    sc = service ("MagneticFieldSvc",pIMF);
    if(sc!=StatusCode::SUCCESS) {
        log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
    }
    BField* bfield = new BField(pIMF);
    log << MSG::INFO << "field z = "<<bfield->bFieldNominal()<< endreq;
    m_context = new TrkContextEv(bfield);
 
    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;
    }
 
    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;
    }
 
    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;
    }
 
    m_mdcDetector = new MdcDetector();
    if(m_tuple) sc = bookTuple();
 
    return StatusCode::SUCCESS;
}

updateTracks()

void CgemMdcFitAlg::updateTracks	(	int	trackId,
		TrkRecoTrk *	trkRecoTrk,
		RecMdcTrack *	recMdcTrack,
		int	trkStat
	)

Definition at line 771 of file CgemMdcFitAlg.cxx.

                                                                                                         {
    
    //get the results of the fit to this track
    const TrkFit* fitresult = trkRecoTrk->fitResult();
    //check the fit worked
    if (fitresult == 0) return;
    
    //new a Rec. Track MdcTrack 
    //RecMdcTrack* recMdcTrack = new RecMdcTrack();
    const TrkHitList* aList = trkRecoTrk->hits();
    //some track info
    const BField& theField = trkRecoTrk->bField();
    double Bz = theField.bFieldZ();
    //Fit info
    int hitId = 0;
    int nHits=0;
    int nSt=0;
    //int nAct=0; int nSeg=0; 
    //int maxlay = 0; int minlay = 43; int seg[11];/* 11 slayer */ int segme;
    //****************************
    //* active flag open this
    //****************************
    //* if (allHit>0){
    //*   nHits = aList->nHit();//add inActive hit also
    //* }else{
    //*   nHits = fitresult->nMdc();// = nActive()
    //* }
    //****************************
    //for 5.1.0 use all hits
    nHits = aList->nHit();
    //****************************
    // use active only
    // nHits = fitresult->nMdc();// = nActive()
    //****************************
 
    int    q        = fitresult->charge();
    double chisq    = fitresult->chisq();
    int    nDof     = fitresult->nDof();
    //track parameters at closest approach to beamline
    double fltLenPoca = 0.0;
    TrkExchangePar helix = fitresult->helix(fltLenPoca);
    //std::cout<< __FILE__ << " phi0  " << helix.phi0() << " omega  " <<helix.omega()<<std::endl;
    double phi0 = helix.phi0();
    double tanDip = helix.tanDip();
    //double z0 = helix.z0();
    double d0 = helix.d0();
    //momenta and position
    //CLHEP::Hep3Vector mom = fitresult->momentum(fltLenPoca);
    HepPoint3D poca = fitresult->position(fltLenPoca);
 
    double helixPar[5];
    //dro =-d0
    helixPar[0] = -d0;
    //phi0 = phi0 - pi/2   [0,2pi)
    double tphi0 = phi0 - Constants::pi/2.;
    if (tphi0 < 0. ) tphi0 += Constants::pi * 2.;
    helixPar[1] = tphi0;
    //kappa = q/pxy
    double pxy = fitresult->pt();
    if (pxy == 0.)  helixPar[2] = 9999.;
    else helixPar[2] = q/fabs(pxy);
    if(pxy>9999.) helixPar[2] = 0.00001;
    //dz = z0
    helixPar[3] = helix.z0();
    //tanl 
    helixPar[4] = tanDip;
    //error
    //V(Y)=S * V(X) * ST , mS = S , mVy = V(Y) , helix.covariance() = V(X)
    HepSymMatrix mS(helix.params().num_row(),0);
    mS[0][0]=-1.;//dr0=-d0
    mS[1][1]=1.;
    mS[2][2]=-333.567/Bz;//pxy -> cpa
    mS[3][3]=1.;
    mS[4][4]=1.;
    HepSymMatrix mVy = helix.covariance().similarity(mS);
    double errorMat[15];
    int k = 0;
    for (int ie = 0 ; ie < 5 ; ie ++){
      for (int je = ie ; je < 5 ; je ++) {
    errorMat[k] = mVy[ie][je];
    k++;
      }
    }
    double p,px,py,pz;
    px = pxy * (-sin(helixPar[1]));
    py = pxy * cos(helixPar[1]);
    pz = pxy * helixPar[4];
    p = sqrt(pxy*pxy + pz*pz);
    //theta, phi
    double theta = acos(pz/p);
    double phi = atan2(py,px);
    recMdcTrack->setTrackId(trackId);
    recMdcTrack->setHelix(helixPar); 
    recMdcTrack->setCharge(q);
    recMdcTrack->setPxy(pxy);
    recMdcTrack->setPx(px);
    recMdcTrack->setPy(py);
    recMdcTrack->setPz(pz);
    recMdcTrack->setP(p);
    recMdcTrack->setTheta(theta);
    recMdcTrack->setPhi(phi);
    recMdcTrack->setPoca(poca);
    recMdcTrack->setX(poca.x());//poca
    recMdcTrack->setY(poca.y());
    recMdcTrack->setZ(poca.z());
    recMdcTrack->setR(sqrt(poca.x()*poca.x() + poca.y()*poca.y()));
    HepPoint3D pivot(0.,0.,0.);
    recMdcTrack->setPivot(pivot);
    recMdcTrack->setVX0(0.);//pivot 
    recMdcTrack->setVY0(0.);
    recMdcTrack->setVZ0(0.);
    recMdcTrack->setError(errorMat);
    recMdcTrack->setError(mVy);
    recMdcTrack->setChi2(chisq);
    recMdcTrack->setNdof(nDof);
    recMdcTrack->setStat(tkStat);
    recMdcTrack->setNhits(nHits);
    //-----------hit list-------------
    HitRefVec  hitRefVec;
    ClusterRefVec clusterRefVec;
    //vector<const RecCgemCluster*> clusterCol;
    vector<int> vecHits;
    //for fiterm
    int maxLayerId = -1;
    int minLayerId = 43;
    double fiTerm = 999.;
    const MdcRecoHitOnTrack* fiTermHot = NULL;
    TrkHitList::hot_iterator hot = aList->begin();
    for (;hot!=aList->end();hot++){
      //cout<<__FILE__ <<"  "<<__LINE__ <<"  "<< typeid(*hot).name()<<endl;
      if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
        const MdcRecoHitOnTrack* recoHot;             
        recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
        bool findhit = false;
        RecMdcHit* recMdcHit;
        Identifier mdcId = recoHot->mdcHit()->digi()->identify();
        HitRefVec vechits = recMdcTrack->getVecHits();
        HitRefVec::iterator itHit = vechits.begin();
        for( ; itHit != vechits.end(); itHit++)
        {
            Identifier mdcid = (*itHit)->getMdcId();
            if(mdcid == mdcId){
                recMdcHit = (*itHit);
                findhit = true;
                break;
            }
        }
        if(!findhit)continue;
        //if (!recoHot->mdcHit()) return;
        //id
        recMdcHit->setId(hitId);
        //---------BESIII----------
        //phiWire - phiHit <0 flagLR=0 left driftleft<0 ;
        //phiWire - phiHit >0 flagLR=1 right driftright>0;
        //flag = 2 ambig;
        //-----Babar wireAmbig----
        //-1 -> right, 0 -> don't know, +1 -> left
        //+1 phiWire-phiHit<0 Left,
        //-1 phiWire-phiHit>0 Right, 
        //0  don't know   
        //ambig() w.r.t track direction
        //wireAmbig() w.r.t. wire location 
        double hotWireAmbig = recoHot->wireAmbig();
        double driftDist = fabs(recoHot->drift());
        double sigma = recoHot->hitRms();       
        double doca = fabs(recoHot->dcaToWire());
        int flagLR=2;
        if ( hotWireAmbig == 1){
          flagLR = 0; //left driftDist <0
          doca *= -1.;//2012-07-19 
        }else if( hotWireAmbig == -1){
          flagLR = 1; //right driftDist >0
        }else if( hotWireAmbig == 0){
          flagLR = 2; //don't know 
        }
        recMdcHit->setFlagLR(flagLR);
        recMdcHit->setDriftDistLeft((-1)*driftDist);
        recMdcHit->setErrDriftDistLeft(sigma);
        recMdcHit->setDriftDistRight(driftDist);
        recMdcHit->setErrDriftDistRight(sigma);
        //Mdc Id
        recMdcHit->setMdcId(mdcId);
        //corrected ADC
 
        //contribution to chisquare
        //fill chisq
        double res=999.,rese=999.;
        if (recoHot->resid(res,rese,false)){
        }else{}
        double deltaChi=0;
        recoHot->getFitStuff(deltaChi);//yzhang open 2010-09-20 
        recMdcHit->setChisqAdd( deltaChi * deltaChi );
        //set tracks containing this hit
        recMdcHit->setTrkId(trackId);
        //doca
        recMdcHit->setDoca(doca);//doca sign left<0
        //entrance angle
        
        recMdcHit->setEntra(recoHot->entranceAngle());
        //z of hit 
        HepPoint3D pos; Hep3Vector dir;
        double fltLen = recoHot->fltLen(); 
        recMdcHit->setFltLen(fltLen);
        recoHot->getParentRep()->traj().getInfo(fltLen,pos,dir);
        recMdcHit->setZhit(pos.z());
        double tof = recoHot->getParentRep()->arrivalTime(fltLen);
        recMdcHit->setTdc(recoHot->mdcHit()->tdcIndex());
        recMdcHit->setAdc(recoHot->mdcHit()->adcIndex()); 
        //drift time      
        recMdcHit->setDriftT(recoHot->mdcHit()->driftTime(tof,pos.z()));
        //for fiterm
        int layerId = recoHot->mdcHit()->layernumber();
        int wireId = recoHot->mdcHit()->wirenumber();
        //std::cout<<" MdcTrack::store() ("<<layerId<<","<<wireId<<") fltLen "<<fltLen<<" wireAmbig  "<<hotWireAmbig<<"    y "<<pos.y()<<std::endl;
    //  <<recoHot->entranceAngle()<<std::endl;
        if (layerId >= maxLayerId){
          maxLayerId = layerId;
          fiTermHot = recoHot;
        } 
        if (layerId < minLayerId){
          minLayerId = layerId;
        }
        // status flag
        if (recoHot->isActive()) {
          recMdcHit->setStat(1);
        }else{
          recMdcHit->setStat(0);    
        }
        // for 5.1.0 use all hits
        if (recoHot->layer()->view()) {
          ++nSt;
        }
        //hitList->push_back(recMdcHit);
        SmartRef<RecMdcHit> refHit(recMdcHit);
        hitRefVec.push_back(refHit);        
        vecHits.push_back(mdcId.get_value());
        ++hitId;          
      }else if(typeid(*hot)==typeid(CgemHitOnTrack)){
        const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));            
        const RecCgemCluster* recCgemCluster = recoHot->baseHit();
        //cout<<"clusterid:"<< recCgemCluster->getclusterid()<<" layer:"<<recCgemCluster->getlayerid()<<endl;
        clusterRefVec.push_back(recCgemCluster);
        //clusterCol.push_back(recCgemCluster);
      }
    }
    std::sort(clusterRefVec.begin(),clusterRefVec.end(),sortCluster);
    //std::sort(clusterCol.begin(),clusterCol.end(),sortCluster);
    //fi terminal angle
    if (fiTermHot!=NULL){
      fiTerm=(1./sqrt(1.+tanDip*tanDip))*fiTermHot->fltLen()*helix.omega();
    }
    recMdcTrack->setFiTerm(fiTerm);
    // number of  stereo hits contained
    recMdcTrack->setNster(nSt);  
    recMdcTrack->setFirstLayer(maxLayerId);  
    recMdcTrack->setLastLayer(minLayerId);  
    recMdcTrack->setVecHits(hitRefVec);
    recMdcTrack->setVecClusters(clusterRefVec);
    recMdcTrack->setNcluster(clusterRefVec.size());
    //trackList->push_back(recMdcTrack);
 
    //cout<<__FILE__ <<"  "<<__LINE__ << endl;
    //for(int i=0;i<clusterRefVec.size();i++)
    //{
        //cout<<clusterRefVec[i]->getlayerid()<<endl;
    //}
    //hot = aList->begin();
    //for (;hot!=aList->end();hot++){
      //if(typeid(*hot)==typeid(MdcRecoHitOnTrack) || typeid(*hot)==typeid(MdcHitOnTrack) /*|| typeid(*hot)==typeid(TrkHitOnTrk)*/ ){
        //const MdcRecoHitOnTrack* recoHot;           
        //recoHot = dynamic_cast<const MdcRecoHitOnTrack*> (&(*hot));
        //int layerId = recoHot->mdcHit()->layernumber();
        //int wireId = recoHot->mdcHit()->wirenumber();
        //cout<<"hit ("<<layerId<<","<<wireId<<")"<<endl;
      //}else if(typeid(*hot)==typeid(CgemHitOnTrack)){
        //const CgemHitOnTrack* recoHot = dynamic_cast<const CgemHitOnTrack*> (&(*hot));          
        //const RecCgemCluster* recCgemCluster = recoHot->baseHit();
        //cout<<"hit ("<<recCgemCluster->getlayerid()<<","<<recCgemCluster->getclusterid()<<")"<<endl;
        //clusterRefVec.push_back(recCgemCluster);
      //}
    //}
}//end of storeTrack

Referenced by execute().

---

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

• CgemMdcFitAlg.h
• CgemMdcFitAlg.cxx