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();
        int                         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, double chiCut=7., int debug=0);// 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(int type=2);
        //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 updateLayerStat();
 
        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();
 
        bool nearStereoHits() {return m_nearStereoHits;};
 
        // --- for V-cluster selection with CGEM alone
        void clearXVhitsCgem();
        int findHelixInCgem();
 
 
    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;
        bool                        m_nearStereoHits;
 
        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; 
 
        // --- for V-cluster selection with CGEM alone
        //int m_nCgemLayerHasV;
        //vector<double> m_s_cgem[3];
        //vector<double> m_z_cgem[3];
        vector<HoughHit*> m_XVhits_cgem[3];
 
        // --- // 
};
 
typedef vector<HoughTrack*>::iterator TrackVector_Iterator;
#endif