MpdParticle.cxx 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476
  1. /// \class MpdParticle
  2. ///
  3. /// Particle object in MPD (to work with decays).
  4. ///
  5. /// \author Alexander Zinchenko (LHEP, JINR, Dubna)
  6. #include "MpdParticle.h"
  7. #include "MpdHelix.h"
  8. #include "MpdKalmanFilter.h"
  9. #include "MpdKalmanHit.h"
  10. #include "MpdKalmanTrack.h"
  11. #include "MpdMotherFitterPart.h"
  12. #include <TDatabasePDG.h>
  13. //__________________________________________________________________________
  14. MpdParticle::MpdParticle()
  15. : TObject(), fIndx(-1), fCharge(0),
  16. fMass(TDatabasePDG::Instance()->GetParticle("pi+")->Mass()), fChi2(0.0), fChi2ver(-1.0),
  17. fFlag(0), fPoint00(kTRUE)
  18. {
  19. /// Default constructor
  20. fMeas.ResizeTo(5,1);
  21. fq.ResizeTo(3,1);
  22. fx.ResizeTo(3,1);
  23. fJ.ResizeTo(3,3);
  24. fJinv.ResizeTo(3,3);
  25. fD.ResizeTo(3,3);
  26. fE.ResizeTo(3,3);
  27. //fQ.ResizeTo(3,3);
  28. fA.ResizeTo(5,3);
  29. fB.ResizeTo(5,3);
  30. fC.ResizeTo(3,3);
  31. fG.ResizeTo(5,5);
  32. fW.ResizeTo(3,3);
  33. fDaughtersInds.clear();
  34. Double_t pos[3] = {0.0}, magf[3] = {0.0};
  35. MpdKalmanFilter::Instance()->GetField(pos,magf);
  36. fieldConst = 0.3 * 0.01 * magf[2] / 10;
  37. fXY0[0] = fXY0[1] = 0.0;
  38. }
  39. //__________________________________________________________________________
  40. MpdParticle::MpdParticle(const MpdKalmanTrack& track, Int_t indx, Double_t mass, Double_t *orig)
  41. : TObject(track), fIndx(indx), fCharge(track.Charge()),
  42. fMass(mass), fChi2(0.0), fChi2ver(-1.0), fFlag(0), fPoint00(kTRUE)
  43. {
  44. /// Constructor from Kalman track
  45. fMeas.ResizeTo(5,1);
  46. fq.ResizeTo(3,1);
  47. fx.ResizeTo(3,1);
  48. fJ.ResizeTo(3,3);
  49. fJinv.ResizeTo(3,3);
  50. fD.ResizeTo(3,3);
  51. fE.ResizeTo(3,3);
  52. //fQ.ResizeTo(3,3);
  53. fA.ResizeTo(5,3);
  54. fB.ResizeTo(5,3);
  55. fC.ResizeTo(3,3);
  56. fG.ResizeTo(5,5);
  57. fW.ResizeTo(3,3);
  58. fDaughtersInds.clear();
  59. AddDaughter(indx);
  60. Double_t pos[3] = {0.0}, magf[3] = {0.0};
  61. MpdKalmanFilter::Instance()->GetField(pos,magf);
  62. fieldConst = 0.3 * 0.01 * magf[2] / 10;
  63. fXY0[0] = fXY0[1] = 0.0;
  64. //AZ if (orig) fPoint00 = kFALSE;
  65. Track2Part(track, kTRUE, orig);
  66. }
  67. //__________________________________________________________________________
  68. MpdParticle::MpdParticle (const MpdParticle& part)
  69. : TObject(part), fIndx(part.fIndx), fPdg(part.fPdg), fCharge(part.fCharge),
  70. fMass(part.fMass), fieldConst(part.fieldConst), fMeas(part.fMeas),
  71. fq(part.fq), fx(part.fx), fJ(part.fJ), fJinv(part.fJinv), fD(part.fD),
  72. fE(part.fE), fA(part.fA), fB(part.fB), fC(part.fC), fG(part.fG),
  73. fW(part.fW), fChi2(part.fChi2), fChi2ver(part.fChi2ver), fFlag(part.fFlag),
  74. fPoint00(part.fPoint00)
  75. {
  76. ///copy constructor
  77. fXY0[0] = part.fXY0[0];
  78. fXY0[1] = part.fXY0[1];
  79. fDaughtersInds = part.fDaughtersInds;
  80. }
  81. //__________________________________________________________________________
  82. MpdParticle & MpdParticle::operator=(const MpdParticle& part)
  83. {
  84. /// Asignment operator
  85. // check assignement to self
  86. if (this == &part) return *this;
  87. // base class assignement
  88. TObject::operator=(part);
  89. fIndx = part.fIndx; fPdg = part.fPdg; fCharge = part.fCharge;
  90. fMass = part.fMass; fieldConst = part.fieldConst; fMeas = part.fMeas;
  91. fq = part.fq; fx = part.fx; fJ = part.fJ; fJinv = part.fJinv;
  92. fD = part.fD; fE = part.fE; fA = part.fA; fB = part.fB;
  93. fC = part.fC; fG = part.fG; fW = part.fW; fChi2 = part.fChi2;
  94. fChi2ver = part.fChi2ver; fFlag = part.fFlag;
  95. fXY0[0] = part.fXY0[0];
  96. fXY0[1] = part.fXY0[1];
  97. fDaughtersInds = part.fDaughtersInds;
  98. return *this;
  99. }
  100. //__________________________________________________________________________
  101. MpdParticle::~MpdParticle()
  102. {
  103. /// Destructor
  104. }
  105. /*
  106. //__________________________________________________________________________
  107. Int_t MpdParticle::Compare(const TObject* part) const
  108. {
  109. /// "Compare" function to sort in descending order in pt
  110. MpdKalmanTrack *trackKF = (MpdKalmanTrack*) track;
  111. Double_t pt = 1. / TMath::Abs(trackKF->GetParam(4));
  112. Double_t ptthis = 1. / TMath::Abs((*fParam)(4,0));
  113. if (ptthis < pt) return 1;
  114. else if (ptthis > pt) return -1;
  115. return 0;
  116. }
  117. */
  118. //__________________________________________________________________________
  119. void MpdParticle::SetMass (Double_t mass)
  120. {
  121. /// Set particle mass (if negative, use PDG table value)
  122. fMass = (mass > -1.0) ? mass : TDatabasePDG::Instance()->GetParticle(fPdg)->Mass();
  123. }
  124. //__________________________________________________________________________
  125. void MpdParticle::Track2Part(const MpdKalmanTrack &track, Bool_t setWeight, Double_t *orig)
  126. {
  127. /// Different parameterization for MpdParticle as compared to MpdKalmanTrack
  128. /// (see R.Luchsinger, Ch.Grab CPC 76 (1993) 263-280), but also different
  129. /// from the described in the paper - different order of parameters in order to be
  130. /// more consistent with MpdKalmanTrack
  131. // !!! Check for modular tracking (cases with local coordinates (nodes)) !!!
  132. Double_t vert[3] = {0};
  133. if (orig == NULL) orig = vert;
  134. fMeas(4,0) = track.GetParam(4) * fieldConst; // k
  135. Double_t phi = track.GetParam(2);
  136. fMeas(2,0) = TMath::ATan2 (TMath::Sin(phi),TMath::Cos(phi)); // phi: -pi:pi
  137. fMeas(3,0) = track.Theta(); // theta
  138. fMeas(0,0) = track.GetPosNew(); // DCA w.r.t. (0,0)
  139. fMeas(1,0) = track.GetParam(1); // z at DCA
  140. if (fMeas(0,0) > 1.e-7) {
  141. // Compute sign of DCA (DCA x Pt)
  142. phi = track.GetParam(0) / fMeas(0,0);
  143. fXY0[0] = fMeas(0,0) * TMath::Cos(phi);
  144. fXY0[1] = fMeas(0,0) * TMath::Sin(phi);
  145. Double_t dx = fXY0[0] - orig[0];
  146. Double_t dy = fXY0[1] - orig[1];
  147. TVector3 dcaV(dx,dy,0.0);
  148. TVector3 momV(TMath::Cos(fMeas(2,0)),TMath::Sin(fMeas(2,0)),0.0);
  149. fMeas(0,0) = TMath::Sqrt(dx*dx+dy*dy) * TMath::Sign(1.0,dcaV.Cross(momV).Z()); // signed DCA w.r.t. "orig"
  150. //AZ
  151. fXY0[0] = dx;
  152. fXY0[1] = dy;
  153. }
  154. //cout << fXY0[0] << " " << -fMeas(0,0)*TMath::Cos(fMeas(2,0)+TMath::PiOver2()) << " "
  155. // << fXY0[1] << " " << -fMeas(0,0)*TMath::Sin(fMeas(2,0)+TMath::PiOver2()) << endl;
  156. if (!setWeight) return; // skip all the rest
  157. fq(0,0) = fMeas(2,0);
  158. fq(1,0) = fMeas(3,0);
  159. fq(2,0) = fMeas(4,0);
  160. // Obtain weight matrix at DCA
  161. MpdKalmanTrack track1 = track;
  162. if (fPoint00) {
  163. // For tracks extrapolated to the point (0,0)
  164. track1.SetParamNew(*track1.GetParam());
  165. track1.SetPos(track1.GetPosNew());
  166. track1.ReSetWeight();
  167. TMatrixD g = *track1.GetWeight(); // track weight matrix
  168. //cout << " track Weight " << endl;
  169. //g.Print();
  170. // Propagate to the point where track covar. matrix was stored
  171. MpdKalmanHit hit;
  172. if (track.GetNode() == "") {
  173. hit.SetType(MpdKalmanHit::kFixedR);
  174. hit.SetPos(track.GetPos());
  175. } else {
  176. //cout << " !!! Not implemented for local coordinates yet !!! " << endl;
  177. //exit(0);
  178. hit.SetType(MpdKalmanHit::kFixedP);
  179. TString detName = track.GetNode();
  180. if (track.GetUniqueID()) {
  181. // ITS
  182. detName = detName(16,detName.Length());
  183. detName += "#0";
  184. }
  185. MpdKalmanGeoScheme *geo = MpdKalmanFilter::Instance()->GetGeo();
  186. hit.SetDetectorID(geo->DetId(detName));
  187. // Find distance from the current track position to the last point (plane) -
  188. // to define direction (mainly for ITS)
  189. TVector3 pos = geo->GlobalPos(&hit);
  190. TVector3 norm = geo->Normal(&hit);
  191. Double_t v7[7] = {0.0};
  192. TString node = track1.GetNode();
  193. track1.SetNode("");
  194. MpdKalmanFilter::Instance()->SetGeantParamB(&track1, v7, 1);
  195. Double_t d = -(pos * norm); // Ax+By+Cz+D=0, A=nx, B=ny, C=nz
  196. TVector3 v3(v7[0], v7[1], v7[2]);
  197. d += v3 * norm;
  198. if (d < 0) track1.SetDirection(MpdKalmanTrack::kOutward);
  199. }
  200. MpdKalmanFilter::Instance()->PropagateToHit(&track1,&hit,kFALSE,kTRUE);
  201. track1.SetWeight(g);
  202. /*
  203. TMatrixD jacob = MpdKalmanFilter::Instance()->GetJacob(); // propagation Jacobian
  204. jacob.Invert();
  205. // Weight propagation to DCA
  206. TMatrixD tmp(g,TMatrixD::kMult,jacob); // WD
  207. TMatrixD weight1(jacob,TMatrixD::kTransposeMult,tmp); // DtWD
  208. g.Print();
  209. weight1.Print();
  210. // Apply Jacobian
  211. */
  212. } // if (fPoint00)
  213. track1.SetPos(track1.GetPosNew());
  214. //MpdMotherFitterPart::Instance()->WeightAtDca(this, track1, orig);
  215. WeightAtDca(track1, orig);
  216. //cout << " Weight at DCA " << endl;
  217. //fG.Print();
  218. }
  219. //__________________________________________________________________________
  220. void MpdParticle::WeightAtDca(MpdKalmanTrack &tr, Double_t *vert)
  221. {
  222. // Obtain MpdParticle weight at DCA from MpdKalmanTrack weight
  223. TMatrixDSym *covar = tr.Weight2Cov(); // get covariance matrix
  224. TMatrixD g = *tr.GetWeight(); // track weight matrix
  225. TMatrixD meas0 = GetMeas();
  226. TMatrixD param0 = *tr.GetParamNew();
  227. if (tr.GetNodeNew() != "") {
  228. // Local coordinates
  229. tr.SetNode(tr.GetNodeNew());
  230. tr.SetNodeNew("");
  231. }
  232. TMatrixD jacob(5,5);
  233. //Double_t vert[3] = {0.0};
  234. Double_t dPar;
  235. // Loop over parameters to find change of the propagated vs initial ones
  236. Bool_t ok = kTRUE;
  237. for (Int_t i = 0; i < 5; ++i) {
  238. dPar = TMath::Sqrt((*covar)(i,i));
  239. if (i < 4) dPar = TMath::Min (dPar, 0.1);
  240. else dPar = TMath::Min (dPar, 0.1*TMath::Abs(param0(4,0)));
  241. tr.SetParam(param0);
  242. if (i == 4) dPar *= TMath::Sign(1.,-param0(4,0)); // 1/p
  243. //else if (i == 2) dPar *= sign;
  244. else if (i == 3) dPar *= TMath::Sign(1.,-param0(3,0)); // dip-angle
  245. tr.SetParam(i,param0(i,0)+dPar);
  246. ok = MpdKalmanFilter::Instance()->FindPca(&tr, vert);
  247. tr.SetParam(*tr.GetParamNew());
  248. MpdParticle tmpPart;
  249. tmpPart.Track2Part(tr, kFALSE, vert);
  250. Double_t meas = 0.0;
  251. for (Int_t j = 0; j < 5; ++j) {
  252. if (j == 2) meas = MpdKalmanFilter::Instance()->Proxim(meas0(j,0),tmpPart.GetMeas(j));
  253. else meas = tmpPart.GetMeas(j);
  254. jacob(j,i) = (meas - meas0(j,0)) / dPar;
  255. //cout << i << " " << j << " " << dPar << " " << track->GetParamNew(j) << " " << paramNew0(j,0)
  256. // << " " << track->GetPos() << " " << track->GetPosNew() << " " << jacob(j,i) << endl;
  257. }
  258. }
  259. //jacob.Print();
  260. jacob.Invert();
  261. TMatrixD tmp(g,TMatrixD::kMult,jacob); // WD
  262. TMatrixD weight1(jacob,TMatrixD::kTransposeMult,tmp); // DtWD
  263. //part->SetG (weight1);
  264. /*AZ
  265. if (part->Point00()) {
  266. // Add multiple scattering
  267. weight1.Invert(); // covar
  268. Double_t step = tr.GetPos(); // radius of TPC inner layer
  269. Double_t x0 = 13363.6; // rad. length - TPCMixture
  270. TString mass2;
  271. mass2 += (part->GetMass() * part->GetMass());
  272. Double_t angle2 = MpdKalmanFilter::Instance()->Scattering(&tr, x0, step, mass2);
  273. Double_t th = tr.GetParamNew(3);
  274. Double_t cosTh = TMath::Cos(th);
  275. weight1(2,2) += (angle2 / cosTh / cosTh);
  276. weight1(3,3) += angle2;
  277. weight1.Invert(); // weight
  278. }
  279. */
  280. SetG (weight1);
  281. }
  282. //__________________________________________________________________________
  283. Double_t MpdParticle::BuildMother(vector<MpdParticle*> &vDaught)
  284. {
  285. /// Build mother particle from daughters which were smoothed
  286. /// according to the decay vertex constraint (after FindVertex)
  287. return MpdMotherFitterPart::Instance()->BuildMother (this, vDaught);
  288. }
  289. //__________________________________________________________________________
  290. Double_t MpdParticle::BuildMother(vector<MpdKalmanTrack*> &vTracks, vector<MpdParticle*> &vDaught)
  291. {
  292. /// Build mother particle from daughters which were smoothed
  293. /// according to the decay vertex constraint (after FindVertex).
  294. /// Daughters are built from tracks and parametrized at their
  295. /// intersection point.
  296. return MpdMotherFitterPart::Instance()->BuildMother(this, vTracks, vDaught);
  297. }
  298. //__________________________________________________________________________
  299. void MpdParticle::FillJ()
  300. {
  301. /// Fill Jacobian matrix fJ
  302. Double_t p0 = Momentum();
  303. Double_t pt = Pt();
  304. Double_t ph0 = Phi();
  305. Double_t th0 = Theta();
  306. Double_t k0 = fq(2,0);
  307. Double_t cosPh = TMath::Cos(ph0);
  308. Double_t sinPh = TMath::Sin(ph0);
  309. Double_t sinTh = TMath::Sin(th0);
  310. Double_t cosTh = TMath::Cos(th0);
  311. fJ = 0.0;
  312. if (fCharge) {
  313. // Charged particle
  314. // d/d(phi)
  315. fJ(0,0) = -pt * sinPh;
  316. fJ(1,0) = pt * cosPh;
  317. // d/d(theta)
  318. fJ(2,1) = -pt / sinTh / sinTh;
  319. // d/dk
  320. fJ(0,2) = -pt / k0 * cosPh;
  321. fJ(1,2) = -pt / k0 * sinPh;
  322. //fJ(2,2) = -pt / k0 / TMath::Tan(th0);
  323. Double_t coTan = cosTh / sinTh;
  324. fJ(2,2) = -pt / k0 * coTan;
  325. } else {
  326. // neutral - check !!!
  327. // d/d(phi)
  328. fJ(0,0) = -p0 * sinTh * sinPh;
  329. fJ(1,0) = p0 * sinTh * cosPh;
  330. // d/d(theta)
  331. fJ(0,1) = p0 * cosTh * cosPh;
  332. fJ(1,1) = p0 * cosTh * sinPh;
  333. fJ(2,1) = -p0 * sinTh;
  334. // d/dp
  335. fJ(0,2) = sinTh * cosPh;
  336. fJ(1,2) = sinTh * sinPh;
  337. fJ(2,2) = cosTh;
  338. }
  339. TMatrixD tmp(TMatrixD::kTransposed,fJ);
  340. //fJ = tmp;
  341. }
  342. //__________________________________________________________________________
  343. void MpdParticle::FillJinv(TVector3& mom3)
  344. {
  345. /// Fill Jacobian matrix fJinv
  346. Double_t px0 = mom3.X();
  347. Double_t py0 = mom3.Y();
  348. Double_t pz0 = mom3.Z();
  349. Double_t p0 = mom3.Mag();
  350. Double_t pt0 = mom3.Pt();
  351. //Double_t ph0 = mom3.Phi();
  352. //Double_t th0 = mom3.Theta();
  353. //Double_t k0 = Theta();
  354. fJinv = 0.0;
  355. if (fCharge) {
  356. // Charged mother - check !!!
  357. // d/d(px)
  358. fJinv(0,0) = -py0 / pt0 / pt0;
  359. fJinv(1,0) = px0 * pz0 / p0 / p0 / pt0;
  360. fJinv(2,0) = fCharge * fieldConst * px0 / pt0 / pt0 / pt0;
  361. // d/d(py)
  362. fJinv(0,1) = px0 / pt0 / pt0;
  363. fJinv(1,1) = py0 * pz0 / p0 / p0 / pt0;
  364. fJinv(2,1) = fCharge * fieldConst * py0 / pt0 / pt0 / pt0;
  365. // d/d(pz)
  366. fJinv(1,2) = -pt0 / p0 / p0;
  367. } else {
  368. // neutral mother
  369. // d/d(px)
  370. fJinv(0,0) = -py0 / pt0 / pt0;
  371. fJinv(1,0) = px0 * pz0 / p0 / p0 / pt0;
  372. fJinv(2,0) = px0 / p0;
  373. // d/d(py)
  374. fJinv(0,1) = px0 / pt0 / pt0;
  375. fJinv(1,1) = py0 * pz0 / p0 / p0 / pt0;
  376. fJinv(2,1) = py0 / p0;
  377. // d/d(pz)
  378. fJinv(1,2) = -pt0 / p0 / p0;
  379. fJinv(2,2) = pz0 / p0;
  380. }
  381. TMatrixD tmp(TMatrixD::kTransposed,fJinv);
  382. //fJinv = tmp;
  383. }
  384. //__________________________________________________________________________
  385. Double_t MpdParticle::Chi2Vertex(MpdVertex *vtx)
  386. {
  387. /// Compute Chi2 w.r.t. vertex
  388. fChi2ver = MpdMotherFitterPart::Instance()->Chi2Vertex(this, vtx);
  389. return fChi2ver;
  390. }
  391. //__________________________________________________________________________
  392. Double_t MpdParticle::Energy() const
  393. {
  394. /// Return particle energy
  395. Double_t mom = Momentum();
  396. return TMath::Sqrt (mom * mom + fMass * fMass);
  397. }
  398. //__________________________________________________________________________
  399. Double_t MpdParticle::Rapidity() const
  400. {
  401. /// Return particle rapidity
  402. TVector3 mom3 = Momentum3();
  403. Double_t mom = Momentum();
  404. Double_t e = TMath::Sqrt (mom * mom + fMass * fMass);
  405. Double_t pz = mom3.Z();
  406. Double_t y = 0.5 * TMath::Log ( (e+pz) / (e-pz) );
  407. return y;
  408. }
  409. //__________________________________________________________________________
  410. void MpdParticle::Print()
  411. {
  412. /// Print particle info
  413. }
  414. ClassImp(MpdParticle)