MpdKalmanFilter.h 6.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117
  1. #ifndef MPDKALMANFILTER_H
  2. #define MPDKALMANFILTER_H
  3. /// \ingroup rec
  4. /// \class MpdKalmanFilter
  5. /// \brief Kalman filter track reconstructor in the MPD detector
  6. ///
  7. /// \author Alexander Zinchenko, LHEP JINR Dubna
  8. #include "MpdKalmanHit.h"
  9. #include "MpdKalmanTrack.h"
  10. #include "MpdKalmanGeoScheme.h"
  11. #include "FairTask.h"
  12. #include "FairField.h"
  13. #include <TMatrixD.h>
  14. #include <TMatrixDSym.h>
  15. #include "TClonesArray.h"
  16. class MpdKalmanFilter : public FairTask
  17. {
  18. public:
  19. // KG public cpnstructor for TStreamer reading
  20. MpdKalmanFilter(const char *name="MpdKalmanFilter", const char *title="MPD Task"); ///< Ctor
  21. static MpdKalmanFilter* Instance(); ///< get singleton instance
  22. static MpdKalmanFilter* Instance(const char *name, const char *title="MPD Task"); ///< get singleton instance
  23. virtual void Exec(Option_t * option);
  24. void Reset();
  25. void Register();
  26. MpdKalmanGeoScheme* GetGeo() { return fGeoScheme; } ///< pointer to geo manager
  27. Bool_t PropagateToHit(MpdKalmanTrack *track, const MpdKalmanHit *hit, Bool_t calcLeng = kTRUE, Bool_t local = kFALSE, Double_t stepBack = -1.0); ///< propagate track
  28. Bool_t PropagateParamToHit(MpdKalmanTrack *track, const MpdKalmanHit *hit, Bool_t calcLeng = kTRUE, Bool_t local = kFALSE, Double_t stepBack = -1.0); ///< propagate track parameters
  29. Double_t FilterHit(MpdKalmanTrack *track, const MpdKalmanHit *hit,
  30. TMatrixDSym &pointWeight, TMatrixD &paramTmp); ///< compute Chi2 from hit addition (filter hit)
  31. Double_t FilterHitR(MpdKalmanTrack *track, const MpdKalmanHit *hit,
  32. TMatrixDSym &pointWeight, TMatrixD &paramTmp); ///< compute Chi2 from hit addition for "barrel" tracks and hits
  33. Double_t FilterHitZ(MpdKalmanTrack *track, const MpdKalmanHit *hit,
  34. TMatrixDSym &pointWeight, TMatrixD &parFilt); ///< compute Chi2 from hit addition for "end-cap" tracks and hits
  35. Double_t FilterStrip(MpdKalmanTrack *track, const MpdKalmanHit *hit, TMatrixDSym &pointWeight,
  36. TMatrixD &parFilt); ///< Kalman filtering for "strip" hits
  37. Double_t FilterStripLocal(MpdKalmanTrack *track, const MpdKalmanHit *hit, TMatrixDSym &pointWeight,
  38. TMatrixD &parFilt, Double_t &posNew); ///< Kalman filtering for "strip" hits in local coordinates
  39. void MnvertLocal(Double_t *a, Int_t l, Int_t, Int_t n, Int_t &ifail);
  40. Double_t Proxim(Double_t phi0, Double_t phi, Double_t scale = 1.0); // adjust angle for continuity
  41. Double_t Proxim(MpdKalmanTrack *track, const MpdKalmanHit *hit, Double_t scale = 1.0); ///< adjust R-Phi coord. for continuity
  42. Bool_t PropagateParamR(MpdKalmanTrack *track, const MpdKalmanHit *hit, Bool_t calcLeng, Double_t *vert = (Double_t*)0x0); ///< propagate params
  43. Bool_t PropagateParamP(MpdKalmanTrack *track, const MpdKalmanHit *hit, Bool_t calcLeng, Bool_t local = kFALSE, Double_t stepBack = -1.0); ///< propagate params
  44. Bool_t PropagateParamP(MpdKalmanTrack *track, const Double_t *plane, Bool_t calcLeng, Bool_t local = kFALSE, Double_t stepBack = -1.0); ///< propagate params to plane
  45. Bool_t PropagateParamZ(MpdKalmanTrack *track, const MpdKalmanHit *hit, Bool_t calcLeng); // propagate params
  46. Bool_t FindPca(MpdKalmanTrack *track, Double_t *vert); ///< find PCA w.r.t. given point
  47. void Convert(const MpdKalmanTrack *track1, MpdKalmanTrack *track2); ///< convert track of diff. types
  48. Bool_t Refit(MpdKalmanTrack *track, Int_t iDir, Bool_t leng=kFALSE); ///< refit track using its points
  49. Double_t Interp(Int_t nDim, const Double_t *mom, const Double_t *dedxm, Double_t p); ///< energy loss interpolation
  50. Double_t Scattering(MpdKalmanTrack *track, Double_t dx, TString mass2="0.0194797849", Int_t charge=1); // compute mult. scat. angle
  51. Double_t Scattering(MpdKalmanTrack *track, Double_t x0, Double_t step, TString mass2="0.0194797849", Int_t charge=1); // compute mult. scat. angle
  52. void SetGeantParamB(MpdKalmanTrack *track, Double_t *v3, Double_t dir); // fill Geant vector
  53. void GetField(Double_t *position, Double_t *field); // get mag. field vector
  54. Int_t IsNumer() const { return fNumer; } // if != 0 - numerical propagator
  55. void SetNumer(Int_t numer) { fNumer = numer; } // if != 0 - numerical propagator
  56. Double_t* ExtrapOneStep(MpdKalmanTrack *track, Double_t step, Int_t flag = 0); // propagate thru one step (use cache if flag != 0)
  57. virtual InitStatus Init();
  58. TMatrixD& GetJacob() { return fJacob; } // return Jacobian
  59. protected:
  60. virtual InitStatus ReInit();
  61. virtual void Finish();
  62. virtual ~MpdKalmanFilter(); ///< Destructor
  63. private:
  64. // automatic deleting when application exit
  65. static void DestroyInstance (){
  66. if (fgKF)
  67. delete fgKF;
  68. }
  69. static MpdKalmanFilter* fgKF; //! pointer to Singleton instance, may be ! required to exclude recursion
  70. Bool_t AnalParamZ(MpdKalmanTrack *track, const MpdKalmanHit *hit, Bool_t calcLeng); // propagate params analytically
  71. Bool_t AnalParamX(MpdKalmanTrack *track, const MpdKalmanHit *hit, Bool_t calcLeng, Bool_t local = kFALSE); // propagate params analytically
  72. Bool_t PropagateWeight(MpdKalmanTrack *track, const MpdKalmanHit *hit, Double_t sign = 1., Double_t stepBack = -1.0); // propagate weight
  73. Bool_t AnalyticJacob(MpdKalmanTrack *track, const MpdKalmanHit *hit, TMatrixD &jacob); // analytical Jacobian computation
  74. Bool_t AnalyticJacobX(MpdKalmanTrack *track, const MpdKalmanHit *hit, TMatrixD &jacob); // analytical Jacobian computation
  75. void SetGeantParamE(MpdKalmanTrack *track, Double_t *v3, Double_t dir); // fill Geant vector
  76. Double_t Proxim(const MpdKalmanHit *hit0, const MpdKalmanHit *hit, Double_t scale = 1.0); // adjust R-Phi coord. for continuity
  77. void GoOutward(MpdKalmanTrack *track); // propagate track outward
  78. void BackTrace(MpdKalmanTrack *track, TMatrixDSym &weight); // propagate track thru found hits
  79. void GetFromGeantParamB(MpdKalmanTrack *track, Double_t *v3, Double_t dir); // fill params
  80. void GetFromGeantParamE(MpdKalmanTrack *track, Double_t *v3, Double_t dir); // fill params
  81. void ExtrapOneStepRungekutta(Double_t charge, Double_t step, Double_t* vect,
  82. Double_t* vout); // RungeKutta track propagation
  83. void ExtrapOneStepHelix(Double_t charge, Double_t step, Double_t* vect,
  84. Double_t* vout); // helix track propagation
  85. void ExtrapOneStepHelix3(Double_t charge, Double_t step, Double_t* vect,
  86. Double_t* vout) const; // helix track propagation
  87. void SetField(FairField *field); // set mag. field
  88. // KG exlude to copy by TStreamer because of get in Init()
  89. FairField *fMagField; //! magnetic field
  90. private:
  91. MpdKalmanGeoScheme *fGeoScheme; // pointer to geometry manager
  92. Int_t fNumer; // if != 0 - numerical propagator
  93. Double_t fVectorG[8]; // cache for propagator
  94. TMatrixD fJacob;
  95. // Some constants
  96. static const Int_t fgkTriesMax; // max number of attempts to find exact position
  97. static const Double_t fgkEpsilon; // tracking precision (cm)
  98. ClassDef(MpdKalmanFilter,0);
  99. };
  100. #endif