MpdMotherFitterSimple.cxx 9.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317
  1. /// \class MpdMotherFitterSimple
  2. ///
  3. /// Simple mother fitter for the MPD detector (using helices and straight lines)
  4. /// \author Alexander Zinchenko (LHEP, JINR, Dubna)
  5. #include "MpdMotherFitterSimple.h"
  6. #include "MpdKalmanFilter.h"
  7. #include "MpdKalmanGeoScheme.h"
  8. #include "MpdKalmanTrack.h"
  9. #include "MpdKalmanHit.h"
  10. #include "MpdHelix.h"
  11. #include "FairField.h"
  12. #include "FairMCTrack.h"
  13. #include "FairRunAna.h"
  14. #include "FairTask.h"
  15. #include <TMath.h>
  16. #include <TMatrixD.h>
  17. #include <TMinuit.h>
  18. #include <TVector3.h>
  19. #include <iostream>
  20. using std::cout;
  21. using std::endl;
  22. MpdMotherFitterSimple* MpdMotherFitterSimple::fgMF = 0x0;
  23. //__________________________________________________________________________
  24. MpdMotherFitterSimple::MpdMotherFitterSimple()
  25. : FairTask()
  26. {
  27. /// Default constructor
  28. for (Int_t i = 0; i < 3; ++i) fTrC[i] = NULL;
  29. }
  30. //__________________________________________________________________________
  31. MpdMotherFitterSimple::MpdMotherFitterSimple(const char *name, const char *title)
  32. : FairTask(name)
  33. {
  34. /// Constructor
  35. for (Int_t i = 0; i < 3; ++i) fTrC[i] = NULL;
  36. fgMF = this;
  37. }
  38. //__________________________________________________________________________
  39. MpdMotherFitterSimple* MpdMotherFitterSimple::Instance()
  40. {
  41. /// Get pointer to the mother fitter singleton object
  42. if (!fgMF){
  43. fgMF = new MpdMotherFitterSimple;
  44. fgMF->Init();
  45. // automatic destruction is forced
  46. std::atexit(DestroyInstance);
  47. }
  48. return fgMF;
  49. }
  50. //__________________________________________________________________________
  51. MpdMotherFitterSimple* MpdMotherFitterSimple::Instance(const char *name, const char *title)
  52. {
  53. /// Get pointer to the mother fitter singleton object
  54. if (!fgMF){
  55. fgMF = new MpdMotherFitterSimple(name, title);
  56. fgMF->Init();
  57. // automatic destruction is forced
  58. std::atexit(DestroyInstance);
  59. }
  60. return fgMF;
  61. }
  62. //__________________________________________________________________________
  63. MpdMotherFitterSimple::~MpdMotherFitterSimple()
  64. {
  65. /// Destructor
  66. //FairRootManager *manager = FairRootManager::Instance();
  67. //manager->Write();
  68. for (Int_t i = 0; i < 3; ++i) delete fTrC[i];
  69. }
  70. //__________________________________________________________________________
  71. InitStatus MpdMotherFitterSimple::Init() {
  72. cout << "InitStatus MpdMotherFitterSimple::Init\n\n";
  73. Double_t bZ = 5.;
  74. FairField * magField = FairRunAna::Instance()->GetField();
  75. if (!magField || TMath::Abs(magField->GetBz(0,0,0)) < 0.01) {
  76. cout << " -I- MpdMotherFitterSimple::Init - Using the default constant magnetic field Bz = 5 kG " << endl;
  77. } else bZ = TMath::Abs(magField->GetBz(0,0,0));
  78. fieldConst = 0.3 * 0.01 * bZ / 10;
  79. for (Int_t i = 0; i < 3; ++i) fTrC[i] = new MpdHelix(0.0, 0.0, 0.0, TVector3(0,0,0));
  80. if (!gMinuit) new TMinuit(2);
  81. return kSUCCESS;
  82. }
  83. //__________________________________________________________________________
  84. InitStatus MpdMotherFitterSimple::ReInit()
  85. {
  86. //fMagField = FairRunAna::Instance()->GetField(); // !!! interim solution !!!
  87. return kSUCCESS;
  88. }
  89. //__________________________________________________________________________
  90. void MpdMotherFitterSimple::Reset()
  91. {
  92. ///
  93. //cout << " MpdMotherFitterSimple::Reset " << endl;
  94. }
  95. //__________________________________________________________________________
  96. void MpdMotherFitterSimple::Register()
  97. {
  98. ///
  99. //FairRootManager::Instance()->
  100. //Register("TpcLheKalmanTrack","TpcLheKalmanTrack", fTracks, kFALSE);
  101. }
  102. //__________________________________________________________________________
  103. void MpdMotherFitterSimple::Finish()
  104. {
  105. ///
  106. }
  107. //__________________________________________________________________________
  108. void MpdMotherFitterSimple::Exec(Option_t * option)
  109. {
  110. ///
  111. }
  112. //__________________________________________________________________________
  113. TLorentzVector MpdMotherFitterSimple::BuildMother(vector<Double_t> &mass, vector<MpdKalmanTrack*> &vDaughtTr,
  114. vector<TLorentzVector*> &vDaughtLv)
  115. {
  116. // Build mother from tracks and/or Lorentz vectors
  117. Int_t nDaughtTr = vDaughtTr.size();
  118. Int_t nDaughtLv = vDaughtLv.size();
  119. TVector3 vtx;
  120. TLorentzVector mother;
  121. MpdKalmanTrack *tr = NULL;
  122. if (nDaughtLv == 0) {
  123. // Only tracks
  124. FindVertex(vDaughtTr);
  125. for (Int_t i = 0; i < nDaughtTr; ++i) {
  126. //tr = vDaughtTr[i];
  127. Double_t pt = Pt(fTrC[i]);
  128. Double_t phi = Phi(fTrC[i]);
  129. Double_t dip = fTrC[i]->dipAngle();
  130. TLorentzVector part;
  131. part.SetXYZM (pt*TMath::Cos(phi), pt*TMath::Sin(phi), pt*TMath::Tan(dip), mass[i]);
  132. mother += part;
  133. }
  134. return mother;
  135. }
  136. // At present this part works only for one neutral (LorentzVector) and one charged track
  137. // (to reconstruct cascades)
  138. // Create MpdHelix
  139. tr = vDaughtTr[0];
  140. Double_t r = tr->GetPosNew();
  141. Double_t phi = tr->GetParam(0) / r;
  142. Double_t x = r * TMath::Cos(phi);
  143. Double_t y = r * TMath::Sin(phi);
  144. Double_t dip = tr->GetParam(3);
  145. Double_t cur = TMath::Abs (tr->GetParam(4)) * fieldConst;
  146. TVector3 o(x, y, tr->GetParam(1));
  147. Int_t h = (Int_t) TMath::Sign(1.1,tr->GetParam(4));
  148. fTrC[0]->setParameters(cur, dip, tr->GetParam(2)-TMath::PiOver2()*h, o, h);
  149. /*
  150. Double_t s = trC.pathLength(primVert);
  151. TVector3 pca = trC.at(s);
  152. pca -= primVert;
  153. // !!! Cut on K-
  154. if (tr->GetChi2Vertex() < gChi2K || pca.Mag() < gDcaK) continue;
  155. */
  156. // Find decay vertex of a neutral and charged track
  157. TLorentzVector lambda = *vDaughtLv.front();
  158. fMomN = lambda.Vect();
  159. fMomN *= (1. / fMomN.Mag());
  160. //vtxN = vtxL0; ///
  161. //void SetFCN(void (*fcn)(Int_t &, Double_t *, Double_t &f, Double_t *, Int_t))
  162. gMinuit->SetFCN(FcnDist);
  163. Double_t arglist[10];
  164. Int_t ierflg = 0;
  165. arglist[0] = 1;
  166. gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
  167. arglist[0] = -1; //1; //-1;
  168. gMinuit->mnexcm("SET PRINT", arglist, 1, ierflg);
  169. Double_t vstart[2] = {-0.1,-0.1};
  170. static Double_t step[2] = {0.1, 0.1};
  171. gMinuit->mnparm(0, "lengN", vstart[0], step[0], 0,0,ierflg);
  172. gMinuit->mnparm(1, "lengC", vstart[1], step[1], 0,0,ierflg);
  173. // Now ready for minimization step
  174. arglist[0] = 500;
  175. arglist[1] = 1.;
  176. gMinuit->mnexcm("MIGRAD", arglist, 2, ierflg);
  177. // Print results
  178. //*
  179. Double_t amin,edm,errdef;
  180. Int_t nvpar,nparx,icstat;
  181. gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
  182. fDist = amin;
  183. //gMinuit->mnprin(3,amin);
  184. //*/
  185. // !!! Cut on distance between daughters
  186. //if (amin > gDistLK) continue;
  187. Double_t lengC, lengN, parErr;
  188. gMinuit->GetParameter(0,lengN,parErr);
  189. gMinuit->GetParameter(1,lengC,parErr);
  190. TVector3 posC = fTrC[0]->at(lengC);
  191. TVector3 posN = TVector3(fVtx) + fMomN * lengN;
  192. TVector3 posOm = (posN + posC) * 0.5;
  193. posOm.GetXYZ(fVtx);
  194. Double_t pt = TMath::Abs (1. / tr->GetParam(4));
  195. Double_t s = fTrC[0]->pathLength(posC); // PCA with decay vertex
  196. fTrC[0]->moveOrigin(s);
  197. phi = Phi(fTrC[0]);
  198. dip = tr->GetParam(3);
  199. TLorentzVector part;
  200. part.SetXYZM (pt*TMath::Cos(phi), pt*TMath::Sin(phi), pt*TMath::Tan(dip), mass[0]);
  201. mother = lambda;
  202. mother += part;
  203. return mother;
  204. }
  205. //____________________________________________________________________________
  206. void FcnDist(Int_t &npar, Double_t *gin, Double_t &rf, Double_t *par, Int_t iflag)
  207. {
  208. // Compute distance between straight line and helix
  209. TVector3 mom = MpdMotherFitterSimple::Instance()->GetMomN();
  210. mom *= par[0];
  211. TVector3 posN = TVector3(MpdMotherFitterSimple::Instance()->GetVertex());
  212. posN += mom;
  213. TVector3 posC = MpdMotherFitterSimple::Instance()->GetHelix(0)->at(par[1]);
  214. posC -= posN;
  215. rf = posC.Mag();
  216. //cout << par[0] << " " << par[1] << " " << f << endl;
  217. }
  218. //__________________________________________________________________________
  219. void MpdMotherFitterSimple::FindVertex(vector<MpdKalmanTrack*> vDaught)
  220. {
  221. /// Find decay vertex position as PCA of helices
  222. // Create MpdHelix's
  223. Int_t nDaught = vDaught.size();
  224. if (nDaught > 2) Fatal("MpdMotherFitterSimple::FindVertex"," Too many tracks !!!");
  225. TVector3 sum, p1, p2;
  226. for (Int_t i = 0; i < nDaught; ++i) {
  227. Double_t r = vDaught[i]->GetPosNew();
  228. Double_t phi = vDaught[i]->GetParam(0) / r;
  229. Double_t x = r * TMath::Cos(phi);
  230. Double_t y = r * TMath::Sin(phi);
  231. Double_t dip = vDaught[i]->GetParam(3);
  232. Double_t cur = TMath::Abs (vDaught[i]->GetParam(4)) * fieldConst;
  233. TVector3 o(x, y, vDaught[i]->GetParam(1));
  234. Int_t h = (Int_t) TMath::Sign(1.1,vDaught[i]->GetParam(4));
  235. //helix[itr++] = MpdHelix(cur, dip, phi, o, h);
  236. //setParameters(double c, double dip, double phase, const TVector3 o, int h);
  237. fTrC[i]->setParameters(cur, dip, vDaught[i]->GetParam(2)-TMath::PiOver2()*h, o, h);
  238. //cout << helix[i]->xcenter() << " " << helix[i]->ycenter() << endl;
  239. //sum += o;
  240. }
  241. Int_t ncombs = 0, nD1 = nDaught - 1;
  242. Double_t pathMax = 0.0;
  243. for (Int_t i = 0; i < nD1; ++i) {
  244. for (Int_t j = 1; j < nDaught; ++j) {
  245. pair<Double_t,Double_t> paths = fTrC[i]->pathLengths(*fTrC[j]);
  246. //cout << " Intersection: " << paths.first << " " << paths.second << endl;
  247. p1 = fTrC[i]->at(paths.first);
  248. p2 = fTrC[j]->at(paths.second);
  249. sum += (p1 + p2);
  250. ncombs += 2;
  251. pathMax = TMath::Max (pathMax, TMath::Abs(paths.first));
  252. pathMax = TMath::Max (pathMax, TMath::Abs(paths.second));
  253. }
  254. }
  255. if (pathMax > 80.) sum.SetXYZ(0,0,0); // check constrains
  256. else sum *= (1./ncombs);
  257. sum.GetXYZ(fVtx);
  258. fDist = (p1 - p2).Mag();
  259. // Move helices to the vertex
  260. for (Int_t i = 0; i < nDaught; ++i) {
  261. Double_t s = fTrC[i]->pathLength(sum); // PCA with decay vertex
  262. fTrC[i]->moveOrigin(s);
  263. }
  264. }
  265. //__________________________________________________________________________
  266. Double_t MpdMotherFitterSimple::Phi(MpdHelix *helix) const
  267. {
  268. // Compute Phi of the track
  269. return helix->phase() + TMath::PiOver2() * helix->h(); // check !!!
  270. }
  271. //__________________________________________________________________________
  272. ClassImp(MpdMotherFitterSimple)