MpdItsKalmanTrack.cxx 9.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289
  1. /// \class MpdItsKalmanTrack
  2. ///
  3. /// Kalman filter track object for the MPD central detector
  4. /// Track parameters: 0: RPhi - coordinate in R-Phi direction
  5. /// 1: Z - longitudinal coordinate
  6. /// 2: Phi - local azimuthal angle
  7. /// 3: Theta - local dip angle (angle w.r.t. the transverse plane)
  8. /// 4: q/Pt - signed inverse Pt
  9. /// \author Alexander Zinchenko (LHEP, JINR, Dubna)
  10. #include "MpdItsKalmanTrack.h"
  11. #include "MpdTpcKalmanTrack.h"
  12. #include "MpdEctKalmanTrack.h"
  13. #include "MpdKalmanFilter.h"
  14. #include "MpdKalmanGeoScheme.h"
  15. #include "MpdKalmanHit.h"
  16. #include "MpdKalmanTrack.h"
  17. ///#include "MpdCellAutomat.h"
  18. #include "MpdVector.h"
  19. ///#include "MpdCellTrack.h"
  20. #include <TMatrixD.h>
  21. #include <TMatrixDSym.h>
  22. #include <TClonesArray.h>
  23. #include <TVector3.h>
  24. #include <TMath.h>
  25. #include <Riostream.h>
  26. //__________________________________________________________________________
  27. MpdItsKalmanTrack::MpdItsKalmanTrack()
  28. : MpdKalmanTrack(),
  29. fTrHits(0x0),
  30. // fKHits(0x0),//attention!
  31. fNofIts(0),
  32. fChi2Its(0.0)
  33. {
  34. /// Default constructor
  35. }
  36. //__________________________________________________________________________
  37. MpdItsKalmanTrack::MpdItsKalmanTrack(MpdKalmanHit *hitOut, MpdKalmanHit *hitIn,
  38. TVector3 &vertex, Double_t pt)
  39. : MpdKalmanTrack(0., vertex),
  40. fTrHits(new TClonesArray("MpdKalmanHit")),
  41. // fKHits(new TClonesArray("MpdKalmanHit")), //attention!!!
  42. fNofIts(0),
  43. fChi2Its(0.0)
  44. {
  45. /// Constructor from 2 hits
  46. Double_t rOut = 0, phiOut = 0, rIn = 0, phiIn = 0;
  47. if (hitIn->GetType() == MpdKalmanHit::kFixedP) {
  48. TVector3 posIn = MpdKalmanFilter::Instance()->GetGeo()->GlobalPos(hitIn);
  49. TVector3 posOut = MpdKalmanFilter::Instance()->GetGeo()->GlobalPos(hitOut);
  50. rOut = posOut.Pt();
  51. phiOut = posOut.Phi();
  52. rIn = posIn.Pt();
  53. phiIn = posIn.Phi();
  54. } else {
  55. rOut = hitOut->GetPos();
  56. phiOut = hitOut->GetMeas(0) / rOut;
  57. rIn = hitIn->GetPos();
  58. phiIn = hitIn->GetMeas(0) / rIn;
  59. }
  60. Double_t xIn = rIn * TMath::Cos(phiIn);
  61. Double_t yIn = rIn * TMath::Sin(phiIn);
  62. Double_t xOut = rOut * TMath::Cos(phiOut);
  63. Double_t yOut = rOut * TMath::Sin(phiOut);
  64. Double_t parOut[4] = {rOut, phiOut, xOut, yOut};
  65. Double_t parIn[4] = {rIn, phiIn, xIn, yIn};
  66. fPos = fPosNew = rOut;
  67. (*fParam)(0,0) = rOut * phiOut; // R-Phi coordinate
  68. //(*fParam)(1,0) = posOut.Z(); // Z-coordinate
  69. (*fParam)(1,0) = hitOut->GetMeas(1); // Z-coordinate
  70. (*fParam)(2,0) = TMath::ATan2 (yOut-yIn, xOut-xIn); // Track Phi
  71. //(*fParam)(3,0) = TMath::ATan2 (posOut.Z()-posIn.Z(), rOut-rIn); // Track Theta
  72. (*fParam)(3,0) = TMath::ATan2 (hitOut->GetMeas(1)-hitIn->GetMeas(1), rOut-rIn); // Track Theta
  73. (*fParam)(4,0) = 1. / pt; // q/Pt
  74. *fParamNew = *fParam;
  75. fWeight->Zero();
  76. EvalCovar(hitOut, hitIn, parOut, parIn);
  77. fHits->Add(hitOut);
  78. fLastLay = hitOut->GetLayer(); // layer number
  79. }
  80. //__________________________________________________________________________
  81. MpdItsKalmanTrack::MpdItsKalmanTrack (const MpdItsKalmanTrack& track)
  82. : MpdKalmanTrack(track),
  83. fTrHits(new TClonesArray("MpdKalmanHit")),
  84. fNofIts(track.fNofIts),
  85. fChi2Its(track.fChi2Its)
  86. {
  87. ///copy constructor
  88. Int_t nHits = track.fTrHits->GetEntriesFast();
  89. for (Int_t i = 0; i < nHits; ++i) {
  90. MpdKalmanHit *hit = (MpdKalmanHit*) (track.fTrHits->UncheckedAt(i));
  91. new ((*fTrHits)[i]) MpdKalmanHit(*hit);
  92. }
  93. }
  94. //__________________________________________________________________________
  95. MpdItsKalmanTrack & MpdItsKalmanTrack::operator=(const MpdItsKalmanTrack& track)
  96. {
  97. /// Asignment operator
  98. // check assignement to self
  99. if (this == &track) return *this;
  100. // base class assignement
  101. MpdKalmanTrack::operator=(track);
  102. if (fTrHits) { fTrHits->Delete(); delete fTrHits; }
  103. fTrHits = new TClonesArray("MpdKalmanHit");
  104. fNofIts = track.fNofIts;
  105. fChi2Its = track.fChi2Its;
  106. Int_t nHits = track.fTrHits->GetEntriesFast();
  107. for (Int_t i = 0; i < nHits; ++i) {
  108. MpdKalmanHit *hit = (MpdKalmanHit*)(track.fTrHits->UncheckedAt(i));
  109. new ((*fTrHits)[i]) MpdKalmanHit(*hit);
  110. }
  111. if (track.fHits == 0x0) {
  112. for (Int_t i = 0; i < nHits; ++i) {
  113. fHits->Add(fTrHits->UncheckedAt(i));
  114. }
  115. }
  116. return *this;
  117. }
  118. //__________________________________________________________________________
  119. MpdItsKalmanTrack::MpdItsKalmanTrack (const MpdTpcKalmanTrack& track)
  120. : MpdKalmanTrack(track),
  121. fTrHits(new TClonesArray("MpdKalmanHit")),
  122. fNofIts(0),
  123. fChi2Its(0.0)
  124. {
  125. /// constructor from TPC track
  126. TClonesArray* hits = track.GetTrHits();
  127. Int_t nHits = hits->GetEntriesFast();
  128. for (Int_t i = 0; i < nHits; ++i) {
  129. MpdKalmanHit *hit = (MpdKalmanHit*) (hits->UncheckedAt(i));
  130. new ((*fTrHits)[i]) MpdKalmanHit(*hit);
  131. }
  132. }
  133. //__________________________________________________________________________
  134. MpdItsKalmanTrack::MpdItsKalmanTrack (const MpdEctKalmanTrack& track)
  135. : MpdKalmanTrack(track),
  136. fTrHits(new TClonesArray("MpdKalmanHit")),
  137. fNofIts(0),
  138. fChi2Its(0.0)
  139. {
  140. /// constructor from TPC track
  141. TClonesArray* hits = track.GetTrHits();
  142. Int_t nHits = hits->GetEntriesFast();
  143. for (Int_t i = 0; i < nHits; ++i) {
  144. MpdKalmanHit *hit = (MpdKalmanHit*) (hits->UncheckedAt(i));
  145. new ((*fTrHits)[i]) MpdKalmanHit(*hit);
  146. }
  147. }
  148. //__________________________________________________________________________
  149. MpdItsKalmanTrack::MpdItsKalmanTrack (const MpdVector& track, TVector3& vec)
  150. : MpdKalmanTrack(0.0, vec),
  151. fTrHits(new TClonesArray("MpdKalmanHit")),
  152. fNofIts(0),
  153. fChi2Its(0.0)
  154. {
  155. /// construct from last layer celltrack
  156. const MpdVector* temp = &track;
  157. for (Int_t i = 0; i < 5; ++i) {
  158. MpdKalmanHit* hit = temp->GetKalmanHit();
  159. temp = temp->GetPrevTrackPointer();
  160. new ((*fTrHits)[i]) MpdKalmanHit(*hit);
  161. if (temp == NULL) break;
  162. }
  163. fTrHits->Sort();
  164. }
  165. //__________________________________________________________________________
  166. MpdItsKalmanTrack::~MpdItsKalmanTrack()
  167. {
  168. /// Destructor
  169. ///if (fTrHits) fTrHits->Clear("C");
  170. if (fTrHits) fTrHits->Delete(); //AZ
  171. delete fTrHits;
  172. fTrHits = NULL;
  173. }
  174. //__________________________________________________________________________
  175. void MpdItsKalmanTrack::Reset()
  176. {
  177. /// Reset track
  178. //if (fTrHits) fTrHits->Clear("C");
  179. if (fTrHits) fTrHits->Delete();
  180. delete fTrHits;
  181. fTrHits = NULL;
  182. Clear();
  183. }
  184. //__________________________________________________________________________
  185. void MpdItsKalmanTrack::EvalCovar(const MpdKalmanHit *hitOut, const MpdKalmanHit *hitIn, Double_t *parOut, Double_t *parIn)
  186. {
  187. /// Evaluate covariance matrix for track seed
  188. (*fWeight)(0,0) = hitOut->GetErr(0) * hitOut->GetErr(0); // <RphiRphi>
  189. (*fWeight)(0,0) *= 4.; // extra factor of 2
  190. (*fWeight)(1,1) = hitOut->GetErr(1) * hitOut->GetErr(1); // <zz>
  191. Double_t phiOut = parOut[1], phiIn = parIn[1];
  192. Double_t dx = parOut[2] - parIn[2], dy = parOut[3] - parIn[3];
  193. Double_t dist2 = dx * dx + dy * dy;
  194. Double_t sinPhi = TMath::Sin ((*fParam)(2,0));
  195. Double_t cosPhi = TMath::Cos ((*fParam)(2,0));
  196. Double_t pOut = TMath::Cos(phiOut) * cosPhi + TMath::Sin(phiOut) * sinPhi;
  197. Double_t pIn = TMath::Cos(phiIn) * cosPhi + TMath::Sin(phiIn) * sinPhi;
  198. (*fWeight)(2,2) = (pOut * pOut + pIn * pIn) / dist2 * (*fWeight)(0,0); // <PhiPhi>
  199. (*fWeight)(2,2) *= 4.; // extra factor of 2
  200. Double_t tanThe = TMath::Tan((*fParam)(3,0));
  201. Double_t dRad = parOut[0] - parIn[0];
  202. Double_t denom = dRad * (1.+tanThe*tanThe);
  203. (*fWeight)(3,3) = (*fWeight)(1,1) * 2. / denom / denom; // <TheThe>
  204. //(*fWeight)(3,3) = (*fWeight)(1,1) * 20. / denom / denom; // <TheThe>
  205. //(*fWeight)(4,4) = ((*fParam)(4,0)*0.5) * ((*fParam)(4,0)*0.5); // error 50%
  206. //(*fWeight)(4,4) = ((*fParam)(4,0)*0.75) * ((*fParam)(4,0)*0.75); // error 75%
  207. (*fWeight)(4,4) = ((*fParam)(4,0)*2.) * ((*fParam)(4,0)*2.); // error 200%
  208. //fWeight->Print();
  209. //fWeight->Invert(); // weight matrix
  210. Int_t iok = 0;
  211. MpdKalmanFilter::Instance()->MnvertLocal(fWeight->GetMatrixArray(), 5, 5, 5, iok);
  212. //fWeight->Print();
  213. }
  214. //__________________________________________________________________________
  215. void MpdItsKalmanTrack::StartBack()
  216. {
  217. /// Prepare for back tracing
  218. fHits->Sort(); // in descending order in radius
  219. fChi2 = 0.;
  220. fTrackDir = kInward;
  221. //if (fLastLay == 50) { cout << " back" << endl; Weight2Cov()->Print(); exit(0); }
  222. Int_t nHits = GetNofHits();
  223. if (nHits == 1) return;
  224. nHits *= nHits;
  225. //*
  226. for (Int_t i = 0; i < 5; ++i) {
  227. for (Int_t j = i; j < 5; ++j) {
  228. //if (j == i) (*fWeight)(i,j) /= 100.;
  229. if (j == i) (*fWeight)(i,j) /= nHits;
  230. else (*fWeight)(i,j) = (*fWeight)(j,i) = 0.;
  231. }
  232. }
  233. //*/
  234. //fWeight->Zero();
  235. //EvalCovar((MpdKalmanHitR*)fHits->UncheckedAt(0),(MpdKalmanHitR*)fHits->UncheckedAt(1));
  236. }
  237. //__________________________________________________________________________
  238. Int_t MpdItsKalmanTrack::Compare(const TObject* track) const
  239. {
  240. /// "Compare" function to sort in descending order in pt
  241. MpdKalmanTrack *trackKF = (MpdKalmanTrack*) track;
  242. Double_t pt = 1. / TMath::Abs(trackKF->GetParam(4));
  243. Double_t ptthis = 1. / TMath::Abs((*fParam)(4,0));
  244. if (ptthis < pt) return 1;
  245. else if (ptthis > pt) return -1;
  246. return 0;
  247. }
  248. ClassImp(MpdItsKalmanTrack)