TpcLheTrackFinder.cxx 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585
  1. #include "TpcLheTrackFinder.h"
  2. #include "TpcLheCMTrack.h"
  3. #include "TpcLheCMPoint.h"
  4. #include "lhe.h"
  5. #include "FairMCApplication.h"
  6. #include "FairTask.h"
  7. #include "FairRunAna.h"
  8. #include "TpcGeoPar.h"
  9. #include "FairGeoNode.h"
  10. #include "FairGeoVector.h"
  11. #include "FairGeoMedium.h"
  12. #include "FairRootManager.h"
  13. #include "TObjectTable.h"
  14. #include "TClonesArray.h"
  15. ClassImp(TpcLheTrackFinder)
  16. //________________________________________________________________
  17. TpcLheTrackFinder::TpcLheTrackFinder() {
  18. //---
  19. fFoundTracks = new TClonesArray("TpcLheTrack");
  20. fCMHits = new TClonesArray("TpcLheCMPoint");
  21. fBench = new TBenchmark();
  22. fVertex = NULL;
  23. }
  24. //_________________________________________________________________
  25. TpcLheTrackFinder::
  26. TpcLheTrackFinder( const char *name, const char *title):FairTask(name) {
  27. //---
  28. fFoundTracks = new TClonesArray("TpcLheTrack");
  29. fCMHits = new TClonesArray("TpcLheCMPoint");
  30. fBench = new TBenchmark();
  31. fVertex = NULL;
  32. }
  33. //_________________________________________________________________
  34. TpcLheTrackFinder::~TpcLheTrackFinder() {
  35. if (fCMTracks) { fCMTracks->Delete(); delete fCMTracks; }
  36. if (fCMHits) { fCMHits->Delete(); delete fCMHits; }
  37. if (fSegments) delete fSegments;
  38. FairRootManager *fManager =FairRootManager::Instance();
  39. fManager->Write();
  40. }
  41. //_________________________________________________________________
  42. void TpcLheTrackFinder::Register() {
  43. //---
  44. FairRootManager::
  45. Instance()->Register("TpcLheTrack",
  46. "Lhe", fFoundTracks, kTRUE);
  47. FairRootManager::
  48. Instance()->Register("TpcLheCMPoint",
  49. "Lhe",fCMHits, kTRUE);
  50. }
  51. //________________________________________________________________
  52. InitStatus TpcLheTrackFinder::Init() {
  53. // ---
  54. FairRootManager *fManager = FairRootManager::Instance();
  55. fTrackCuts = TpcLheTrackCuts::Instance();
  56. if(fOption.Contains("geant")) {
  57. fGeantTracks = (TClonesArray *)fManager->GetObject("LheGeantTrack");
  58. }
  59. else {
  60. fLheHits = (TClonesArray *)fManager->GetObject("LheHit");
  61. }
  62. Register();
  63. // create TObjArrays
  64. fCMTracks = new TObjArray(64);
  65. fVertex = new TpcLhePoint(0.0, 0.0, 0.0);
  66. fSegment = new TpcLheSegments();
  67. fSegment->Init();
  68. fSegments = fSegment->GetSegments();
  69. return kSUCCESS;
  70. }
  71. //_________________________________________________________________
  72. void TpcLheTrackFinder::Exec(Option_t * option) {
  73. if(fOption.Contains("geant")) {
  74. cout << " Copy Geant tracks \n" << endl;
  75. CopyClones(fGeantTracks, fFoundTracks);
  76. return;
  77. }
  78. if (fBench) {
  79. fBench->Start("finder");
  80. }
  81. Reset();
  82. Int_t n_hits = fLheHits->GetEntriesFast(); // number of hits
  83. TpcLheHit *ghit = NULL;
  84. Int_t good_hits = 0;
  85. for (Int_t ih = 0; ih < n_hits; ih++) {
  86. ghit = (TpcLheHit *) fLheHits->UncheckedAt(ih);
  87. TClonesArray &cmhits = *fCMHits;
  88. TpcLheCMPoint *cmhit = new(cmhits[good_hits++]) TpcLheCMPoint(ghit);
  89. cmhit->SetHitNumber(ghit->GetHitNumber());
  90. cmhit->Setup(fVertex);
  91. cmhit->SetUsage(kFALSE);
  92. }
  93. cout << " Working with " << good_hits << " hits\n";
  94. fSegment->FillSegments(fCMHits);
  95. DoTracking();
  96. if(fBench) {
  97. cout << endl;
  98. fBench->Show("finder");
  99. }
  100. }
  101. //________________________________________________________________
  102. void TpcLheTrackFinder::AddTrackForFit(TpcLheCMTrack *track_in) {
  103. //--- Add a new track to the list of tracks for this event.
  104. TClonesArray &tracks = *fFoundTracks;
  105. Int_t size = tracks.GetEntriesFast();
  106. TObjArray *rhits = (TObjArray* )track_in->GetCMHits();
  107. cout << "TpcLheTrackFinder::AddTrackForFit "<< rhits << endl;
  108. if (rhits){
  109. Int_t NHits = rhits->GetEntriesFast();
  110. TpcLheTrack* track = new(tracks[size]) TpcLheTrack(size);
  111. track->SetCharge(track_in->GetCharge());
  112. for (Int_t lh = 0; lh < NHits; lh++) {
  113. TpcLheHit* hit = dynamic_cast <TpcLheHit*> (rhits->At(lh));
  114. track->AddHit(hit);
  115. }
  116. }
  117. }
  118. //_________________________________________________________________
  119. void TpcLheTrackFinder::Finish() {
  120. cout << " found "<< fFoundTracks->GetEntries() << " tracks\n";
  121. FairRootManager *fManager =FairRootManager::Instance();
  122. fManager->Fill();
  123. }
  124. //________________________________________________________________
  125. void TpcLheTrackFinder::Reset() {
  126. //---
  127. if (fCMTracks->GetEntriesFast() != 0) fCMTracks->Clear("C");
  128. if (fCMHits->GetEntriesFast() != 0) fCMHits->Clear("C");
  129. if (fFoundTracks->GetEntriesFast() != 0) fFoundTracks->Clear("C");
  130. }
  131. //________________________________________________________________
  132. void TpcLheTrackFinder::DoTracking() {
  133. //--- Tracking
  134. for (Int_t im = HIGH; im <= HIGH; im++ ) {
  135. switch (im) {
  136. case HIGH:
  137. fTrackCuts->SetHighMomTrackCuts();
  138. break;
  139. case LOW:
  140. fTrackCuts->SetLowMomTrackCuts();
  141. break;
  142. default:
  143. fTrackCuts->SetHighMomTrackCuts();
  144. }
  145. LoopOverHits();
  146. char info[2] = {' ','\0'};
  147. TrackingInfo(&info[0]);
  148. }
  149. }
  150. //________________________________________________________________
  151. void TpcLheTrackFinder::CheckClones() {
  152. //---
  153. TIter next(fCMTracks); next.Reset();
  154. TpcLheCMTrack* trk;
  155. Int_t ntracks = fCMTracks->GetEntriesFast();
  156. next.Reset();
  157. while((trk = (TpcLheCMTrack*) next())) {
  158. trk->SetGood(kTRUE);
  159. }
  160. for (Int_t it = 0;it < ntracks; it++) {
  161. trk = (TpcLheCMTrack*)fCMTracks->At(it);
  162. TpcLhePoint hel_ex = trk->GetCircle();
  163. Double_t xc_ex = hel_ex.GetX();
  164. Double_t yc_ex = hel_ex.GetY();
  165. Double_t rad_ex = hel_ex.GetZ();
  166. for (Int_t jt = it + 1; jt < ntracks; jt++) {
  167. TpcLheCMTrack* trk_in = (TpcLheCMTrack*) fCMTracks->At(jt);
  168. TpcLhePoint hel_in = trk_in->GetCircle();
  169. Double_t xc_in = hel_in.GetX();
  170. Double_t yc_in = hel_in.GetY();
  171. Double_t rad_in = hel_in.GetZ();
  172. if (TMath::Abs(xc_ex - xc_in) < 1. &&
  173. TMath::Abs(yc_ex - yc_in) < 1. &&
  174. TMath::Abs(rad_ex - rad_in) < 1.)
  175. trk_in->SetGood(kFALSE);
  176. // cout << "\n Clones\n"; trk->Print(); trk_in->Print();
  177. }
  178. }
  179. next.Reset();
  180. while((trk = (TpcLheCMTrack*) next())) {
  181. if (trk->IsGood()) AddTrackForFit(trk);
  182. }
  183. }
  184. //________________________________________________________________
  185. Bool_t TpcLheTrackFinder::TrackExtension(TpcLheCMTrack *track) {
  186. TpcLheCMPoint *hit = (TpcLheCMPoint *)(track->GetRHits())->First();
  187. Bool_t back;
  188. if (hit->GetZ() < 0.) {
  189. back = kTRUE;
  190. }
  191. else {
  192. back = kFALSE;
  193. }
  194. GetClosestHit(track, 0, back);
  195. if(!fTrackCuts->IsGoodFoundTrack(track)) {
  196. // cleanup
  197. track->Clear();
  198. return kFALSE;
  199. }
  200. else {
  201. return kTRUE;
  202. }
  203. }
  204. //________________________________________________________________
  205. void TpcLheTrackFinder::TrackingInfo(char *info) {
  206. // Information about the tracking process.
  207. cout << " found ";
  208. cout.width(5); cout << fFoundTracks->GetEntriesFast() << " tracks ";
  209. }
  210. //________________________________________________________________
  211. void TpcLheTrackFinder::LoopOverHits() {
  212. //--- loops over all hits
  213. Int_t entries = fCMHits->GetEntriesFast();
  214. for (Int_t nhit = 0; nhit < entries; nhit++) {
  215. TpcLheCMPoint *hit = (TpcLheCMPoint *)fCMHits->At(nhit);
  216. hit->Print();
  217. // start hit was used before
  218. if (!hit->GetUsage() ) {
  219. CreateTrack(hit);
  220. }
  221. }
  222. CheckClones();
  223. }
  224. //_______________________________________________________________________
  225. void TpcLheTrackFinder::CreateTrack(TpcLheCMPoint *seed_hit) {
  226. //---
  227. Double_t theta_cut = 3.*.0002;
  228. Int_t theta_low = fSegment->GetThetaSegm(seed_hit->GetTheta() - theta_cut);
  229. Int_t theta_high = fSegment->GetThetaSegm(seed_hit->GetTheta() + theta_cut);
  230. Double_t phi_cut = 3.*.004;
  231. Int_t phi_low = fSegment->GetPhiSegm(seed_hit->GetPhi() - phi_cut);
  232. Int_t phi_high = fSegment->GetPhiSegm(seed_hit->GetPhi() + phi_cut);
  233. for (Int_t theta_segm = theta_low; theta_segm <= theta_high;
  234. theta_segm++) {
  235. // phi
  236. for (Int_t phi_segm = phi_low; phi_segm <= phi_high;
  237. phi_segm++) {
  238. // loop over entries in one sub segment
  239. TObjArray *cell = (TObjArray *)fSegments->At
  240. // (fSegment->GetSegm(theta_segm, phi_segm));
  241. (fSegment->GetSegm(theta_segm, phi_segm));
  242. Int_t cell_entries = cell->GetEntriesFast();
  243. if (cell_entries) {
  244. for (Int_t hit_num = 0; hit_num < cell_entries; hit_num++) {
  245. TpcLheCMPoint *hit = (TpcLheCMPoint *)cell->At(hit_num);
  246. if( TMath::Abs(seed_hit->GetTheta() - hit->GetTheta()) >
  247. theta_cut) continue;
  248. if( TMath::Abs(seed_hit->GetPhi() - hit->GetPhi()) >
  249. phi_cut) continue;
  250. if (!hit->GetUsage() &&
  251. hit->GetHitNumber() != seed_hit->GetHitNumber()) {
  252. Int_t ntracks = fCMTracks->GetEntries();
  253. fCMTracks->AddAtAndExpand(new TpcLheCMTrack(ntracks), ntracks);
  254. TpcLheCMTrack *new_track =
  255. (TpcLheCMTrack *)fCMTracks->At(ntracks);
  256. // fNTracks++;
  257. new_track->AddPoint(seed_hit, kFALSE);
  258. new_track->AddPoint(hit, kFALSE);
  259. seed_hit->SetUsage(kTRUE);
  260. hit->SetUsage(kTRUE);
  261. TrackExtension(new_track);
  262. } // if (!hit->GetUsage())
  263. } // over entries in segment
  264. } // over segments
  265. } // for (phi)
  266. } // for (theta)
  267. }
  268. //________________________________________________________________
  269. void TpcLheTrackFinder::
  270. GetClosestHit(TpcLheCMTrack *cmtrack, Int_t cur_station, Bool_t back) {
  271. //--- Returns the nearest hit in the next station
  272. TpcLheCMPoint *hit = NULL;
  273. TArrayI *philist = new TArrayI();
  274. TArrayI *thetalist = new TArrayI();
  275. if (cmtrack->GetNumberOfHits() < 3 ) {
  276. GetPhiRange(philist, cmtrack);
  277. GetThetaRange(thetalist, cmtrack);
  278. }
  279. Int_t* jphi = philist->GetArray();
  280. Int_t* jtheta = thetalist->GetArray();
  281. for (Int_t it = 0; it < thetalist->GetSize(); it++) {
  282. Int_t theta_segm = jtheta[it];
  283. for(Int_t ip = 0; ip < philist->GetSize(); ip++) {
  284. Int_t phi_segm = jphi[ip];
  285. // loop over entries in one sub segment
  286. TObjArray *cell = (TObjArray *)fSegments->At
  287. (fSegment->GetSegm(theta_segm, phi_segm));
  288. Int_t entries = cell->GetEntriesFast();
  289. if (entries) {
  290. for (Int_t sub_hit_num = 0; sub_hit_num < entries; sub_hit_num++) {
  291. hit = (TpcLheCMPoint *)cell->At(sub_hit_num);
  292. if (!hit->GetUsage()) { // hit was not used before
  293. if(fTrackCuts->VerifyTrack(cmtrack, hit, back)) {
  294. cmtrack->AddPoint(hit, kFALSE);
  295. hit->SetUsage(kTRUE);
  296. GetClosestHit(cmtrack, 0, back);
  297. }
  298. else {
  299. }
  300. } // if (!hit->GetUsage())
  301. } // over entries
  302. } // // over segments
  303. else {
  304. }
  305. } // for (phi)
  306. } // for (theta)
  307. delete philist;
  308. delete thetalist;
  309. }
  310. //______________________________________________________________________
  311. void TpcLheTrackFinder::GetThetaRange(TArrayI* list,TpcLheCMTrack *cmtrack) {
  312. //---
  313. TObjArray *trackpoint = cmtrack->GetCMHits();
  314. TpcLheCMPoint *hit = (TpcLheCMPoint *)trackpoint->First();
  315. Double_t thet = hit->GetTheta();
  316. Int_t theta_high, theta_low;
  317. if (fSegment->GetThetaSegm(thet) < fSegment->GetNumThetaSegments()/2) {
  318. theta_high = fSegment->GetThetaSegm(thet);
  319. theta_low = 0;
  320. } else {
  321. theta_high = fSegment->GetNumThetaSegments() - 1;
  322. theta_low = fSegment->GetThetaSegm(thet);
  323. }
  324. Int_t nTheta = theta_high - theta_low;
  325. list->Set(nTheta);
  326. for (Int_t s = 0; s < nTheta; s++) {
  327. list->AddAt(theta_low + s, s);
  328. }
  329. }
  330. //______________________________________________________________________
  331. void TpcLheTrackFinder::GetPhiRange(TArrayI* cl1, TpcLheCMTrack *cmtrack ) {
  332. //---
  333. TObjArray *trackpoint = cmtrack->GetCMHits();
  334. TpcLheCMPoint *hit = (TpcLheCMPoint *)trackpoint->First();
  335. TpcLheCMPoint *hit1 = (TpcLheCMPoint *)trackpoint->Last();
  336. Double_t phi = hit->GetPhi();
  337. Double_t phi_c = fTrackCuts->GetPhiPrediction(cmtrack);
  338. Double_t phi_m = phi_c - TMath::Pi()/2.;
  339. Double_t phi_p = phi_c + TMath::Pi()/2.;
  340. Double_t OutRadTpc = .42; // m
  341. Double_t Rad = cmtrack->GetRadius()/100.; // m
  342. Double_t MagField = 2; // Tesla
  343. Double_t pt = .2998 * Rad * MagField; // PR(pt);
  344. TpcLhePoint vert = cmtrack->GetVertex();
  345. Double_t xtr_c = vert.GetX()/100.; // m
  346. Double_t ytr_c = vert.GetY()/100.; // m
  347. Double_t xcross_1 = 0.;
  348. Double_t ycross_1 = 0.;
  349. Double_t xcross_2 = 0.;
  350. Double_t ycross_2 = 0.;
  351. Int_t crossing = fTrackCuts->Circle_Circle_Intersection(0., 0., OutRadTpc,
  352. xtr_c, ytr_c, Rad, &xcross_1, &ycross_1, &xcross_2, &ycross_2);
  353. Int_t range, phi_low, phi_high;
  354. Double_t phi_cr;
  355. // distance from 1st point to line (vertex - 2nd point)
  356. Double_t dist = hit->GetX()*hit1->GetY() - hit->GetY()*hit1->GetX();
  357. if (dist > 0.) {
  358. // negative track -- counterclockwize
  359. cmtrack->SetCharge(-1);
  360. Double_t phi_r = phi_c - TMath::Pi()/2.;
  361. if (phi_r < 0.) phi_r = TMath::TwoPi() + phi_r;
  362. phi_low = fSegment->GetPhiSegm(phi_r);
  363. if(crossing) {
  364. TVector2 v(xcross_2, ycross_2);
  365. phi_cr = (v.Phi() >= 0.) ? phi_cr = v.Phi() :
  366. phi_cr = v.Phi() + 2.*TMath::Pi();
  367. phi_high = fSegment->GetPhiSegm(phi_cr);
  368. }
  369. else {
  370. phi_high = fSegment->GetPhiSegm(phi_c);
  371. }
  372. range = phi_high - phi_low + 1;
  373. if (phi_high < phi_low) {
  374. range = fSegment->GetNumPhiSegments() - phi_low + phi_high + 1;
  375. }
  376. cl1->Set(range);
  377. for (Int_t s = 0; s < range; s++) {
  378. // cout << "range " << range << " cont " << phi_low + s << endl;
  379. Int_t segm = phi_low + s;
  380. if(segm >= fSegment->GetNumPhiSegments())
  381. segm = segm - fSegment->GetNumPhiSegments();
  382. cl1->AddAt(segm, s);
  383. }
  384. } else {
  385. cmtrack->SetCharge(1);
  386. TVector2 v(xcross_1, ycross_1);
  387. phi_cr = (v.Phi() >= 0.) ? phi_cr = v.Phi() :
  388. phi_cr = v.Phi() + 2.*TMath::Pi();
  389. phi_low = fSegment->GetPhiSegm(phi_cr);
  390. phi_high = fSegment->GetPhiSegm(phi_c + TMath::Pi()/2.);
  391. if(crossing) {
  392. phi_low = fSegment->GetPhiSegm(phi_cr);
  393. }
  394. else {
  395. phi_low = fSegment->GetPhiSegm(phi_c);
  396. }
  397. range = phi_high - phi_low + 1;
  398. if (phi_low > phi_high) {
  399. range = fSegment->GetNumPhiSegments() - phi_low + phi_high + 1;
  400. }
  401. cl1->Set(range);
  402. for (Int_t s = 0; s < range; s++) {
  403. Int_t segm = phi_high - s;
  404. if (segm < 0) segm = fSegment->GetNumPhiSegments() + segm;
  405. cl1->AddAt(segm, s);
  406. }
  407. }
  408. }
  409. //______________________________________________________________________
  410. void TpcLheTrackFinder::CopyClones(TClonesArray* cl1,
  411. TClonesArray* cl2) {
  412. //---
  413. Int_t nEntries = cl1->GetEntriesFast();
  414. TClonesArray& clref = *cl2;
  415. TpcLheTrack* old = NULL;
  416. for (Int_t i=0; i < nEntries; i++) {
  417. old = (TpcLheTrack*) cl1->At(i);
  418. new (clref[i]) TpcLheTrack(*old);
  419. }
  420. cout << " - " << cl2->GetEntriesFast() << " Geant tracks are copied.\n";
  421. } //&:)=
  422. /* ============================================================================
  423. ===============================================================================
  424. ===============================================================================
  425. ===============================================================================
  426. ===============================================================================
  427. */