RecMdcTrack.cxx

Go to the documentation of this file.

/*
* @class : RecMdcTrack 
*  this class models "Track" used in Mdc Reconstruction                        
*                                                                              
*1.add constructor information of m_ncluster,m_vecclusters,m_matchChi2. 2016-09-28
* ********************************************************/
 
#include "MdcRecEvent/RecMdcTrack.h"
#include "Identifier/MdcID.h"
#include <set>
 
  const CLID& RecMdcTrack::clID() const {
    return RecMdcTrack::classID();
  }
                                                                               
  const CLID& RecMdcTrack::classID() {
    return CLID_RecMdcTrack;
  }
 
  RecMdcTrack::RecMdcTrack() 
  : m_pivot(0,0,0),
    m_fiTerm(0.0),           
    m_nhits(0),
    m_ncluster(0),
    m_matchChi2(9999.0)
   {   
   }
 
  RecMdcTrack::RecMdcTrack( const RecMdcTrack& rectrk)
  : DstMdcTrack(rectrk),
    m_pivot(rectrk.m_pivot),
    m_fiTerm(rectrk.m_fiTerm),
    m_nhits(rectrk.m_nhits),
    m_vechits(rectrk.m_vechits),
    m_ncluster(rectrk.m_ncluster),
    m_matchChi2(rectrk.m_matchChi2),
    m_vecclusters(rectrk.m_vecclusters)
  {
  }
 
  RecMdcTrack::RecMdcTrack( const DstMdcTrack& dstrk) 
  : DstMdcTrack(dstrk),
    m_pivot(-99.,-99.,-99.),
    m_fiTerm(-99.),
    m_nhits(-99)
  {   
  }
 
  RecMdcTrack& RecMdcTrack::operator=(const RecMdcTrack& rectrk) 
  {
    if( &rectrk != this ){
      DstMdcTrack::operator=(rectrk);
      m_pivot = rectrk.m_pivot;
      m_fiTerm = rectrk. m_fiTerm;
      m_nhits = rectrk. m_nhits;
      m_vechits = rectrk.m_vechits;
      m_ncluster = rectrk.m_ncluster;
      m_matchChi2 = rectrk.m_matchChi2;
      m_vecclusters = rectrk.m_vecclusters;
    }
    return  *this;
  }
 
  RecMdcTrack& RecMdcTrack::operator=(const DstMdcTrack& dstrk)
  {
    if( &dstrk != this ){
      DstMdcTrack::operator=(dstrk);  
      m_pivot[0] = -99.;
      m_pivot[1] = -99.;
      m_pivot[2] = -99.;
      m_fiTerm   = -99.;
      m_nhits    = -99;
    }
    return  *this;
  }  
 
  RecMdcTrack::~RecMdcTrack() {
  }               
 
void RecMdcTrack::setVecHits(HitRefVec vechits)    
{
    m_vechits = vechits;
    updateStatistics();
}
 
void RecMdcTrack::setVecClusters(ClusterRefVec vecclusters)     
{
    m_vecclusters = vecclusters;
    updateStatistics();
}
 
void RecMdcTrack::setVecClusters(ClusterRefVec vecclusters,map<int,int> clusterFitStat)
{
    m_vecclusters = vecclusters;
    m_clusterFitStat = clusterFitStat;
    updateStatistics();
}
 
void RecMdcTrack::updateStatistics()
{
    set<int> mdc_layers;
    set<int> cgem_layers;
 
    int nDof = 0;
    int nSter = 0;
    int maxLayerId = -1;
    int minLayerId = 43;
    int nLayer = 0;
    double fiTerm = 999.;
 
    MdcID* mdcID = new MdcID();
    int nHits = m_vechits.size();
    HitRefVec::iterator iterHit = m_vechits.begin();
    for( ; iterHit != m_vechits.end(); iterHit++){
        Identifier mdcid = (*iterHit)->getMdcId();
        int layerId = MdcID::layer(mdcid);
        //count nDof
        //cout<< (*iterHit)->getStat()<<endl;
        if ((*iterHit)->getStat()==1){
            nDof++;
        }
        //count nSter
        //cout<<layerId<<"  "<<mdcID->is_axial(mdcid)<<endl;
        if (!(mdcID->is_axial(mdcid))){
        //if (layerId >=20 &&layerId <36){
            nSter++;
        }
        // find maxLayerId and minLayerId
        if (layerId >= maxLayerId){
            maxLayerId = layerId;
            //fiTerm=(1./sqrt(1.+m_helix[4]*m_helix[4]))*(*iterHit)->getFltLen()*m_helix[3]/333.567;
        }
        if (layerId < minLayerId){
            minLayerId = layerId;
        }
        //count nlayer
        mdc_layers.insert(layerId);
        //count nHits
        //nHits++;
    }
 
    //cout<<__FILE__ <<"  "<<__LINE__ << endl;
    //cout<<"nDof "<<nDof-5<<endl;
    //cout<<"nSter "<<nSter<<endl;
    //cout<<"maxLayerId "<<maxLayerId<<endl;
    //cout<<"minLayerId "<<minLayerId<<endl;
    //cout<<"nHits "<<nHits<<endl;
    //cout<<endl;
 
    int nCgemXClusters=0;
    int nCgemVClusters=0; 
    int nCluster=m_vecclusters.size();
    ClusterRefVec::iterator iterCluster = m_vecclusters.begin();
    for( ; iterCluster != m_vecclusters.end(); iterCluster++){
        (*iterCluster)->setTrkId(trackId());
        int flag = (*iterCluster)->getflag();
        int clusterid = (*iterCluster)->getclusterid();
        map<int,int>::iterator iter = m_clusterFitStat.find(clusterid);
        int stat =1;
        if(iter != m_clusterFitStat.end())stat = iter->second;
        switch(flag){
            case 0:
                nCgemXClusters++;
                if(stat)nDof++;
                break;
            case 1:
                nCgemVClusters++;
                if(stat)nDof++;
                nSter++;
                break;
            case 2:
                nCgemXClusters++;
                nCgemVClusters++;
                if(stat)nDof+=2;
                nSter++;
                break;
            case 3:
                nCgemXClusters++;
                nCgemVClusters++;
                if(stat)nDof+=2;
                nSter++;
                break;
            default: break;
        }
        int layerId = (*iterCluster)->getlayerid();
        if (layerId >= maxLayerId){
            maxLayerId = layerId;
            //fiTerm=(1./sqrt(1.+m_helix[4]*m_helix[4]))*(*iterHit)->getFltLen()*m_helix[3]/333.567;
        }
        if (layerId < minLayerId){
            minLayerId = layerId;
        }
        cgem_layers.insert(layerId);
        nCluster++;
    }
 
    nLayer = mdc_layers.size() + cgem_layers.size();
 
    //cout<<__FILE__ <<"  "<<__LINE__ << endl;
    //cout<<"nDof "<<nDof-5<<endl;
    //cout<<"nSter "<<nSter<<endl;
    //cout<<"maxLayerId "<<maxLayerId<<endl;
    //cout<<"minLayerId "<<minLayerId<<endl;
    //cout<<"nCluster "<<nCluster<<endl;
    //cout<<"nLayer "<<nLayer<<endl;
    //cout<<endl;
 
    setNdof(nDof-5);
    setNster(nSter);
    setFirstLayer(minLayerId);
    setLastLayer(maxLayerId);
    setNlayer(nLayer);
    setNCgemXCluster(nCgemXClusters);
    setNCgemVCluster(nCgemVClusters);
    setNhits(nHits);
    setNcluster(nCluster);
 
    //cout<<"mdc nlayer:"<<mdc_layers.size()<<endl;
    //cout<<"cgem nlayer:"<<cgem_layers.size()<<endl;
}