DotsConnection.h

Go to the documentation of this file.

#ifndef  DotsConnection_H
#define  DotsConnection_H
 
#include "GaudiKernel/Algorithm.h"
#include "GaudiKernel/NTuple.h"
#include "MdcGeomSvc/MdcGeomSvc.h" 
#include "CgemGeomSvc/CgemGeomSvc.h"
#include "KalFitAlg/helix/Helix.h"
#include "McTruth/CgemMcHit.h"
#include "MdcRawEvent/MdcDigi.h"
#include "RawDataProviderSvc/IRawDataProviderSvc.h"
#include "RawDataProviderSvc/RawDataProviderSvc.h"
#include "MdcRecEvent/RecMdcTrack.h"
#include "DotsConnection/DotsHelixFitter.h"
#include "DotsConnection/MdcDigiGroup.h"
//#include "MdcUtilitySvc/IMdcUtilitySvc.h"
#include "MdcCalibFunSvc/MdcCalibFunSvc.h"
#include "TStopwatch.h"
#include "TH2D.h"
 
//using namespace KalmanFit;
 
// track candidate
// wire index (idx) >0 used in fit, <0 not used in fit, as in CGEM+ODC case, default idx >0
struct trkCandi
{
    int trkIdx;
    int isGood;// 1: good candidates, 0: not good, -1: combined to other candidate
 
    vector<int>    mdcHitIdx[43];
    //vector<vector<int> > mdcHitIdx(43);
 
    vector<int>    mdcXDigiWireIdx;// axial wires, negative if not associated
    vector<double> mdcXDigiChi;// fitted chi for axial wires
    vector<int>    mdcXDigiFitted;//0: not sure, 1: used in fit, -1: not used in fit
    int xLayerMin;
    int xLayerMax;
 
    vector<int>    mdcV1DigiWireIdx;// stereo 1 wires
    vector<int>    mdcV2DigiWireIdx;// stereo 2 wires
    vector<int>    mdcVDigiWireIdx; // stereo wires
    vector<double> mdcVDigiChi;     // fitted chi for stereo wires
    vector<int>    mdcVDigiFitted;  // if used in fit for stereo wires
    
    vector<RecMdcHit> RecMdcHitVec; // for RecMdcHit storage
 
    int layerMin;
    int layerMax;
    
    vector<int>    CgemXClusterID;// <0: inactive in fit
    vector<int>    CgemVClusterID;// <0: inactive in fit
    vector<double> CgemXClusterChi; // fitted chi 
    vector<double> CgemVClusterChi; // fitted chi 
 
    double      Smax;// max arc length (up to broken point)
    double      chi2;
    double      a_helix[5];
    //vector<double> a_helix(5);
    double phi_inner;
    double phi_outer;
    
    // --- for (fat) path finding
 
 
    bool operator==(const trkCandi& rhs) const{
 
 
        const trkCandi& lfs = *this;
        bool eq = true;
        bool eq_thisline ;
        #define COMPARE_A_ITEM__(NAME)                                                          \
        {                                                                                       \
            eq_thisline= (lfs.NAME == rhs.NAME);                                        \
            cout<<"compare ."<<#NAME<< " status:" <<(eq_thisline?"true":"false"); \
            cout << "\n" ;   \
            eq = eq && eq_thisline;                                                         \
        }  
        #define PRINT_IF_FAIL(NAME)\
        if  (!eq_thisline) cout<<"  because {lfs." #NAME ":rhs."#NAME "}: {"<<lfs.NAME<<":"<<rhs.NAME<<"}"<<endl;                                                                       
 
        COMPARE_A_ITEM__(trkIdx);
        COMPARE_A_ITEM__(isGood);
        COMPARE_A_ITEM__(mdcXDigiWireIdx);
        COMPARE_A_ITEM__(mdcXDigiChi);
        COMPARE_A_ITEM__(mdcXDigiFitted);
        COMPARE_A_ITEM__(xLayerMin);
        COMPARE_A_ITEM__(xLayerMax);
        COMPARE_A_ITEM__(mdcV1DigiWireIdx);
        COMPARE_A_ITEM__(mdcV2DigiWireIdx);
        COMPARE_A_ITEM__(mdcVDigiWireIdx); 
        COMPARE_A_ITEM__(mdcVDigiChi);     
        COMPARE_A_ITEM__(mdcVDigiFitted);  
        //COMPARE_A_ITEM__(RecMdcHitVec); //because RecMdcHitVec cannot be easily compared
        COMPARE_A_ITEM__(layerMin);
        COMPARE_A_ITEM__(layerMax);
        COMPARE_A_ITEM__(CgemXClusterID);
        COMPARE_A_ITEM__(CgemVClusterID);
        COMPARE_A_ITEM__(CgemXClusterChi); 
        COMPARE_A_ITEM__(CgemVClusterChi); 
        COMPARE_A_ITEM__(Smax);PRINT_IF_FAIL(Smax);
        COMPARE_A_ITEM__(chi2);PRINT_IF_FAIL(chi2);
 
        for (int i = 0; i < 5; i++){
            COMPARE_A_ITEM__(a_helix[i]);
            PRINT_IF_FAIL(a_helix[i]);
            
        }
        for (int i = 0; i < 43; i++){
            COMPARE_A_ITEM__(mdcHitIdx[i]);
            
        }
        return eq;
 
        #undef COMPARE_A_ITEM__
 
    }
 
    const trkCandi& operator=(const trkCandi& rhs){
        trkCandi& lfs = *this;
        #define COPY_A_ITEM__(NAME)                                                             \
        {                                                                                       \
            lfs.NAME = rhs.NAME;                                                                \
        }                                                                                       
        
        COPY_A_ITEM__(trkIdx);
        COPY_A_ITEM__(isGood);
        COPY_A_ITEM__(mdcXDigiWireIdx);
        COPY_A_ITEM__(mdcXDigiChi);
        COPY_A_ITEM__(mdcXDigiFitted);
        COPY_A_ITEM__(xLayerMin);
        COPY_A_ITEM__(xLayerMax);
        COPY_A_ITEM__(mdcV1DigiWireIdx);
        COPY_A_ITEM__(mdcV2DigiWireIdx);
        COPY_A_ITEM__(mdcVDigiWireIdx); 
        COPY_A_ITEM__(mdcVDigiChi);     
        COPY_A_ITEM__(mdcVDigiFitted);  
        COPY_A_ITEM__(RecMdcHitVec); 
        COPY_A_ITEM__(layerMin);
        COPY_A_ITEM__(layerMax);
        COPY_A_ITEM__(CgemXClusterID);
        COPY_A_ITEM__(CgemVClusterID);
        COPY_A_ITEM__(CgemXClusterChi); 
        COPY_A_ITEM__(CgemVClusterChi); 
        COPY_A_ITEM__(Smax);
        COPY_A_ITEM__(chi2);
        COPY_A_ITEM__(phi_inner);
        COPY_A_ITEM__(phi_outer);
 
        for (int i = 0; i < 5; i++){
            COPY_A_ITEM__(a_helix[i]);
            
        }
        for (int i = 0; i < 43; i++){
            COPY_A_ITEM__(mdcHitIdx[i]);
            
        }
        return lfs;
        #undef COPY_A_ITEM__
    }
 
};
 
struct dot{
    enum DetectorEnum{
        MDC =0, 
        CGEM=1,
        NUL_DETECTOR=-3
    };
    union DigiPtr
    { 
        const MdcDigi* mdcDigi;
        const RecCgemCluster* cgemDigi;
    };
 
    DetectorEnum detector;// 0: MDC, 1: CGEM
    DigiPtr digi;
    int stereo;// 0: axial, 1: stereo 1, -1: stereo 2
    int layer; // layer.
    int id ;// wire idx for MDC, cluster id for CGEM
    // double dd;// rough drift distance
    // std::vector<dot*> neighbors; // id+ or layer+
    // std::vector<dot*> backward_neighbors;  // id- or layer -
    // double phiLR;//left-right ambiguity for MDC, default 99
 
    vector<int> vecCellIdx;
    //vector<int> vecCellIdx_2ndDot;
    vector<int> vecGoodCellIdx[2];
    
    vector<int> trkId;
 
    int cluId;// hit-cluster index in an area
};// dot
 
static const int NUMBER_OF_DOTS_PER_CELL=3;
 
struct Line{
    double a, b; // ax + by + c =0
    double rho, theta; // rho=x*cos(theta)+y*sin(theta)
    //bool LR[NUMBER_OF_DOTS_PER_CELL];
    vector<int> LR;
    vector<double> S;
    vector<int> order;
    double chi2_hit[NUMBER_OF_DOTS_PER_CELL];
    double chi2;
};
 
 
class DotsConnection:public Algorithm {
    public:
        DotsConnection (const std::string& name, ISvcLocator* pSvcLocator);
        StatusCode initialize();
        StatusCode execute();
        StatusCode finalize();
        
        bool dotComp(dot* dot1, dot* dot2) {
            bool less=false;
            int layer1=myMdcGeomSvc->Wire(dot1->id)->Layer();
            int layer2=myMdcGeomSvc->Wire(dot2->id)->Layer();
            if(layer1<layer2) less=true;
            else if(layer1==layer2) {
                int wCell1=myMdcGeomSvc->Wire(dot1->id)->Cell();
                int wCell2=myMdcGeomSvc->Wire(dot2->id)->Cell();
                if(wCell1<wCell2) {
                    if(fabs(wCell1-wCell2)<0.5*myNWire[layer1]) less=true;
                }
                else if(fabs(wCell1-wCell2)>0.5*myNWire[layer1]) less=true;
            }
            return less;
        };
 
        struct cell{
            unsigned int id;
            unsigned int idInGood;
 
            dot* dots[NUMBER_OF_DOTS_PER_CELL];
            //double chi2;
 
            //Line fit_line[];
            vector<Line> vecFittedLine;
            int best_line_index;
 
            //double fitLine();
 
            set<int> set_dotIdx;// for dot-triplet
            
            set<int> set_inEdge; // for the index of in-edges  (good/valid edges)
            set<int> set_outEdge;// for the index of out-edges (good/valid edges)
            set<int> set_fwdCell;// for the index of forward-cells (possible edges tried)
 
            double state;
            int color;// 0: undiscovered,  1: discovered but under processing, 2: finished
 
            // --- get i-th dot idx
            int dotIdx(int i) {
                if(i>=0&&i<set_dotIdx.size()) {
                    set<int>::iterator it=set_dotIdx.begin();
                    while(i>0) { it++; i--; }
                    return *it;
                }
                return -9999;
            };
        };// cell
        vector<const MdcDigi*>  getMdcDigiVec(cell& aCell);
        int numSharedDots(cell& a, cell& b);
        bool sameCell(cell& a, cell& b);
        bool forward(cell& a, cell& b);
        bool isNeighbour(cell& a, cell& b);
        void print(cell& aCell);
        
        bool makeOrder(Line& l)
        {
            if(l.S.size()==0) return false;
            map<double, int> map_S_i;
            int n_minus(0), n_plus(0);
            vector<double>::iterator it_S=l.S.begin();
            //cout<<"  S= ";
            for(int i=0; it_S!=l.S.end(); it_S++,i++)
            {
                //cout<<setw(5)<<*it_S<<"  ";
                if((*it_S)<0) n_minus++;
                else n_plus++;
                map_S_i[fabs(*it_S)]=i;
            }
            //cout<<endl;
            if(n_minus==l.S.size()||n_plus==l.S.size()) 
            {
                map<double, int>::iterator it=map_S_i.begin();
                for(; it!=map_S_i.end(); it++) l.order.push_back((*it).second);
                return true;
            }
            else return false;
        };
 
        struct edge{
            unsigned id;
            pair<int,int> twoCells;// directed edge, with two cell indices in the good-cell-vector
            pair<int,int> twoLines;// index of the fitted lines
            set<int> set_dotIdx;// for dot-quadruplet
            enum dotStatus{
                Prev=0,
                Next=1,
                Common=2
            }; 
            vector<dot*> dots;
            vector<dotStatus> vDotStatus;
 
            vector<Line> vecFittedLine;
            int best_line_index;
        };
        bool makeEdge(cell& a, cell& b, edge& e);
        vector<const MdcDigi*>  getMdcDigiVec(edge& aEdge);
 
        int  myNCycle;
        vector<int> myCellMaxState;// cell index with the maximum state
        vector<vector<int> > CellAutomaton(vector<int>& v_cellIdx, vector<cell>& v_cell, vector<edge>& v_edge);// calculate states, find the (longest) path(s)
        void stateAutomaton(cell& aCell, vector<cell>& v_cell, vector<edge>& v_edge);
        void nextCellFinding(cell& aCell, vector<cell>& v_cell, vector<edge>& v_edge, map<int, pair<double,double> >& setWireHitPos, map<int, int>& setWireEntry);
        map<int, dot> myMapDots;
        vector<cell>  myVecCells;
        vector<int>   myVecGoodCellIdx;
        vector<edge>  myVecEdges;
        void makeDotCellEdge(vector<int>& vecWireIdx);
        void longestPathFromCA(vector<int>& vecWireIdx, vector<vector<int> >& v_path);
        void longestPathFromCA(set<int>& setWireIdx, vector<vector<int> >& v_path);
 
 
    private:
 
        // --- properties ----------
        int myDebug; // flag for debug output
        int myDebugNb;
        int myNtProd; // if produce Ntuple
        int myIniHelix; // initial helix in ideal track fitting, 1: from MC particle at IP, 2: from MC hit innermost
        double myDriftTimeUpLimit;
        double myMdcHitChi2Cut;
        bool myUseWireCrossPoint;
        //bool myUseIdealTracking;
        bool myWireCrossPointAvgMode;
        bool myWireCrossPointPersistent;
 
        bool myRunIdealTracking;
        bool myUseInTrkerHits;
        bool myRunTrkFollow;
 
        // --- event T0
        double getEventStartTime();// in ns
        double myEvtT0;
        
        // --- get MC truth
        //bool  getMCTruth();
        KalmanFit::Helix getMCHelix();
        int myMCparID;
 
        // --- get MC charged final states (e, mu, pi, K, p)
        void getMcFinalChargedStates();
        vector<int>        myVecMCTrkId;
        vector<int>        myVecPDG;
        vector<double>     myVecTrkLenFirstHalf;
        vector< vector<double> >  myVecHelix;
        vector<Event::CgemMcHit*> myVecCgemMcHit;
        vector<const RecCgemCluster*> myVecCgemXcluster;
        vector<const RecCgemCluster*> myVecCgemVcluster;
        vector<const RecCgemCluster*> myVecCgem1DCluster;
        vector<int> myVecCgemXCluIdx[3][2];
        vector<int> myVecCgemVCluIdx[3][2];
        vector<double> myVecCgemXCluChi[3][2];
        vector<double> myVecCgemVCluChi[3][2];
        void clearVecCgemClu() {
            for(int l=0; l<3; l++) {
                for(int m=0; m<2; m++) {
                    myVecCgemXCluIdx[l][m].clear();
                    myVecCgemVCluIdx[l][m].clear();
                    myVecCgemXCluChi[l][m].clear();
                    myVecCgemVCluChi[l][m].clear();
                }
            }
        };
        vector<const MdcDigi*>        myVecMdcDigi;
        void resetFCSVec();
        void associateDigisToMcParticles();
 
        // --- get digi
        RawDataProviderSvc*           myRawDataProviderSvc;
        const MdcDigi*                myMdcDigiPointer[43][288];// for ideal tracking
        void  clearMdcDigiPointer();
        vector<const MdcDigi*>        getMdcDigiVec();
        vector<const RecCgemCluster*> getCgemClusterVec(int view=0);// view = 0(X-view), 1(V-view), 2(X or V-view)
        double getRoughDD(const MdcDigi* aMdcDigi);
 
        // --- fill neighbours
        map<int, const MdcDigi*> myMapMdcDigi;// MDC digi map <index, digi> (all MDC digis)
        map<int, double>         myMapMdcDigiDd;// MDC digi map <index, dd> (all MDC digis)
        MdcDigiGroup  myMdcDigiGroup[6796];// for neighbours
        //MdcDigiGroup* myMdcDigiGroupPnt[43][288];
        vector<int>   myOuterEnds;
        vector<int>   myNeighbourSeeds;
        void clearMdcNeighbours();
        void fillMdcDigiBuildNeighbors(vector<const MdcDigi*> vecMdcDigiPnt, int SameLRange=1,int DiffLRange=1);
        void fillMdcDigiMap(vector<const MdcDigi*>& vecMdcDigiPnt);
        void buildMdcDigiNeighbors(map<int, const MdcDigi*>& aMapMdcDigi, int SameLRange=1,int DiffLRange=1);
        void findOuterEnds();
        //vector<vector<const MdcDigi*> > getMdcHitCluster();
        void getMdcHitCluster();
        //void buildMdcNeighbors();
        //myMdcDigiNeighbours;
        void mdcHitClustering(const vector<int>& vecWireIdx, vector<vector<int> >& vecHitClu);
 
        // --- calibration svc
        MdcCalibFunSvc*           myMdcCalibFunSvc;
        //ICgemCalibFunSvc*         myCgemCalibSvc;
 
        // --- geometry of MDC
        MdcGeomSvc* myMdcGeomSvc;
        int    myNWire[43];
        int    myWireFlag[43];
        int    myOuterWire[43][288][2]; 
        int    myInnerWire[43][288][2]; 
        double myWirePhi[43][288]; 
        double myRLayer[43];
 
        // --- useful functions
        //IMdcUtilitySvc*  myMdcUtilitySvc;
        double dPhi(double phi1, double phi2);// --- delta phi (-pi, pi)
 
        // --- for CGEM
        CgemGeomSvc*              myCgemGeomSvc;
        // --- for CGEM clusters
        //bool sortCluster(const RecCgemCluster* clusterA , const RecCgemCluster* clusterB);
        //SmartDataPtr<RecCgemClusterCol> myPntRecCgemClusterCol;
        RecCgemClusterCol::iterator     myIterCgemClusterBegin;
        RecCgemClusterCol::iterator     myIterCgemClusterEnd;
        RecCgemClusterCol::iterator     myIterClu;
        vector<int> myVecCgemXClusterIdx[3];// cluster index in RecCgemClusterCol [3 layers]
        vector<int> myVecCgemVClusterIdx[3];// [3 layers]
        vector<int> myVecCgemXClusterTrkIdx[3];// trk index in myVecTrkCandidates [3 layers]
        vector<int> myVecCgemVClusterTrkIdx[3];// [3 layers]
        bool getCgemClusters();// fill the above data
        bool getCgemClusterIntersection(int idx_Xclu, int idx_Vclu, double& z);
        void findCgemTrkSeg();
 
        // --- area index
        int myAreaLayer[43];// total 8 (7) areas: 4, 4, 12, 4, 4, 4, 4, 7 (3, 12, 4, 4, 4, 4, 7) for MDC (MDC+CGEM)
 
        // --- global circle/Helix fitter
        DotsHelixFitter myDotsHelixFitter;
        double myChi2CutDiverge;
        void   testDotsHelixFitterAllHits();
        void   testDotsHelixFitterPartHits();
        double myChiCut_circle;
        int    myNmaxDeact_circle;
        double myChiCut_helix;
        int    myNmaxDeact_helix;
        int    myNOuterXHits;
 
        // --- vector for track candidates
        vector<struct trkCandi> myVecTrkCandidates;// for a event
        vector<struct trkCandi> myVecTrkPaths;// for a track candiates with multi-paths
        vector<pair<int, int> > myVecNhitPathId;// for a track candiates with multi-paths [number-of-act-hit, tracknum]
        // --- tmp vector for track candidates fitting
        vector<int> myVecDigiIdx;
        vector<int> myVecCluIdx;
        // --- functions for track candidates
        void printTrkCandi(struct trkCandi& aTrkCandi);
        int  nXHitsActive(struct trkCandi& aTrkCandi);
        int  setFitFlagUncertain(struct trkCandi& aTrkCandi, int layerMin=-1, int layerMax=50);
        int  setXHitsInBigSetFitFlagUncertain(struct trkCandi& aTrkCandi, bool set=true);
        vector<const MdcDigi*>  getXDigiVec(struct trkCandi& aTrkCandi, int sel=0, vector<int>* vecFitFlag=NULL);// sel=0: all, sel=1: only associated hits, sel=2: only fitted hits & uncertain hits
        //vector<const MdcDigi*>  getXDigiVec(struct trkCandi& aTrkCandi, int sel, vector<int>& vecFitFlag);// sel=0: all, sel=1: only associated hits, sel=2: only fitted hits & uncertain hits
        vector<const MdcDigi*>  getVDigiVec(struct trkCandi& aTrkCandi, int sel=0);// sel=0: all, sel=1: only associated hits, sel=2: only fitted hits
        vector<const MdcDigi*>  getDigiVec(struct trkCandi& aTrkCandi, int sel=0);// sel=0: all, sel=1: only associated hits, sel=2: only fitted hits
        vector<const RecCgemCluster*>  getXCluVec(struct trkCandi& aTrkCandi, int sel=0);// sel=0: all, sel=1: only associated hits, sel=2: only fitted hits
        vector<const RecCgemCluster*>  getVCluVec(struct trkCandi& aTrkCandi, int sel=0);// sel=0: all, sel=1: only associated hits, sel=2: only fitted hits
        vector<const RecCgemCluster*>  getCluVec(struct trkCandi& aTrkCandi, int sel=0);// sel=0: all, sel=1: only associated hits, sel=2: only fitted hits
        void tagRedudantHits(struct trkCandi& aTrkCandi);
        bool split(struct trkCandi& aTrkCandi, int idxCandi);
        void tagClusterHits(struct trkCandi& aTrkCandi);// set fit flag = 0 for hits in cluster (> 3*3 hits)
        bool isInOutNeighbours(struct trkCandi& innerTrkCandi, struct trkCandi& outerTrkCandi, int gap=1);
        bool isInOutNeighbours(const vector<int>& innerSet, const vector<int>& outerSet);
        void combineTrkCandi(struct trkCandi& trkCandi_1, struct trkCandi& trkCandi_2);
        void rmPathInTrkCandi(struct trkCandi& aTrkCandi, struct trkCandi& path);
        void resetXDigiFitItem(struct trkCandi& aTrkCandi);
        void resetVDigiFitItem(struct trkCandi& aTrkCandi);
        void updateDigiFitItem(struct trkCandi& aTrkCandi);
        void updateCgemFitItem(struct trkCandi& aTrkCandi, int assignCgemClu=0);
        void dropRedundantMdcXDigi(struct trkCandi& aTrkCandi);
        void calcuCiclePar3CgemClu(struct trkCandi& aTrkCandi);
        std::vector<std::pair<double, double> > findWireCrossPoint(struct trkCandi& aTrkCandi, vector<double>* ref_track_param=nullptr);
        int  circleFitTaubin(struct trkCandi& aTrkCandi, double chiCut=9999.0, bool useIniHelx=false, int hitSel=1, bool setHitFitFlag=false, int layerMin=-1, int layerMax=50, bool use_additional_points=false, bool usePhiAmbi=false, const vector<double>& vPhiAmb=vector<double>() );
        void circleFit(struct trkCandi& aTrkCandi);
        void matchCgemClusterMdcDigi();
        void processTrkCandi();
        void processMdcHitCluter();
        bool completeTrkCandi(struct trkCandi& aTrkCandi, int idxCandi);
        bool matchMdcDigi(struct trkCandi& aTrkCandi, int idx);
        bool matchCgemXCluster(struct trkCandi& aTrkCandi);
        bool selectVHits(struct trkCandi& aTrkCandi);
        int  helixFit(struct trkCandi& aTrkCandi, double chiCut=9999.0, int layerMin=-1, int layerMax=50);
        int  saveGoodTrkCandi();
        void fillVecCgemClu(struct trkCandi& aTrkCandi);
        int  NXcluToDrop(struct trkCandi& aTrkCandi);
        bool isGoodCgemCircleCandi(struct trkCandi& aTrkCandi);
        struct trkCandi getTrkCandi(const vector<int>& vecWireIdx);
 
        vector<int> getUncommonVec(vector<int> v1, vector<int>& v2);
        void setMinusVec(set<int>& s, vector<int>& v);
 
        // --- Hough maps for sz vs tanL
        int    myNBinTanl;
        int    myNBinDz;
        double myTanlRange;
        double myDzRange;
        TH2D   myRoughTanlDzMap;
        int    fillTanlDzMap(double s, double z, int vote=3);
        void   getWeightedPeak(TH2D& h, double& x_peak, double& y_peak, double& x_weight, double&y_weight, int x_ext=2, int y_ext=1);
 
        // --- hough for line rho vs theta
        double myRhoRange;
        double myThetaRange;
        void   getRhoTheta_bisection(vector<const MdcDigi*>& aVecMdcDigi);// bisection method
        TH2D   myRoughRhoThetaMap;
 
        // --- total least line fit for MDC hits
        bool   LineFit_TLS(vector<const MdcDigi*>& aVecMdcDigi, vector<double>& vRho, vector<double>& vTheta, vector<vector<int> >& vvLR, vector<vector<double> >& vvS, bool ini=false, double rho_ini=0, double theta_ini=0, const vector<int>& LR_ini=vector<int>() );
        bool   isGood_Drho2_line(double drho2_tls, int layer=40) {
            if(drho2_tls<0.025) return true;
            else return false;
        }
        double getS(double rho, double theta, double x, double y) {
            if(rho<0) {
                rho*=-1.0;
                theta+=acos(-1.);
            }
            double r=sqrt(x*x+y*y);
            double phi=atan2(y,x);
            double s=r*sin(theta-phi);
            return s;
        }
        void getXY(double rho, double theta, double s, double& x, double& y) {
            if(rho<0) {
                rho*=-1.0;
                theta+=acos(-1.);
            }
            x=rho*cos(theta)+s*sin(theta);
            y=rho*sin(theta)-s*cos(theta);
        }
 
        // --- caculation of the 8 circles tangent to three given drift circles
        vector<vector<double> > getTangentCircles(vector<const MdcDigi*>& aVecMdcDigi);
 
        // --- the tangent circle with the maximum probability
        vector<double> getBestTangentCircle(vector<const MdcDigi*>& aVecMdcDigi);
 
        // --- RecMdcTrackCol
        RecMdcTrackCol* myRecMdcTrackCol;
        RecMdcHitCol*   myRecMdcHitCol;
        bool registerRecMdcTrack();
 
        // --- store RecMdcTrack from myDotsHelixFitter
        bool saveARecMdcTrack(int iTrk);
        bool saveARecMdcTrack(struct trkCandi& aTrkCandi);
 
        // --- event information
        int myRun;
        int myEvt;
 
        // --- timing tools
        Bool_t     myWatch_reset;
        TStopwatch myWatch_tot;
        TStopwatch myWatch_finding;
        TStopwatch myWatch_fitting;
        
        // --- root items for myDotsHelixFitter
        NTuple::Tuple*        myNtHelixFitter;
        NTuple::Item<int>     myRUN;
        NTuple::Item<int>     myEVT;
        NTuple::Item<int>     myPID;
        NTuple::Item<int>     myNPar;
        NTuple::Item<int>     myNXHits;
        NTuple::Item<int>     myNVHits;
        NTuple::Item<int>     myNXCluster;
        NTuple::Item<int>     myNVCluster;
        //NTuple::Item<int>     myInnerVHits;
        NTuple::Array<double> myArrayHelixMC;
        NTuple::Array<double> myArrayHelixMCInDec;
        NTuple::Array<double> myArrayHelixFitted;
        NTuple::Item<int>     myNHitsCircle;
        NTuple::Array<double> myLayerHitsCircle;
        NTuple::Array<double> myChiHitsCircle;
        NTuple::Item<int>     myNHits;
        NTuple::Array<double> myLayerHits;
        NTuple::Array<double> myChiHits;
        NTuple::Item<int>     myTrkIdBest;
        NTuple::Item<int>     myNHitsBestTrk;
        NTuple::Item<int>     myNSameHitsBestTrk;
 
        // --- root items for myNtTrkSeg
        NTuple::Tuple*        myNtTrkSeg;
        NTuple::Item<int>     myNt_nhits;
        NTuple::Array<int>    myNt_layer;
        NTuple::Array<int>    myNt_wire;
        NTuple::Array<double> myNt_Xhit;
        NTuple::Array<double> myNt_Yhit;
        NTuple::Array<double> myNt_DDhit;
        NTuple::Item<double>  myNt_rho;
        NTuple::Item<double>  myNt_theta;
};
 
#endif