Hough3D.cxx

Go to the documentation of this file.

#include "MdcHoughFinder/Hough3D.h"
#include "CgemData/CgemHitOnTrack.h"
int Hough3D::m_debug=0;
vector<MdcHit*> vec_for_clean;
vector<TrkRecoTrk*> vectrk_for_clean;
TrkContextEv* Hough3D::_context=NULL;
const MdcCalibFunSvc* Hough3D::_mdcCalibFunSvc=NULL;
const CgemCalibFunSvc* Hough3D::_cgemCalibFunSvc=NULL;
const CgemGeomSvc* Hough3D::_cgemGeomSvc=NULL;
int Hough3D::m_qualityFactor=3;
//double Hough3D::m_dropTrkDrCut=1;
//double Hough3D::m_dropTrkDzCut=10;
double Hough3D::m_dropTrkDrCut=1;        //no limit of dr dz
double Hough3D::m_dropTrkDzCut=10;
double Hough3D::m_dropTrkNhitCut=9;
double Hough3D::m_dropTrkChi2Cut=99999;
double Hough3D::m_dropTrkChi2NdfCut=99999;
//int    m_cgemfit = 1;
 
Hough3D::Hough3D(){
  //_m_gm = MdcDetector::instance(0);
}
Hough3D::Hough3D(const Hough3D& other){
  _circleR=other._circleR;
  _circleX=other._circleX;
  _circleY=other._circleY;
  _charge=other._charge;
  _d0=other._d0;
  _phi0=other._phi0;
  _omega=other._omega;
  _tanl=other._tanl;
  _z0=other._z0;
 
  _pt=other._pt;
  _pz=other._pz;
  _p=other._p;
  _nfit=other._nfit;
  _bunchT0=other._bunchT0;
  _chi2_aver=other._chi2_aver;
  _m_gm = other._m_gm;
  _recHitVec=other._recHitVec;
  vec_mdcHit=other.vec_mdcHit;
}
 
Hough3D::Hough3D(const Hough2D& hough2D,recHitCol hitCol, double bunchtime,double tanl,double z0):_tanl(tanl),_z0(z0){
  //_m_gm = MdcDetector::instance(0);
  _circleR=hough2D.getCirR();
  _circleX=hough2D.getCirX();
  _circleY=hough2D.getCirY();
  _d0 = hough2D.getD0();
  _phi0 = hough2D.getPhi0();
  _omega= hough2D.getOmega();
  _charge=hough2D.getCharge();
  _bunchT0=bunchtime;
  _recHitVec=hitCol;
  _pt = 0;
  _p = 0;
  _pz = 0;
  _nfit=0;
  newTrack=NULL;
}
Hough3D::Hough3D(const Hough2D& hough2D,recHitCol hitCol, double bunchtime,double tanl,double z0,const vector<MdcHit*>* mdchit):_tanl(tanl),_z0(z0),vec_mdcHit(mdchit){
  //_m_gm = MdcDetector::instance(0);
  _circleR=hough2D.getCirR();
  _circleX=hough2D.getCirX();
  _circleY=hough2D.getCirY();
  _d0 = hough2D.getD0();
  _phi0 = hough2D.getPhi0();
  _omega= hough2D.getOmega();
  _charge=hough2D.getCharge();
  _bunchT0=bunchtime;
  _recHitVec=hitCol;
  _pt = 0;
  _p = 0;
  _pz = 0;
  _nfit=0;
  newTrack=NULL;
}
 
int Hough3D::fit(){
  if(_charge==-1){                                                                      
//  cout<<" charge -1 "<<endl;
//  _d0=-_d0;                                                                               
//  _omega=-1.*fabs(_omega);                                                                   
  }                                                                                       
  if(_charge==1){                                                                       
//  cout<<" charge +1 "<<endl;
//  _d0=_d0;                                                                                
//  _omega=1.*fabs(_omega);
    //_phi0=_phi0-M_PI;                                                                       
  } 
 
  //outerHit();                // delete 4 hit in a wire  for kalman
 
  TrkExchangePar tt(_d0,_phi0,_omega,_z0,_tanl);
  //if (m_debug>0) cout<<"q d0 phi0 omega: "<<_charge<<" "<<_d0<<" "<<_phi0<<" "<<_omega<<" "<<_tanl<<" "<<_z0<<endl;
  double  dr    = -_d0;
  double  phi0  = _phi0-M_PI/2;
  double  kappa = _omega*333.567;
  double  dz    = _z0;
  double  tanl  = _tanl;
  while(phi0>2*M_PI)phi0-=2*M_PI;
  while(phi0<0)phi0+=2*M_PI;
  if (m_debug>0)cout<<"rec:"<<setw(12)<<dr<<"  "<<setw(12)<<phi0<<"  "<<setw(12)<<kappa<<"  "<<setw(12)<<dz<<"  "<<setw(12)<<tanl<<endl;
  //cout<<"rec barbar:"<<setw(12)<<_d0<<"  "<<setw(12)<<_phi0<<"  "<<setw(12)<<_omega<<"  "<<setw(12)<<_z0<<"  "<<setw(12)<<_tanl<<endl;
  //cout<<"rec besIII:"<<setw(12)<<dr<<"  "<<setw(12)<<phi0<<"  "<<setw(12)<<kappa<<"  "<<setw(12)<<dz<<"  "<<setw(12)<<tanl<<endl;
  TrkHelixMaker helixfactory;
  float chisum =0.;
  newTrack = helixfactory.makeTrack(tt, chisum, *_context, _bunchT0);
  //vectrk_for_clean.push_back(newTrack);
  bool permitFlips = true;
  bool lPickHits = true;
  helixfactory.setFlipAndDrop(*newTrack, permitFlips, lPickHits);
  TrkHitList* trkHitList = newTrack->hits();
  int digiId=0;
  std::vector<HoughRecHit>::iterator iter = _recHitVec.begin();
  for (;iter != _recHitVec.end(); iter++, digiId++) {
//cout<<__FILE__<<"  layer: "<<iter->getLayerId()<<"    wire: "<<iter->getWireId()<<"    flag:"<<iter->getflag()<<endl;
    if( (*iter).getflag()!=0 ) continue;
    //if( (*iter)->calDriftDist(_bunchT0,0)>0.8) continue; 
    //if( (*iter)->driftTime(_bunchT0,0)>800) continue;
    //cout<<"rechit dist time  "<<(*iter).calDriftDist(_bunchT0,0)<<" "<<(*iter).driftTime(_bunchT0,0)<<endl;
    if((*iter).getDetectorType()==MDC) {
        const MdcDigi* aDigi = (*iter).digi();
        MdcHit* hit;
        // compare
        vector<MdcHit*>::const_iterator iter_digi = (*vec_mdcHit).begin();
        for (;iter_digi != (*vec_mdcHit).end(); iter_digi++) {
          if( (*iter_digi)->digi() == (*iter).digi() )  {
            hit = (*iter_digi);
          }
        }
        Identifier id = aDigi->identify();
        int layer = MdcID::layer(id);
        int wire  = MdcID::wire(id);
        if(m_debug>0)cout<<"layer: "<<layer<<"    wire: "<<wire<<endl;
//cout<<"layer: "<<layer<<"    wire: "<<wire<<endl;
        hit->setCalibSvc(_mdcCalibFunSvc);
        hit->setCountPropTime(true); 
        // new fltLen and ambig
        int newAmbig = 0;
        // calc. position of this point
        MdcRecoHitOnTrack temp( *hit, newAmbig, _bunchT0*1.e9);//m_bunchT0 nano second here
        MdcHitOnTrack* newhot = &temp;
        double flight;
        double z_flight;
        //if( layer<8 ) flight = sqrt( ((*iter).getsPos())*((*iter).getsPos())+((*iter).getzPos())*((*iter).getzPos())) ;
        //  if( layer<8 ) flight =(*iter).getsPos();
        flight = (*iter).getRet().second;
        double distoTrack= (*iter).getDisToTrack();
        //  if( layer<8 ) newhot->setFltLen( (*iter).getsPos());             //caculate by mc
        //  else newhot->setFltLen( (*iter).getRet().second);             //caculate by mc
        newhot->setFltLen( flight );
        //  if(layer==18) continue;
 
        int use_in3d=1;
        //  if(hit->driftDist(_bunchT0,0)>1.0)  use_in3d=0; 
        if(hit->driftTime(_bunchT0,0)>1000) use_in3d=0; 
        //if(m_debug>0) cout<<"flt : "<<(*iter).getsPos()<<endl;
        //if(m_debug>0) cout<<"flt : "<<flight<<endl;
        //if(m_debug>0) std::cout<<" ("<<layer<<","<<wire<<",rt "<<hit->rawTime()<<",dt "<<hit->driftTime(_bunchT0,0)<<") T0 " << _mdcCalibFunSvc->getT0(layer,wire) <<" timewalk "<<  _mdcCalibFunSvc->getTimeWalk(layer, aDigi->getChargeChannel())<< "dist: "<<hit->driftDist(_bunchT0,0)<<" disToTrk "<< distoTrack<<"  use?: "<<use_in3d<<endl;
        //  if(hit->driftTime(_bunchT0,0)>800) continue; 
        if(use_in3d==0) continue;
        trkHitList->appendHot(newhot);
    }else if( (*iter).getDetectorType()==CGEM && (*iter).getCluster()->getflag()==2){
        if(m_debug>0) std::cout<<__FILE__<<"   "<<__LINE__<<" add CGEM hit to fit  "<<std::endl;
        CgemHitOnTrack* newhot = new CgemHitOnTrack(*((*iter).getCluster()),_bunchT0*1.e9);
        newhot->setCgemGeomSvc(_cgemGeomSvc);
        newhot->setCgemCalibFunSvc(_cgemCalibFunSvc);
        trkHitList->appendHot(newhot);
        //if(m_cgemfit==1 || m_cgemfit==3)trkHitList->appendHot(newhot);
    }
  }
  if( m_debug>0) newTrack->printAll(std::cout);
  /*
  TrkHitList* hitList = newTrack->hits();
  for (TrkHitList::hot_iterator ihit(hitList->begin()); ihit != hitList->end(); ++ihit){
      if(typeid(*ihit)==typeid(CgemHitOnTrack))continue;
      int nUsed = ihit->hit()->nUsedHits();
      std::cout<<"nUsed="<<nUsed<<":"<<endl;
      //ihit->hit()->printAll(std::cout);
      ihit->printAll(std::cout);
      std::pair<TrkFundHit::hot_iterator,TrkFundHit::hot_iterator> q = ihit->hit()->getUsedHits();
      for (TrkFundHit::hot_iterator i = q.first; i != q.second; ++i) {
          TrkRecoTrk * recoTrk=i->parentTrack();
          int id = recoTrk->id();
          cout<<"trkid:"<<id<<endl;
      }
      cout<<endl;
  }
  */
  TrkHitList* qhits = newTrack->hits();
  TrkErrCode err=qhits->fit();
  const TrkFit* newFit = newTrack->fitResult();
  int nActiveHit = newTrack->hots()->nActive();
  //cout<<err<<endl;
  //if(newFit) cout <<" nAct "<<newFit->nActive()<<"  "<<nActiveHit<<endl;
  int fit_stat=false;
  double chi2=-999.;
  if (!newFit || (newFit->nActive()<5)||err.failure()!=0){    
  //if (!newFit )
    if(m_debug>0){
      cout << " global 3d fit failed ";
      if(newFit) cout <<" nAct "<<newFit->nActive()<<endl;
      cout<<"ERR1 failure ="<<err.failure()<<endl;
      fit_stat=0;
    }
  }else{
    TrkExchangePar par = newFit->helix(0.);
    if( abs( 1/(par.omega()) )>0.03) {
      fit_stat = 1;
      chi2=newFit->chisq();
    }
    else {
      fit_stat = 0;
      chi2=-999;
    }
 
    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) {
      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 newTrack;
    //vectrk_for_clean.pop_back();
  }
    if( fit_stat==1 ){
    if(m_debug>0){
      cout << " global 3d fit success"<<endl;
      cout<<__FILE__<<__LINE__<<"AFTER fit 3d "<<endl;
      newTrack->print(std::cout);
    }
    TrkExchangePar par = newFit->helix(0.);
    _d0=par.d0();
    _phi0=par.phi0();
    _omega=par.omega();
    _tanl=par.tanDip();
    _z0=par.z0();
 
    double dr    = -par.d0();
    double phi0  = par.phi0()-M_PI/2;
    if(phi0<0)phi0+=2*M_PI;
    double kappa = par.omega()*333.567;
    double dz=par.z0();
    double tanl=par.tanDip();
    if(m_debug>0)cout<<"Fit:"<<setw(12)<<dr<<"  "<<setw(12)<<phi0<<"  "<<setw(12)<<kappa<<"  "<<setw(12)<<dz<<"  "<<setw(12)<<tanl<<endl<<endl;
    _circleR=_charge/par.omega();
    _circleX= sin(par.phi0()) *(par.d0()+1/par.omega()) * -1.; //z axis verse,x_babar = -x_belle
    _circleY= -1.*cos(par.phi0()) *(par.d0()+1/par.omega()) * -1;//???
    double pxy=fabs(_circleR/333.567);
    double pz=pxy*par.tanDip();
    double pxyz=_charge*sqrt(pxy*pxy+pz*pz);
    _pt=pxy; 
    _pz=pz;
    _p=pxyz;
    if(m_debug>0){
      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;
    if(m_debug>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)){
        if(m_debug>0){ 
          cout <<"(" <<((MdcHit*)(hotIter->hit()))->layernumber()
               <<","<<((MdcHit*)(hotIter->hit()))->wirenumber() 
               <<":"<<hotIter->isActive()<<")  ";
        }
      }else if(typeid(*hotIter)==typeid(CgemHitOnTrack)){
        if(m_debug>0){ 
          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++;
    }
    _nfit=nfit3d;
  }
  _chi2_aver=chi2/_nfit;
  //  for(int i=0;i<t_hitCol.size();i++){
  //    delete t_hitCol[i];
  //  }
  if(m_debug>0) cout<<" in 3D fit number of Active hits : "<<_nfit<<endl;
  return fit_stat;
}
 
void Hough3D::print(){
  cout<<" print hough3d "<<endl;
  cout<<"par: "<<_d0<<" "<<_phi0<<" "<<_omega<<" "<<_tanl<<" "<<_z0<<endl;
}
 
bool less_layer_3D(const int& a,const int& b){
  return a < b;
}
 
void Hough3D::outerHit(){
  for(int ihit=0;ihit<_recHitVec.size();ihit++){
    if( _recHitVec[ihit].calDriftDist(_bunchT0,0)>0.8 || _recHitVec[ihit].driftTime(_bunchT0,0)>800 ) _recHitVec[ihit].setflag(-5);   //flag -5 time or drift
  }
 
  vector<int> vec_layer_num[43];
  for(int ilay=0;ilay<43;ilay++){
    for(int ihit=0;ihit<_recHitVec.size();ihit++){
      if( _recHitVec[ihit].getLayerId()==ilay && _recHitVec[ihit].getflag()==0 ) {
        vec_layer_num[ilay].push_back( _recHitVec[ihit].getWireId() );
      }
    }
    std::sort( vec_layer_num[ilay].begin(),vec_layer_num[ilay].end(),less_layer_3D);
    if( vec_layer_num[ilay].size() < 4 ) continue;
    for(int j=0;j<vec_layer_num[ilay].size();j++) {
      //      if( (vec_layer_num[ilay][j]+1 == vec_layer_num[ilay][j+1]) &&  (vec_layer_num[ilay][j+1]+1 == vec_layer_num[ilay][j+2]) &&  (vec_layer_num[ilay][j+2]+1 == vec_layer_num[ilay][j+3]) ) {
      for(int jhit=0;jhit<_recHitVec.size();jhit++) {
        if( (ilay==_recHitVec[jhit].getLayerId()) && (vec_layer_num[ilay][j]==_recHitVec[jhit].getWireId() ) )   {_recHitVec[jhit].setflag(-1);}
      }
      //}
    }
  }
}