MpdEctKalmanTrack.cxx 9.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306
  1. /// \class MpdEctKalmanTrack
  2. ///
  3. /// Kalman filter track object for the MPD straw end-cap detector
  4. /// Track parameters: 0: X - X-coordinate @ fixed Z
  5. /// 1: Y - Y-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 (LHE, JINR, Dubna)
  10. #include "MpdEctKalmanTrack.h"
  11. #include "FairHit.h"
  12. #include "MpdKalmanFilter.h"
  13. #include "MpdKalmanTrack.h"
  14. //#include "MpdKalmanHitZ.h"
  15. #include "MpdKalmanHit.h"
  16. //#include "TpcLheKalmanTrack.h"
  17. #include "MpdTpcKalmanTrack.h"
  18. #include <TMatrixD.h>
  19. #include <TMatrixDSym.h>
  20. #include <TClonesArray.h>
  21. #include <TMath.h>
  22. #include <Riostream.h>
  23. using namespace std;
  24. //__________________________________________________________________________
  25. MpdEctKalmanTrack::MpdEctKalmanTrack()
  26. : MpdKalmanTrack(),
  27. fLastLay(0),
  28. fTpcIndex(-1),
  29. fTofIndex(-1),
  30. fFlag(0),
  31. fMirrors(0),
  32. fMisses(0),
  33. fTrHits(0x0)
  34. {
  35. /// Default constructor
  36. //SetType(MpdKalmanTrack::kFixedZ);
  37. SetType(MpdKalmanTrack::kEndcap);
  38. }
  39. //__________________________________________________________________________
  40. MpdEctKalmanTrack::MpdEctKalmanTrack(Int_t tpcIndx, const MpdTpcKalmanTrack &track)
  41. : MpdKalmanTrack(track),
  42. fLastLay(0),
  43. fTpcIndex(tpcIndx),
  44. fTofIndex(-1),
  45. fFlag(0),
  46. fMirrors(0),
  47. fMisses(0),
  48. fTrHits(new TClonesArray("MpdKalmanHit"))
  49. {
  50. /// Constructor from TPC track
  51. //MpdKalmanFilter::Instance()->Convert(track, this);
  52. SetDirection(MpdKalmanTrack::kOutward);
  53. /*
  54. SetWeight(*fCovar);
  55. Int_t iok;
  56. MpdKalmanFilter::Instance()->MnvertLocal(fWeight->GetMatrixArray(), 5, 5, 5, iok);
  57. fPos = fPosNew;
  58. *fParamNew = *fParam;
  59. */
  60. SetWeight(*fWeightAtHit);
  61. fPosNew = fPos = fPosAtHit;
  62. *fParamNew = *fParamAtHit;
  63. *fParam = *fParamNew;
  64. fLength = fLengAtHit;
  65. //fHits->Clear();
  66. fNhits = 0;
  67. Int_t nHits = track.GetNofTrHits();
  68. for (Int_t i = 0; i < nHits; ++i) fHits->Add(track.GetTrHits()->UncheckedAt(i));
  69. SetFromTpc();
  70. }
  71. //__________________________________________________________________________
  72. MpdEctKalmanTrack::MpdEctKalmanTrack(Int_t tofIndx, Int_t tpcIndx,
  73. FairHit *tof, TpcLheHit *tpc, TVector3 &vert)
  74. : MpdKalmanTrack(tof->GetZ(),vert),
  75. fLastLay(0),
  76. fTpcIndex(tpcIndx),
  77. fTofIndex(tofIndx),
  78. fFlag(0),
  79. fMirrors(0),
  80. fMisses(0),
  81. //fTrHits(new TClonesArray("MpdKalmanHitZ"))
  82. fTrHits(new TClonesArray("MpdKalmanHit"))
  83. {
  84. /// Constructor from ETOF and TPC hits
  85. //SetTrackID(tof->GetTrackID());
  86. SetTrackID(-1);
  87. SetType(MpdKalmanTrack::kEndcap);
  88. //SetDirection(MpdKalmanTrack::kInward);
  89. }
  90. //__________________________________________________________________________
  91. MpdEctKalmanTrack::MpdEctKalmanTrack(Int_t tofIndx, Int_t ectIndx,
  92. //MpdEtofPoint *tof, MpdKalmanHitZ *ect, TVector3 &vert)
  93. FairHit *tof, MpdKalmanHit *ect, TVector3 &vert)
  94. : MpdKalmanTrack(tof->GetZ(),vert),
  95. fLastLay(-1),
  96. fTpcIndex(ectIndx),
  97. fTofIndex(tofIndx),
  98. fFlag(0),
  99. fMirrors(0),
  100. fMisses(0),
  101. fTrHits(new TClonesArray("MpdKalmanHit"))
  102. {
  103. /// Constructor from ETOF and ECT (CPC) hits
  104. // Add first hit
  105. if (ect) {
  106. fLastLay = ect->GetLayer();
  107. fHits->Add(ect);
  108. }
  109. SetType(MpdKalmanTrack::kEndcap);
  110. SetDirection(MpdKalmanTrack::kInward);
  111. SetFromEtof();
  112. }
  113. //__________________________________________________________________________
  114. MpdEctKalmanTrack::MpdEctKalmanTrack (const MpdEctKalmanTrack& track)
  115. : MpdKalmanTrack(track),
  116. fLastLay(track.fLastLay),
  117. fTpcIndex(track.fTpcIndex),
  118. fTofIndex(track.fTofIndex),
  119. fFlag(track.fFlag),
  120. fMirrors(track.fMirrors),
  121. fMisses(track.fMisses),
  122. //fTrHits(new TClonesArray("MpdKalmanHitZ"))
  123. fTrHits(new TClonesArray("MpdKalmanHit"))
  124. {
  125. ///copy constructor
  126. Int_t nHits = track.fTrHits->GetEntriesFast();
  127. for (Int_t i = 0; i < nHits; ++i) {
  128. //MpdKalmanHitZ *hit = (MpdKalmanHitZ*)(track.fTrHits->UncheckedAt(i));
  129. //new ((*fTrHits)[i]) MpdKalmanHitZ(*hit);
  130. MpdKalmanHit *hit = (MpdKalmanHit*)(track.fTrHits->UncheckedAt(i));
  131. new ((*fTrHits)[i]) MpdKalmanHit(*hit);
  132. }
  133. }
  134. //__________________________________________________________________________
  135. MpdEctKalmanTrack & MpdEctKalmanTrack::operator=(const MpdEctKalmanTrack& track)
  136. {
  137. /// Asignment operator
  138. // check assignement to self
  139. if (this == &track) return *this;
  140. // base class assignement
  141. MpdKalmanTrack::operator=(track);
  142. fLastLay = track.fLastLay;
  143. fTpcIndex = track.fTpcIndex;
  144. fTofIndex = track.fTofIndex;
  145. fFlag = track.fFlag;
  146. fMirrors = track.fMirrors;
  147. fMisses = track.fMisses;
  148. //fTrHits = new TClonesArray("MpdKalmanHitZ");
  149. fTrHits = new TClonesArray("MpdKalmanHit");
  150. //fLengAtHit = track.fLengAtHit;
  151. Int_t nHits = track.fTrHits->GetEntriesFast();
  152. for (Int_t i = 0; i < nHits; ++i) {
  153. //MpdKalmanHitZ *hit = (MpdKalmanHitZ*)(track.fTrHits->UncheckedAt(i));
  154. //new ((*fTrHits)[i]) MpdKalmanHitZ(*hit);
  155. MpdKalmanHit *hit = (MpdKalmanHit*)(track.fTrHits->UncheckedAt(i));
  156. new ((*fTrHits)[i]) MpdKalmanHit(*hit);
  157. }
  158. return *this;
  159. }
  160. //__________________________________________________________________________
  161. MpdEctKalmanTrack::~MpdEctKalmanTrack()
  162. {
  163. /// Destructor
  164. //if (fTrHits) fTrHits->Clear("C");
  165. if (fTrHits) fTrHits->Delete();
  166. delete fTrHits;
  167. }
  168. //__________________________________________________________________________
  169. void MpdEctKalmanTrack::Reset()
  170. {
  171. /// Reset track
  172. //if (fTrHits) fTrHits->Clear("C");
  173. if (fTrHits) fTrHits->Delete();
  174. delete fTrHits;
  175. Clear();
  176. }
  177. //__________________________________________________________________________
  178. Double_t MpdEctKalmanTrack::GetRadNew()
  179. {
  180. /// Get track new radial position
  181. //if (fTrackType == kFixedZ) {
  182. if (fTrackType == kEndcap) {
  183. Double_t x = GetParamNew(0), y = GetParamNew(1);
  184. Double_t r = x * x + y * y;
  185. return TMath::Sqrt(r);
  186. }
  187. return fPosNew;
  188. }
  189. //__________________________________________________________________________
  190. /*
  191. void TpcLheKalmanTrack::EvalCovar(const MpdKalmanHitR *hitOut, const MpdKalmanHitR *hitIn)
  192. {
  193. /// Evaluate covariance matrix for track seed
  194. (*fWeight)(0,0) = hitOut->GetRphiErr() * hitOut->GetRphiErr(); // <RphiRphi>
  195. (*fWeight)(0,0) *= 4.; // extra factor of 2
  196. (*fWeight)(1,1) = hitOut->GetZerr() * hitOut->GetZerr(); // <zz>
  197. Double_t phiOut = hitOut->GetRphi() / hitOut->GetR();
  198. Double_t phiIn = hitIn->GetRphi() / hitIn->GetR();
  199. Double_t xIn = hitIn->GetR() * TMath::Cos(phiIn);
  200. Double_t yIn = hitIn->GetR() * TMath::Sin(phiIn);
  201. Double_t xOut = hitOut->GetR() * TMath::Cos(phiOut);
  202. Double_t yOut = hitOut->GetR() * TMath::Sin(phiOut);
  203. Double_t dist2 = (xOut-xIn)*(xOut-xIn) + (yOut-yIn)*(yOut-yIn);
  204. Double_t sinPhi = TMath::Sin ((*fParam)(2,0));
  205. Double_t cosPhi = TMath::Cos ((*fParam)(2,0));
  206. Double_t pOut = TMath::Cos(phiOut) * cosPhi + TMath::Sin(phiOut) * sinPhi;
  207. Double_t pIn = TMath::Cos(phiIn) * cosPhi + TMath::Sin(phiIn) * sinPhi;
  208. (*fWeight)(2,2) = (pOut * pOut + pIn * pIn) / dist2 * (*fWeight)(0,0); // <PhiPhi>
  209. (*fWeight)(2,2) *= 4.; // extra factor of 2
  210. Double_t tanThe = TMath::Tan((*fParam)(3,0));
  211. Double_t dRad = hitOut->GetR() - hitIn->GetR();
  212. Double_t denom = dRad * (1.+tanThe*tanThe);
  213. (*fWeight)(3,3) = (*fWeight)(1,1) * 2. / denom / denom; // <TheThe>
  214. //(*fWeight)(4,4) = ((*fParam)(4,0)*0.5) * ((*fParam)(4,0)*0.5); // error 50%
  215. //(*fWeight)(4,4) = ((*fParam)(4,0)*0.75) * ((*fParam)(4,0)*0.75); // error 75%
  216. (*fWeight)(4,4) = ((*fParam)(4,0)*2.) * ((*fParam)(4,0)*2.); // error 200%
  217. //fWeight->Print();
  218. //fWeight->Invert(); // weight matrix
  219. Int_t iok = 0;
  220. MpdKalmanFilter::Instance()->MnvertLocal(fWeight->GetMatrixArray(), 5, 5, 5, iok);
  221. //fWeight->Print();
  222. }
  223. */
  224. //__________________________________________________________________________
  225. void MpdEctKalmanTrack::StartBack()
  226. {
  227. /// Prepare for back tracing
  228. fHits->Sort(); // in descending order in radius
  229. fChi2 = 0.;
  230. fTrackDir = kInward;
  231. //if (fLastLay == 50) { cout << " back" << endl; Weight2Cov()->Print(); exit(0); }
  232. if (GetNofHits() == 1) return;
  233. /*
  234. for (Int_t i = 0; i < 5; ++i) {
  235. for (Int_t j = i; j < 5; ++j) {
  236. if (j == i) (*fWeight)(i,j) /= 100.;
  237. else (*fWeight)(i,j) = (*fWeight)(j,i) = 0.;
  238. }
  239. }
  240. */
  241. fWeight->Zero();
  242. //EvalCovar((MpdKalmanHitR*)fHits->UncheckedAt(0),(MpdKalmanHitR*)fHits->UncheckedAt(1));
  243. }
  244. //__________________________________________________________________________
  245. void MpdEctKalmanTrack::Print(Option_t *opt)
  246. {
  247. /// Print track info
  248. cout << " -------------------------------------------------------------------- " << endl;
  249. cout << " Track ID: " << GetTrackID() << "; type: " << GetType() << "; direction: " << GetDirection() << endl;
  250. cout << " Track position: ";
  251. if (GetType() == MpdKalmanTrack::kBarrel)
  252. cout << "(R-Phi, Z, R): " << GetParam(0) << " " << GetParam(1) << " " << GetPos() << endl;
  253. else cout << "(X, Y, Z): " << GetParam(0) << " " << GetParam(1) << " " << GetPos() << endl;
  254. cout << " Track Pt: " << 1. / GetParam(4) << endl;
  255. cout << " Track new position ";
  256. if (fParamNew && GetType() == MpdKalmanTrack::kBarrel)
  257. cout << "(R-Phi, Z, R): " << GetParamNew(0) << " " << GetParamNew(1) << " " << GetPosNew() << endl;
  258. else cout << "(X, Y, Z): " << GetParamNew(0) << " " << GetParamNew(1) << " " << GetPosNew() << endl;
  259. cout << " Track Pt: " << 1. / GetParamNew(4) << endl;
  260. cout << " Number of hits: " << GetHits()->GetEntriesFast() << " " << GetNofTrHits()
  261. << ", Chi2: " << GetChi2() << endl;
  262. cout << " -------------------------------------------------------------------- " << endl;
  263. }
  264. ClassImp(MpdEctKalmanTrack)