TpcLheTrackFitter.cxx 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517
  1. #include "TpcLheTrackFitter.h"
  2. #include "TpcLheTrack.h"
  3. #include "lhe.h"
  4. //#include "FairTrackParH.h"
  5. #include "FairField.h"
  6. #include "FairMCApplication.h"
  7. #include "FairRunAna.h"
  8. #include "FairRootManager.h"
  9. #include "TGeoTrack.h"
  10. #include "TH2F.h"
  11. using namespace std;
  12. // ---------------------------------------------------------------
  13. // --- Interface with TrackFinder and output ---
  14. TpcLheTrackFitter* TpcLheTrackFitter::ftInstance = NULL;
  15. //___________________________________________________________
  16. TpcLheTrackFitter* TpcLheTrackFitter::Instance() {
  17. //---
  18. return ftInstance;
  19. }
  20. //___________________________________________________________
  21. TpcLheTrackFitter::~TpcLheTrackFitter() {
  22. //
  23. FairRootManager *fManger =FairRootManager::Instance();
  24. fManger->Write();
  25. }
  26. //___________________________________________________________
  27. TpcLheTrackFitter::TpcLheTrackFitter() {
  28. //---
  29. if( !ftInstance ) ftInstance = this;
  30. }
  31. //___________________________________________________________
  32. TpcLheTrackFitter::TpcLheTrackFitter(const char *name, const char *title)
  33. :FairTask(name) {
  34. //---
  35. if( !ftInstance ) ftInstance = this;
  36. // fOutputTracks = new TClonesArray("TpcFFTrack");
  37. }
  38. //___________________________________________________________
  39. InitStatus TpcLheTrackFitter::Init() {
  40. // cout << "InitStatus TpcLheTrackFitter::Init()" << endl;
  41. FairRootManager *fManager =FairRootManager::Instance();
  42. fTpcHits = (TClonesArray *)fManager->GetObject("LheHit");
  43. fTpcPoints = (TClonesArray *)fManager->GetObject("TpcPoint");
  44. if ( ! fTpcPoints ) {
  45. cout << "-I- "<< GetName() << "::Init: No TpcPoint array!" << endl;
  46. return kERROR;
  47. }
  48. fTpcTracks = (TClonesArray *)fManager->GetObject("TpcLheTrack");
  49. if ( ! fTpcTracks ) {
  50. cout << "-I- "<< GetName() << "::Init: No TpcLheTrack array!" << endl;
  51. return kERROR;
  52. }
  53. fTpcTrCand = new TClonesArray("TGeoTrack");
  54. fTpcTrFit = new TClonesArray("TGeoTrack");
  55. fManager->Register("TrackCand", "Cand",fTpcTrCand, kTRUE);
  56. fManager->Register("TrackFit", "Fit", fTpcTrFit, kTRUE);
  57. // get the field in Memory:
  58. FairRunAna *fRun=FairRunAna::Instance();
  59. fMagField = (FairField*) fRun->GetField();
  60. fTrackCuts =TpcLheTrackCuts::Instance();
  61. fXYG= new TH2F("XY Geant","XY",100,0,100,100,0,100);
  62. fXYF= new TH2F("XY Fit","XY",100,0,100,100,0,100);
  63. fYZG= new TH2F("YZ Geant","YZ",100,0,100,100,0,100);
  64. fYZF= new TH2F("YZ Fit","YZ",100,0,100,100,0,100);
  65. fXZG= new TH2F("XZ Geant","XZ",100,0,100,100,0,100);
  66. fXZF= new TH2F("XZ Fit","XZ",100,0,100,100,0,100);
  67. return kSUCCESS;
  68. }
  69. //______________________________________________________
  70. void TpcLheTrackFitter::Exec(Option_t * option) {
  71. cout << " ===== TpcLheTrackFitter ===== " << endl;
  72. fTpcTrCand->Delete();
  73. fTpcTrFit->Delete();
  74. if (fTpcTracks) {
  75. cout << " Number of tracks for fitting " <<
  76. fTpcTracks->GetEntriesFast() << endl;
  77. }
  78. Int_t nTracks = fTpcTracks->GetEntriesFast();
  79. for (Int_t i = 0; i < nTracks; i++) {
  80. TpcLheTrack* track = (TpcLheTrack*) fTpcTracks->At(i);
  81. cout << "\n\n Track " << i << "\n";
  82. // cout << "\n\n Track " << track->GetTrackNumber() << "\n";
  83. // track->Print();
  84. TClonesArray& clref = *fTpcTrCand;
  85. Int_t size = clref.GetEntriesFast();
  86. fTrCan= new(clref[size]) TGeoTrack();
  87. TClonesArray& clref1 = *fTpcTrFit;
  88. Int_t size1 = clref1.GetEntriesFast();
  89. fTrFit= new(clref1[size1]) TGeoTrack();
  90. HelixFit(track);
  91. Info4Fit(track);
  92. }
  93. #if 0
  94. /** Constructor with all track variables (x,y,z in SC) **/
  95. FairTrackParH(Double_t x, Double_t y, Double_t z,
  96. Double_t lambda, Double_t phi, Double_t qp,
  97. Double_t CovMatrix[15]);
  98. /** Constructor track parameters with position (LAB) momentum **/
  99. FairTrackParH(TVector3 pos, TVector3 Mom, TVector3 posErr,
  100. TVector3 MomErr, Double_t q);
  101. #endif
  102. }
  103. //_____________________________________________________________________________
  104. void TpcLheTrackFitter::Info4Fit(TpcLheTrack *track) {
  105. //---
  106. TObjArray *rhits = (TObjArray* )track->GetRHits();
  107. Int_t nHits = rhits->GetEntries();
  108. // alpha = .2998 * magfield / 100 <- for obtain moment in GeV/c
  109. Double_t alpha = .2998 * .02 ;
  110. Double_t Q = double(TMath::Sign(1, track->GetCharge()));
  111. TpcLhePoint hel = track->GetCircle();
  112. Double_t xc = hel.GetX();
  113. Double_t yc = hel.GetY();
  114. Double_t rad = hel.GetZ();
  115. TVector2 centre(xc, yc);
  116. Double_t fi0 = centre.Phi(); // 0 - 2Pi
  117. Double_t d0 = Q*(TMath::Sqrt(xc*xc + yc*yc) - hel.GetZ());
  118. TVector2 pivot(-xc, -yc);
  119. Double_t lam = track->GetTanDipAngle();
  120. // now suppose that all tracks are primary
  121. Double_t z0 = 0; //track->GetZ0();
  122. // loop over hits
  123. for (int ih = 0; ih < nHits; ih++) {
  124. TpcLheHit* hit = (TpcLheHit*)rhits->At(ih);
  125. TpcPoint *point =(TpcPoint *) fTpcPoints->At(hit->GetRefIndex());
  126. cout << "\n Geant x, y, z: " << point->GetX() << " " << point->GetY()<< " " <<
  127. point->GetZ()
  128. << " px, py, pz: " << point->GetPx() << " " << point->GetPy() <<
  129. " " << point->GetPz() << endl;
  130. fTrCan->AddPoint(point->GetX(), point->GetY(), point->GetZ(),0 );
  131. // calculate deflection angle between point and phi0
  132. //
  133. fXYG->Fill(point->GetX(),point->GetY());
  134. fYZG->Fill(point->GetY(),point->GetZ());
  135. fXZG->Fill(point->GetX(),point->GetZ());
  136. TVector2 pnt(hit->GetX() - xc, hit->GetY() - yc);
  137. Double_t phi = pnt.DeltaPhi(pivot);
  138. // Double_t phi = pnt.DeltaPhi(centre);
  139. Double_t x = d0*TMath::Cos(fi0) +
  140. rad *(TMath::Cos(fi0) - TMath::Cos(fi0 + phi));
  141. Double_t y = d0*TMath::Sin(fi0) +
  142. rad *(TMath::Sin(fi0) - TMath::Sin(fi0 + phi));
  143. Double_t z = z0 - rad*lam*(phi)*TMath::Sign(-1.,phi);
  144. Double_t px = -Q*alpha*rad*TMath::Sin(fi0 + phi);
  145. Double_t py = Q*alpha*rad*TMath::Cos(fi0 + phi);
  146. Double_t pz = Q*alpha*rad*lam*TMath::Sign(-1.,phi);
  147. // write momentum into output tree, /sl/ 29.08.07
  148. track->SetPx(px);
  149. track->SetPy(py);
  150. track->SetPz(pz);
  151. cout << " Fitted x, y, z: " << x << " " << y << " " << z
  152. << " px, py, pz: " << px << " " << py << " " << pz << endl;
  153. fTrFit->AddPoint(x, y, z, 0 );
  154. fXYF->Fill(x,y);
  155. fYZF->Fill(y,z);
  156. fXZF->Fill(x,z);
  157. }
  158. }
  159. //_____________________________________________________________________________
  160. Int_t TpcLheTrackFitter::CircleFit(TpcLheTrack *track) {
  161. //--- From Ososkov CircleCOP() Comp.Phys.Com 33, p.329
  162. TObjArray *rhits = (TObjArray* )track->GetRHits();
  163. Int_t NHits = rhits->GetEntries();
  164. double M0,Mx,My;
  165. M0 = NHits;
  166. Mx=My=0.;
  167. for (Int_t lh = 0; lh < NHits; lh++) {
  168. TpcLheHit* hit = (TpcLheHit*) rhits->At(lh);
  169. Mx += hit->GetX();
  170. My += hit->GetY();
  171. }
  172. Mx /= M0;
  173. My /= M0;
  174. // computing moments (note: all moments are normed, i.e. divided
  175. // by N)
  176. double Xi,Yi,Zi;
  177. double Mxy, Mxx, Myy, Mxz, Myz, Mzz;
  178. Mxx = Myy = Mxy = Mxz = Myz = Mzz = 0.;
  179. for (Int_t lh = 0; lh < NHits; lh++) {
  180. TpcLheHit* hit = (TpcLheHit*) rhits->At(lh);
  181. // for (i=0; i<data.n; i++) {
  182. Xi = hit->GetX() - Mx;
  183. Yi = hit->GetY() - My;
  184. Zi = Xi*Xi + Yi*Yi;
  185. Mxy += Xi*Yi;
  186. Mxx += Xi*Xi;
  187. Myy += Yi*Yi;
  188. Mxz += Xi*Zi;
  189. Myz += Yi*Zi;
  190. Mzz += Zi*Zi;
  191. }
  192. Mxx /= M0;
  193. Myy /= M0;
  194. Mxy /= M0;
  195. Mxz /= M0;
  196. Myz /= M0;
  197. Mzz /= M0;
  198. // computing the coefficients of the characteristic polynomial
  199. double Mz,Mxz2,Myz2,Cov_xy; //,temp;
  200. double A0, A1, A2, A22, epsilon = 0.000000000001;
  201. double Dy, xnew, xold, ynew, yold = 100000000000.;
  202. Mz = Mxx + Myy;
  203. Cov_xy = Mxx*Myy - Mxy*Mxy;
  204. Mxz2 = Mxz*Mxz;
  205. Myz2 = Myz*Myz;
  206. A2 = 4.*Cov_xy - 3.*Mz*Mz - Mzz;
  207. A1 = Mzz*Mz + 4.*Cov_xy*Mz - Mxz2 - Myz2 - Mz*Mz*Mz;
  208. A0 = Mxz2*Myy + Myz2*Mxx - Mzz*Cov_xy - 2.*Mxz*Myz*Mxy + Mz*Mz*Cov_xy;
  209. A22 = A2 + A2;
  210. // iter = 0;
  211. xnew = 0.;
  212. // Newton's method starting at x=0
  213. int iter, iterMax = 20;
  214. for (iter=0; iter < iterMax; iter++) {
  215. ynew = A0 + xnew*(A1 + xnew*(A2 + 4.*xnew*xnew));
  216. if (fabs(ynew)>fabs(yold)) {
  217. // printf("Newton2 goes wrong direction: ynew=%f
  218. // yold=%f\n",ynew,yold);
  219. xnew = 0.;
  220. break;
  221. }
  222. Dy = A1 + xnew*(A22 + 16.*xnew*xnew);
  223. xold = xnew;
  224. xnew = xold - ynew/Dy;
  225. if (fabs((xnew-xold)/xnew) < epsilon) break;
  226. }
  227. if (iter == iterMax-1) {
  228. // printf("Newton2 does not converge in %d
  229. // iterations\n",iterMax);
  230. xnew = 0.;
  231. }
  232. if (xnew < 0.) {
  233. iter=30;
  234. // printf("Negative root: x=%f\n",xnew);
  235. }
  236. // computing the circle parameters
  237. double GAM,DET;
  238. double Xcenter,Ycenter,Radius;
  239. GAM = - Mz - xnew - xnew;
  240. DET = xnew*xnew - xnew*Mz + Cov_xy;
  241. Xcenter = (Mxz*(Myy-xnew) - Myz*Mxy)/DET/2.;
  242. Ycenter = (Myz*(Mxx-xnew) - Mxz*Mxy)/DET/2.;
  243. Radius = sqrt(Xcenter*Xcenter+Ycenter*Ycenter-GAM);
  244. // cout << " x center " << Xcenter + Mx <<
  245. // " y center " << Ycenter + My <<
  246. // " radius " << Radius << endl;
  247. track->SetCircle(Xcenter + Mx, Ycenter + My, Radius);
  248. return 1;
  249. }
  250. //_____________________________________________________________________________
  251. Int_t TpcLheTrackFitter::DeepFit(TpcLheTrack *track) {
  252. //--- line fit of array of points
  253. TObjArray *rHits = (TObjArray* )track->GetRHits();
  254. TpcLheHit *first = (TpcLheHit *)rHits->First();
  255. Double_t dx, dy ;
  256. dx = first->GetX() - track->GetVertex().GetX();
  257. dy = first->GetY() - track->GetVertex().GetY();;
  258. TpcLhePoint circ = track->GetCircle();
  259. Double_t radius = circ.GetZ();
  260. Double_t localPsi = 0.5F * sqrt ( dx*dx + dy*dy ) / radius ;
  261. Double_t total_s ;
  262. if ( fabs(localPsi) < 1. ) {
  263. total_s = 2.0 * radius * asin ( localPsi ) ;
  264. }
  265. else {
  266. total_s = 2.0 * radius * TMath::Pi() ;
  267. }
  268. Double_t sum = 0.F ;
  269. Double_t ss = 0.F ;
  270. Double_t sz = 0.F ;
  271. Double_t sss = 0.F ;
  272. Double_t ssz = 0.F ;
  273. Int_t nHits = rHits->GetEntries();
  274. Double_t dpsi, s;
  275. Double_t fS[nHits];
  276. Double_t fZWeight[nHits];
  277. for(Int_t i = 0; i < nHits; i++) {
  278. TpcLheHit *hit = (TpcLheHit *)rHits->At(i);
  279. fZWeight[i] = 1./(hit->GetZerr()*hit->GetZerr());
  280. if(i > 0) {
  281. TpcLheHit *last = (TpcLheHit *)rHits->At(i-1);
  282. dx = hit->GetX() - last->GetX();
  283. dy = hit->GetY() - last->GetY();
  284. dpsi = 0.5 * (Double_t)sqrt ( dx*dx + dy*dy ) / radius ;
  285. if(fabs(dpsi) > 1) return 1;
  286. // track->SetPsierr(dpsi);
  287. s = fS[i-1] - 2.0 * radius * (Double_t)asin ( dpsi ) ;
  288. fS[i]=s;
  289. }
  290. else
  291. fS[i] = total_s;
  292. sum += fZWeight[i];
  293. ss += fZWeight[i] * fS[i];
  294. sz += fZWeight[i] * hit->GetZ();
  295. sss += fZWeight[i] * fS[i] * fS[i];
  296. ssz += fZWeight[i] * fS[i] * hit->GetZ();
  297. }
  298. Double_t chi2;
  299. Double_t det = sum * sss - ss * ss;
  300. if ( fabs(det) < 1e-20) {
  301. chi2 = 99999.F ;
  302. //track->SetChiSq2(chi2);
  303. return 0 ;
  304. }
  305. // tan \lambda = dSdZ
  306. Double_t tanl = (Double_t)((sum * ssz - ss * sz ) / det );
  307. Double_t z0 = (Double_t)((sz * sss - ssz * ss ) / det );
  308. track->SetTanDipAngle(tanl);
  309. track->SetZ0(z0);
  310. // cout << "tanl " << tanl << " z0 " << z0 << endl;
  311. // calculate chi-square
  312. chi2 = 0.;
  313. Double_t r1 ;
  314. for(Int_t i=0; i < nHits; i++) {
  315. TpcLheHit *hit = (TpcLheHit *)rHits->At(i);
  316. r1 = hit->GetZ() - tanl * fS[i] - z0 ;
  317. chi2 += (fZWeight[i]) * (r1 * r1);
  318. }
  319. track->SetDipChi2(chi2);
  320. // should be done later
  321. // calculate estimated variance
  322. // varsq=chi/(double(n)-2.)
  323. // calculate covariance matrix
  324. // siga=sqrt(varsq*sxx/det)
  325. // sigb=sqrt(varsq*sum/det)
  326. //
  327. Double_t dtanl = (Double_t) ( sum / det );
  328. Double_t dz0 = (Double_t) ( sss / det );
  329. track->SetTanDipAngleErr(dtanl);
  330. track->SetZ0Err(dz0);
  331. cout << "soeren track " << track->GetPx() << "\n";
  332. return 0 ;
  333. }
  334. //_____________________________________________________________________________
  335. Int_t TpcLheTrackFitter::HelixFit(TpcLheTrack *track) {
  336. //--- Create helix as fit of array of points
  337. CircleFit(track);
  338. DeepFit(track);
  339. return 1;
  340. }
  341. //_____________________________________________________
  342. void TpcLheTrackFitter::Finish() {
  343. //---
  344. // FairRootManager *fManger =FairRootManager::Instance();
  345. // fManger->Fill();
  346. /* fXYF->Write();
  347. fYZF->Write();
  348. fXZF->Write();
  349. fXYG->Write();
  350. fYZG->Write();
  351. fXZG->Write();
  352. */
  353. printf("\n\n *** Finish ***");
  354. }
  355. ClassImp(TpcLheTrackFitter)