///////////////////////////////////////////////////////////////// // // Track Finding based on the Conformal Mapping + Hough // Author: Yuki Fujii (yfujii@ihep.ac.cn) // ///////////////////////////////////////////////////////////////// #ifndef THoughTrackFinder_hxx_seen #define THoughTrackFinder_hxx_seen #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include class IHoughTrackFinder : public IVTrackFinder { public: IHoughTrackFinder(const char*name, const char* title); virtual ~IHoughTrackFinder(); typedef std::tr1::array HelixTrackParam; IHoughTrackFinder* GetTrackFinder() {return this;} /// main task called in event-by-event void FindTrack(std::vector< COMET::IHandle > &track); /// This method is filled in derived classes COMET::IHandle Process(const COMET::IAlgorithmResult& input); /// load the event Int_t Load(COMET::ICOMETEvent& anEvent); /// save the event Int_t Save(COMET::ICOMETEvent* anEvent); /// called at the begin of run or else (should not be in event-by-event) Int_t Init(); /// called at the end of run or else (should not be in event-by-event) Int_t Finish(); /// Apply Conformal Mapping void ConformalMapping(const Double_t x, const Double_t y, Double_t &cmU, Double_t &cmV) {cmU = x/(x*x+y*y); cmV = y/(x*x+y*y); return;} /// Inverse Transformation from Hough domain parameters(theta, rho) to circle parameters(x0,y0,r) in Cartesian coordinate void HoughToCircleParam(const Double_t theta, const Double_t rho, std::vector& cparam) { cparam.push_back(cos(theta)/2./rho); cparam.push_back(sin(theta)/2./rho); cparam.push_back(fabs(1./2./rho)); return; } /// Fill 2D Hough Histogram and return combi hit positions in DS local coordinate Int_t FillHistogramPhiR(const COMET::IHitSelection* hits, std::vector& combiHitsLoc); /// 2-dimensional Peak Search Int_t PeakSearch2D(std::vector< std::vector >& circles); /// Make a track candidate std::vector MakeTrackCandidatesFromHits(const COMET::IHandle& hits); /// Remove hits that are isolated COMET::IHitSelection* RemoveIsolatedHits(const COMET::IHandle& hits); /// Filter the hits COMET::IHitSelection* HitFiltering(const std::vector& circle, const COMET::IHitSelection* hits, TVector3& posMean, TVector3& posRMS); /// Re-fine circle/helix parameters Bool_t RefineCircle(std::vector& circle, std::vector& cerror, const COMET::IHitSelection* hits, Bool_t helixMode = false); /// Create the combined hits Bool_t MakeCombiHits(const COMET::IHitSelection* hits, std::vector& combiHits); /// Calculate the distance between two wires Double_t DistanceBetweenWires(const TVector3& wireAEnd0, const TVector3& wireAEnd1, const TVector3& wireBEnd0, const TVector3& wireBEnd1, TVector3& crossPoint); private: /// Not in use now... TVector2 CalcPerpendicular2D(const TVector2 vec0, const TVector2 vec1, const TVector2 circle) { TVector2 norm = (vec1-vec0).Unit(); Double_t scale = ((vec0-circle).Mod2() - (vec0-circle)*(vec1-circle)) /(vec0-vec1).Mod(); return (vec0 + scale * norm); } /// Fit for combi hits void CircleCombiFitting(std::vector< std::vector >& circles); /// Do fitting using 2-D circle function Bool_t CircleFitting(std::vector& circle, std::vector& cerror, const COMET::IHitSelection* hits); /// Do fitting using 3-D hexix function Bool_t HelixFitting (std::vector& circle, std::vector& cerror, const COMET::IHitSelection* hits); /// 2D circle fitting for combi hits (fClusters) Double_t CircleFitCombiFunction(const Double_t *par) const; /// 2D circle fitting for hits (fHitSelection) Double_t CircleFitFunction(const Double_t *par) const; /// Helix fitting for hits (fHitSelection) Double_t HelixFitFunction(const Double_t *par) const; /// Circle filtering Bool_t CircleFiltering(const std::vector& circle, const TVector3& end0, const TVector3& end1, const Double_t& dDist, TVector3& wirePos); /// Helix filtering Bool_t HelixFiltering (const std::vector& circle, const TVector3& end0, const TVector3& end1, const Double_t& dDist, TVector3& wirePos); /// Get the xy position of wire at given z TVector3 GetWirePositionAtZ(const TVector3& end0, const TVector3& end1, const Double_t zpos); /// Calculate the helix direction at the starting point TVector3 CalcDirection(HelixTrackParam& helix); COMET::IHitSelection* GetHitSelection() const { return fHitSelection;} void PrepareOutputs(); void WriteOutputs(); private: /// parameters for Hough transformation & peak search TH2F* fHistPhiR; /// 2D histogram in Hough domain Bool_t fWeightedPS; /// Enable weighted peak search in Hough domain using number of neighboring hits TF1* fHough; /// Hough curve TF1* fFuncDeltaR; /// Function to check the delta R TF1* fFuncHelixDeltaR; /// Function to check the delta R for helix COMET::IHitSelection* fHitSelection; /// HitSelection used for the fitting Double_t fRadiusThreForCircleFit; /// Threshold to determine "good hits" in the circle fitting Double_t fRadiusThreForHelixFit; /// Threshold to determine "good hits" in the helix fitting Double_t fRadiusResoForFit; /// Resolution to calculate the chi-square in the fitting Int_t fNBinTheta; /// Number of bins for theta in Hough domain Int_t fNBinR; /// Number of bins for R in Hough domain Double_t fRangeR; /// R range in Hough domain [-fRangeR, fRangeR] Int_t fNMaxPeak; /// Number of maximum peak in Hough peak search Double_t fSigmaHough; /// Peak separation in Hough peak search Double_t fPeakThreHough; /// Peak threshold in Hough peak search /// parameters for clustering Bool_t fMakeCombiHits; /// Flag to make combination hits Int_t fMinCombiHitsInCircle; /// Minimum number of combi hits in a circle Double_t fMaxDeltaXYcls; /// Maximum distance in XY plane between two hits to make a cluster Double_t fMaxDeltaZcls; /// Maximum distance in Z direction between two hits to make a cluster Double_t fMaxAbsZcls; /// Maximum allowed |Z| to make a cluster Double_t fMinNumberCls; /// Minimum number of clusters in an above region Double_t fMaxDeltaRPhiclsW; /// Maximum distance in R-Phi plane between two hits to make a cluster (for wider region) Double_t fMaxDeltaZclsW; /// Maximum distance in Z direction between two hits to make a cluster (for wider region) Double_t fMinNumberClsW; /// Minimum number of clusters in an above region (for wider region) /// Threshold for curvature Double_t fMinimumRadius; /// Minimum radius in cm Double_t fMaximumRadius; /// Maximum radius in cm /// Parameters for isolated hits removal Bool_t fRemoveIsolatedHits; /// Flag to remove isolated hits Double_t fNeighborR; /// R threshold to define neighboring hits in cm Double_t fNeighborPhi; /// Phi threshold to define neighboring hits in radian /// Parameters for the selection criteria to choose the track candidate UInt_t fMinHitsPerCircle1st; /// Minimum number of hits beloniging to the circle after FirstFiltering UInt_t fMinHitsPerCircle2nd; /// Minimum number of hits beloniging to the circle after SecondFiltering /// Threshold for filtering std::vector fMinDeltaRFilter; /// Threshold of delta R for i-th filtering in cm Int_t fNumFilter; //! Count numbe of filtering /// Debugging purpose Bool_t fOutputHists; /// Flag to output the histograms into output file Int_t fNumOfEvents; /// Number of events to be output TFile* fOutFile; /// File to put the histograms TH1F* fDrawFrame; /// Draw Frame TH1F* fDrawFrameZR; /// Draw Frame in Z-R plane Bool_t fIsMC; /// Flag to check whether input IHitSelection is "mchits" or not TGraph* fMCHits; /// MC hits TGraph* fMCHitsZR; /// MC hits TGraph* fAllHitWires; /// All wires contain "hits" TGraph* fClusters; /// Cluster positions TGraph* fFilteredWires1st; /// Hit wires after 1st filtering TGraph* fFilteredWires2nd; /// Hit wires after 2nd filtering TGraph* fFilteredWires3rd; /// Hit wires after 3rd filtering std::vector fCirclesIni; /// Circles std::vector fCircles1st; /// Circles std::vector fCircles2nd; /// Circles std::vector fCircles3rd; /// Circles std::vector fHelixTrackParamArray; ROOT::Math::Functor fFunctorCircleCombi; ROOT::Math::Functor fFunctorCircle; ROOT::Math::Functor fFunctorHelix; protected: }; #endif