MpdTpcKalmanFilter.h 7.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151
  1. #ifndef MPDTPCKALMANFILTER_H
  2. #define MPDTPCKALMANFILTER_H
  3. /// \ingroup rec
  4. /// \class MpdTpcKalmanFilter
  5. /// \brief Kalman filter track reconstructor in MPD central detector
  6. ///
  7. /// \author Alexander Zinchenko, LHEP JINR Dubna
  8. //#include "TpcPadPlane.h"
  9. #include "TpcPoint.h"
  10. #include "MpdTpcHit.h"
  11. #include "MpdTpcSectorGeo.h"
  12. #include "MpdTpcKalmanTrack.h"
  13. #include "MpdKalmanHit.h"
  14. #include "FairTask.h"
  15. #include "FairField.h"
  16. #include <TH1F.h>
  17. #include <TMatrixD.h>
  18. #include <TMatrixDSym.h>
  19. #include "TClonesArray.h"
  20. #include <map>
  21. //#include <tuple>
  22. #include <vector>
  23. class MpdTpcKalmanFilter : public FairTask
  24. {
  25. public:
  26. struct matrix4 {
  27. TMatrixD xextr, xfilt, jacob, wfilt;
  28. matrix4() {resize();}
  29. matrix4(TMatrixD a, TMatrixD b, TMatrixD c, TMatrixD d) {
  30. resize();
  31. xextr = a;
  32. xfilt = b;
  33. jacob = c;
  34. wfilt = d;
  35. }
  36. void resize() {
  37. xextr.ResizeTo(5,1);
  38. xfilt.ResizeTo(5,1);
  39. jacob.ResizeTo(5,5);
  40. wfilt.ResizeTo(5,5);
  41. }
  42. };
  43. MpdTpcKalmanFilter(); ///< Default ctor
  44. MpdTpcKalmanFilter(const char *name, const char *title="TPC Kalman filter"); ///< Ctor
  45. virtual ~MpdTpcKalmanFilter(); ///< Destructor
  46. virtual InitStatus Init();
  47. virtual InitStatus ReInit();
  48. virtual void Exec(Option_t * option);
  49. virtual void Finish();
  50. void Reset();
  51. void Register();
  52. //virtual InitStatus ReInit();
  53. Int_t GetNofHitsInLayer(Int_t lay) { return (Int_t)fhLays->GetCellContent(lay+1,0); }
  54. Int_t GetHitsInLayer(Int_t lay) { return fLayPointers[lay]; } ///< first index of hits in layer
  55. TClonesArray* GetHits() { return fHits; } ///< get array of all hits
  56. TClonesArray* GetTracks() { return fTracks; } ///< get array of tracks
  57. Int_t RunKalmanFilter(MpdTpcKalmanTrack *track); ///< run Kalman filter
  58. Int_t GetParticleId(Int_t id); ///< particle ID for track id
  59. void SetModular(Int_t modular) { fModular = modular; } ///< set != 0 if modular geom. of r/out chambers
  60. void SetSectorGeo(MpdTpcSectorGeo *secGeo) { fSecGeo = secGeo; } ///< set sector geometry
  61. //Bool_t Refit(MpdKalmanTrack *track, Double_t mass = 0.13957, Int_t charge = 1, Bool_t skip = kFALSE, Int_t iDir = 1, Bool_t exclude = kFALSE, std::map<Double_t,std::tuple<TMatrixD,TMatrixD,TMatrixD,TMatrixD> > *cache = NULL); ///< refit track using its points for given particle mass and charge
  62. Bool_t Refit(MpdKalmanTrack *track, Double_t mass = 0.13957, Int_t charge = 1, Bool_t skip = kFALSE, Int_t iDir = 1, Bool_t exclude = kFALSE, std::map<Double_t,matrix4> *cache = NULL); ///< refit track using its points for given particle mass and charge
  63. Int_t GetHitID(MpdKalmanHit *hit); // get hit ID from MC point ID
  64. void UseTpcHit(Bool_t useMC = kTRUE) { fUseMCHit = useMC; }; // to use TpcHit branch instead of MpdTpcHit
  65. void FillGeoScheme(); // fill Kalman filter geometry manager info (for modular geometry of r/out chambers)
  66. Bool_t IsTpcHit() const { return fUseMCHit; } // tracking with hits or clusters ?
  67. void Smooth(MpdTpcKalmanTrack *track, std::vector<std::pair<Double_t,TMatrixD> >& vecSmooth);
  68. private:
  69. // Some constants
  70. static const Int_t fgkLays = 150; // number of padrows (with some margin)
  71. static const Int_t fgkSecs = 24; // number of sectors (with some margin)
  72. private:
  73. virtual void SetParContainers();
  74. void GetTrackSeeds(Int_t iPass); // build track seeds
  75. void GetTrackSeedsEndCaps(); // build track seeds from end-cap regions
  76. void DoTracking(Int_t iPass); // run tracking
  77. Double_t EvalPt(const MpdKalmanHit *hit1, const MpdKalmanHit *hit2, Double_t *params); // evaluate Pt
  78. void Cluster2KalmanHits(); // create Kalman hits from clusters
  79. void MakeKalmanHits(); // create Kalman hits
  80. void MakeKalmanHitsModul(); // create Kalman hits for modular geom. of r/out chambers
  81. void GoOutward(MpdTpcKalmanTrack *track); // propagate track outward
  82. Bool_t BackTrace(MpdTpcKalmanTrack *track, TMatrixDSym &weight, Int_t iDir = 1, Bool_t corDedx = kFALSE); // propagate track thru found hits
  83. void GoOut(); // backpropagate tracks outward
  84. void RemoveDoubles(TClonesArray *tracks); // remove double tracks
  85. Int_t GetNofCommonHits(MpdKalmanTrack *tr1, MpdKalmanTrack *tr2); // get common hits
  86. // Are compared tracks similar each other (if common hits in 2 tracks is greater limit)
  87. Bool_t AreTracksDoubles(MpdKalmanTrack *tr1, MpdKalmanTrack *tr2);
  88. Double_t CorrectForLoss(Double_t pt, Double_t the, Int_t id); // correct for dE loss in pipe
  89. //Bool_t SameOrigin(TpcLheHit *hit, Int_t idKF, Int_t *mcTracks); // check hit origin
  90. Bool_t SameOrigin(TpcPoint *hit, Int_t idKF, Int_t *mcTracks); // check hit origin
  91. void StoreTracks(); // transfer tracks from fTrackCand to fTracks
  92. void ExcludeHits(); // exclude used hits
  93. void RemoveShorts(); // remove short tracks (Nhits < 4)
  94. void GoToBeamLine(); // propagate tracks to the beam line
  95. void AddHits(Int_t indx0 = -1); // add hit objects to tracks
  96. void SetTrackID(MpdTpcKalmanTrack *track); // set track ID from IDs of its hits
  97. //TpcPoint* GetPoint(MpdKalmanHit *hit); // get MCPoint pointer for the hit
  98. MpdTpcHit* GetTpcHit(MpdKalmanHit *hit); // get TpcHit pointer for the Kalman hit
  99. Int_t SectorNo(const char* cpath); // extract sector number from TGeo path
  100. Double_t CorrectForLoss(Double_t pt, Double_t the, Double_t mass, Int_t charge = 1); ///< energy loss correction
  101. Double_t CorrectForLossFluct(Double_t pt, Double_t the, Double_t mass, Int_t charge = 1); ///< energy loss fluct. correction
  102. void MergeTracks(Int_t ipass); // merge tracks
  103. Double_t Interp2d(const Double_t *moms, const Double_t *thes, std::vector<std::vector<Double_t> > &sigmas,
  104. Double_t p, Double_t dip); // 2-d linear interpolation
  105. Int_t fNofEvents; // number of events processed
  106. Int_t fNTracks; // number of found tracks
  107. Int_t fNPass; // number of reco passes to run
  108. TClonesArray *fHits; // array of hits
  109. TClonesArray *fKHits; // array of Kalman hits
  110. TClonesArray *fTracks; // array of tracks
  111. TClonesArray *fTrackCand; // array of track candidates
  112. TClonesArray *fLHEtracks; // array of "LHE" MC tracks
  113. TClonesArray *fMCtracks; // array of MC tracks
  114. TClonesArray *fTpcPoints; // array of TPC points
  115. Int_t *fLayPointers; //! locations of hits from different layers
  116. Int_t fLaySecBegPointers[fgkLays][fgkSecs]; //! locations of hits from different layers and sectors (first hit)
  117. Int_t fLaySecEndPointers[fgkLays][fgkSecs]; //! locations of hits from different layers and sectors (last hit)
  118. TH1F *fhLays; // histo with layer hit multiplicities
  119. //std::multiset<Int_t> fLayset; // layer set of hits
  120. Double_t fVertZ; // primary vertex position estimate
  121. Int_t fZflag; // primary vertex position estimate quality flag
  122. Int_t fModular; // not equal 0 if modular geometry of readout chambers
  123. //const TpcPadPlane *fPadPlane; //! pointer to pad plane
  124. MpdTpcSectorGeo *fSecGeo; //! pointer to sector geometry
  125. Bool_t fUseMCHit; // to use TpcHit branch (hit producer) instead of MpdTpc (clusters)
  126. std::map<Double_t,matrix4> *fCache; // cached track parameters for smoother
  127. private:
  128. // Some constants
  129. static const Double_t fgkChi2Cut; // max accepted Chi2 of hit for track
  130. ClassDef(MpdTpcKalmanFilter,0);
  131. };
  132. #endif