HoughTrack.h

Go to the documentation of this file.

#ifndef HOUGHTRACK2_H
#define HOUGHTRACK2_H
 
#include "HoughTransAlg/HoughHit.h"
#include "TrackUtil/Helix.h"
//#include "KalFitAlg/helix/Helix.h"
#include "TrkBase/TrkErrCode.h"
#include "MdcGeom/MdcDetector.h"
#include "TrkFitter/TrkContextEv.h"
#include "MdcData/MdcHit.h"
#include "CgemData/CgemHitOnTrack.h"
#include "TrkBase/TrkRecoTrk.h"
#include "TGraph.h"
#include <set>
#include <map>
#include "TMinuit.h"
#include "DotsConnection/DotsHelixFitter.h"
 
using namespace std;
 
class HoughHit;
class HoughTrack:public Helix
{
  public:
    HoughTrack(int charge, double angle, double rho, double dAngle, double dRho, int trkID);
    HoughTrack(int charge, const HepPoint3D & position, const Hep3Vector & momentum, int trkID);
    HoughTrack(HepPoint3D& pivot, HepVector& a, int trkID);
    //HoughTrack(Helix& helix, int trkID);
    HoughTrack(const HoughTrack& other);
    HoughTrack& operator=(const HoughTrack& other);
    HoughTrack();
    //~HoughTrack();
            
    void                        setTrkID(int trkID) {m_trkID = trkID;}
    void                        setFlag(int flag) {m_flag = flag;}
    void                        setCharge(int charge) {m_charge = charge;}
    void                        setAngle(double angle) {m_angle = angle;}
    void                        setRho(double rho ) {m_rho = rho;}
    void                        setDAngle(double dAngle){m_dAngle = dAngle;}
    void                        setDRho(double dRho){m_dRho = dRho;}
    void                        setDTanl(double dTanl){m_dTanl = dTanl;}
    void                        setDDz(double dDz){m_dDz = dDz;}
    void                        setChi2(double chi2){m_chi2 = chi2;}
    //void                        setHelix(Helix helix) {m_helix = helix;}
    //void                        setHot(vector<HoughHit> hot) {m_hot = hot;}
    void                        setDz(double dz){m_dz = dz;}
    void                        setTanl(double tanl){m_tanl = tanl;}
 
    int                         getTrkID() const {return m_trkID;}
    int                         getCharge() const {return m_charge;}
    int                         getFlag() const {return m_flag;}
    double                      getAngle() const {return m_angle;}
    double                      getRho() const {return m_rho;}
    double                      getDAngle() const {return m_dAngle;}
    double                      getDRho() const {return m_dRho;}
    double                      getDTanl() const {return m_dTanl;}
    double                      getDDz() const {return m_dDz;}
    double                      getChi2() const {return m_chi2;}
    int                         getCircleFitStat() const {return m_circleFitStat;};
    //Helix*                      getHelix() {return &m_helix;}
    //vector<HoughHit>*           getHot() {return &m_hot;}
    TrkRecoTrk*                 getTrkRecoTrk() {return m_trkRecoTrk;}
    double                      getDr() const {return m_dr;}
    double                      getPhi0() const {return m_phi0;}
    double                      getKappa() const {return m_kappa;}
    double                      getDz() const {return m_dz;}
    double                      getTanl() const {return m_tanl;}
    //double                      angle() const {return m_helix.center().phi();}
    //double                      rho() const {return 1./m_helix.radius();}
    //double                      angle() const {return center().phi();}
    //double                      rho() const {return 1./radius();}
 
    //void                        addHot(HoughHit& hit){m_hot.push_back(hit);}
    //void                        dropHot(HoughHit& hit);
    //void                        dropHot(int hitID);
    //int                         findHot(vector<HoughHit*>& hitList, int flag);
    //int                         findHot(vector<HoughHit*>& hitList, vector<HoughHit*>& hot, int flag);
    //int                         findHot(vector<HoughHit>& hitList, vector<HoughHit*>& hitPntList, int flag, int charge=0);// llwang
    int                         findXHot(vector<HoughHit*>& hitList, int charge);// llwang
    //void                        sortHot();
    void                        sortHot(vector<HoughHit*>& hotList);
    
    vector<HoughHit*>           getVecHitPnt() {return m_vecHitPnt;}; // llwang 2018-12-26
    void                        addHitPnt(HoughHit* aHitPnt) {m_vecHitPnt.push_back(aHitPnt);};// llwang 2018-12-26
    void                        dropHitPnt(HoughHit* aHitPnt); 
    void                        dropVHitPnt(HoughHit* aHitPnt); 
    //void                        sortHitPnt();// sort hits, llwang 2019-01-02
    int                         getNHitsShared();
    void                        dropRedundentCgemXHits();
    void                        dropRedundentCgemVHits();
 
    int                         getNhitFirstHalf() {return m_nHitFirstHalf;};
    int                         getNhitSecondHalf() {return m_nHitSecondHalf;};
    int                         getNhitUnusedFirstHalf() {return m_nHitUnused_FirstHalf;};
    int                         getNhitUnusedSecondHalf() {return m_nHitUnused_SecondHalf;};
    void                        resetNhitHalf() {m_nHitFirstHalf=0; m_nHitSecondHalf=0; m_nHitUnused_FirstHalf=0; m_nHitUnused_SecondHalf=0;};
 
    double                      driftDistRes(HoughHit* hit);
    int                         judgeHalf(HoughHit* hit);
    int                         judgeCharge(HoughHit* hit);
    int                         judgeCharge(double xHit, double yHit);
    //double                      flightLength(HoughHit& hit);
    TrkErrCode                  fitCircle(const MdcDetector* mdcDetector, TrkContextEv* trkContextEv, double bunchT0);
 
    //void                        setDotsHelixFitter(DotsHelixFitter* fitter);
    int                         fitCircle(DotsHelixFitter* fitter, double bunchT0);// circle fit with DotsHelixFitter
    //vector<MdcDigi*>            getVecMdcDigi(vector<HoughHit*> aHoughHitVec);
    //vector<RecCgemCluster*>     getVecCgemCluster(vector<HoughHit*> aHoughHitVec);
    int             fitHelix(DotsHelixFitter* fitter, double bunchT0, RecCgemClusterCol::iterator recCgemClusterColBegin, double averageChi2cut=25);// helix fit with DotsHelixFitter
 
 
    //TrkErrCode                  fitHelix(const MdcDetector* mdcDetector, TrkContextEv* trkContextEv, double bunchT0, vector<MdcHit*> mdcHitCol);
    int                         calculateZ_S(HoughHit* hit);
    void                        updateHelix();
    void                        update(double angle, double rho);//llwang
    void                        updateCirclePar(double dr, double phi0, double kappa);//llwang
    void                        markUsedHot(vector<HoughHit*>& hitPntList, int use=1);
    void                        markUsedHot(int use=1);
    void                        print();
    void                        printHot();
    //HepVector                   calCircle(vector<HoughHit*>& aVecHoughHit, double& chi2);
    
    //vector<double>              CircleFitByTaubin(vector< pair<double, double> > pos_hits);// llwang 2018-12-22
    //vector<double>              fitCircle_Taubin(vector<HoughHit*>& aVecHoughHit);// llwang 2018-12-22
    //double                      chi2_circle(vector<HoughHit*>& aVecHoughHit);// llwang 2018-12-22
 
    void                        clearHits();
    bool                        isAGoodCircle();
    int                         getNTried() {return m_untried;};// llwang 2018-12-23
 
    void                        releaseSelHits(); // release m_vecHitPnt, llwang 2018-12-27
 
    vector<double>              getVecHitRes(){return m_vecHitResidual;}
    vector<HoughHit*>           getVecStereoHitPnt() {return m_vecStereoHitPnt;}; // llwang 2018-12-26
    void                        addVecStereoHitPnt(HoughHit* aHitPnt) {m_vecStereoHitPnt.push_back(aHitPnt);};
    vector<double>              getVecStereoHitRes(){return m_vecStereoHitRes;}
    TrkErrCode                  fitHelix(const MdcDetector* mdcDetector, const BField* bField, double bunchT0, vector<HoughHit*>& hot,int Layer);
    TrkErrCode                  fitHelix(const MdcDetector* mdcDetector, TrkContextEv* trkContextEv, double bunchT0, vector<MdcHit*>& mdcHitCol, vector<HoughHit*>& hot);
    int                         findVHot(vector<HoughHit*>& hitList, int charge, int maxLayer, double cutFactor=1.0);
    vector<HoughHit*>           getHotList(int type=2);
    HoughTrack*                 getMcTrack()const{return m_mcTrack;}
    void                        setMcTrack(HoughTrack* mcTrack){m_mcTrack = mcTrack;}
    void                        clearMemory();
 
    static int                  m_clearTrack;
    static int                  m_useCgemInGlobalFit;//0:not at all, 1:circle fit, 2:circle fit + final helix fit, >=3:circle fit + all helix fit
    static TGraph*              m_cut[2][43];
    //void setMdcHit(const vector<MdcHit*>* mdchit)  {vec_mdcHit= mdchit;}
 
    void setRecMdcHitVec(vector<RecMdcHit>& aRecMdcHitVec) { m_vecRecMdcHit=aRecMdcHitVec; };
    vector<RecMdcHit>* getRecMdcHitVec() {return &m_vecRecMdcHit;};
    void printRecMdcHitVec();
            
  private:
    int                         m_trkID;
    int                         m_charge;
    double                      m_angle;// beta or phi
    double                      m_rho;
    int                         m_flag;
 
    double                      m_dAngle;// bin width
    double                      m_dRho;// bin width
    double                      m_dTanl;// bin width
    double                      m_dDz;// bin width
 
    double                      m_dr;
    double                      m_phi0;
    double                      m_kappa;
    double                      m_dz;
    double                      m_tanl;
 
    double                      m_chi2;
 
    //Helix m_helix;
    //vector<HoughHit>            m_hot;
    TrkRecoTrk*                 m_trkRecoTrk;
    
    // --- added by llwang //
    int                         m_circleFitStat;
    vector<HoughHit*>           m_vecHitPnt;
    vector<double>              m_vecHitResidual;
    vector<double>              m_vecHitChi2;
    map<int,double>             m_map_lay_d; 
    map<int,HoughHit*>          m_map_lay_hit;
    std::set<int>               m_setLayer;
    int XGapSize(std::set<int> aLaySet, int& gapMax);// gap statistics, remove the last few hits after a big gap
    void XhitCutWindow(double rho, int ilayer, double charge, double& cut1, double &cut2);
    static TGraph*              m_cut1_cgem, *m_cut2_cgem;                                                        
    static TGraph*              m_cut1_ODC1, *m_cut2_ODC1;                                                        
    static TGraph*              m_cut1_ODC2, *m_cut2_ODC2;
    
    int                         m_unused;
    int                         m_untried;
    int                         m_nHitFirstHalf;
    int                         m_nHitUnused_FirstHalf;
    int                         m_nHitSecondHalf;
    int                         m_nHitUnused_SecondHalf;
    int                         m_nGap;
    int                         m_maxGap;
    //bool    m_;
 
    //TMinuit*                   m_pMnTrk;
 
    vector<HoughHit*>           m_vecStereoHitPnt;
    vector<double>              m_vecStereoHitRes;
    HoughTrack*                 m_mcTrack;
    vector<CgemHitOnTrack*>     m_cgemHitVector;//store data in that created by "new"
    //const vector<MdcHit*>*      vec_mdcHit;
 
 
    // --- DotsHelixFitter
    //DotsHelixFitter*            m_DotsHelixFitterPnt;
    // --- input for DotsHelixFitter
    vector<const MdcDigi*>        m_vecMdcDigiPnt;
    vector<const RecCgemCluster*> m_vecCgemClusterPnt;
    // --- RecMdcHits for digis active in DotsHelixFitter
    vector<RecMdcHit>             m_vecRecMdcHit; 
 
 
    // --- // 
};
 
typedef vector<HoughTrack*>::iterator TrackVector_Iterator;
#endif