TpcLheHitsMaker.cxx 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477
  1. #include "TpcLheHitsMaker.h"
  2. #include "TpcLheHit.h"
  3. #include "TpcLheTrack.h"
  4. #include "FairRunAna.h"
  5. #include "TpcPoint.h"
  6. #include "MpdMCTrack.h"
  7. #include "FairMCApplication.h"
  8. #include "FairRootManager.h"
  9. #include "TGeoManager.h"
  10. #include "TGeoVolume.h"
  11. #include "FairVolume.h"
  12. #include "FairTask.h"
  13. #include "TParticle.h"
  14. #include "TRandom.h"
  15. #include "Riostream.h"
  16. #include "TH1F.h"
  17. #include <iomanip>
  18. //#include "lhe.h"
  19. using namespace std;
  20. //#define TST
  21. //#define TPC
  22. #define TRACKS
  23. #define GEANT
  24. #define MASSPROD
  25. //_________________________________________________________________
  26. TpcLheHitsMaker::TpcLheHitsMaker() {
  27. //---
  28. fLheHits = new TClonesArray("TpcLheHit");
  29. fGeantTracks = new TClonesArray("TpcLheTrack");
  30. }
  31. //_________________________________________________________________
  32. TpcLheHitsMaker::TpcLheHitsMaker(const char *name,
  33. const char *title):FairTask(name) {
  34. //---
  35. fLheHits = new TClonesArray("TpcLheHit");
  36. fGeantTracks = new TClonesArray("TpcLheTrack");
  37. }
  38. //_________________________________________________________________
  39. TpcLheHitsMaker::~TpcLheHitsMaker() {
  40. //
  41. FairRootManager *fManager =FairRootManager::Instance();
  42. fManager->Write();
  43. delete fGeantTracks;
  44. }
  45. //_________________________________________________________________
  46. InitStatus TpcLheHitsMaker::ReInit() {
  47. // --- Update the pointers after reinitialisation
  48. // cout << "TpcLheHitsMaker::ReInit() \n";
  49. #if 0
  50. FairRunAna* ana = FairRunAna::Instance();
  51. FairRuntimeDb* rtdb=ana->GetRuntimeDb();
  52. fGeoPar=(FairGeoStsPar*)(rtdb->getContainer("FairGeoStsPar"));
  53. #endif
  54. return kSUCCESS;
  55. }
  56. //_________________________________________________________________
  57. InitStatus TpcLheHitsMaker::Init() {
  58. // cout << "InitStatus TpcLheHitsMaker::Init\n\n";
  59. fCuts = TpcLheTrackCuts::Instance();
  60. FairRootManager *fManager =FairRootManager::Instance();
  61. fTpcPoints = (TClonesArray *)fManager->GetObject("TpcPoint");
  62. if ( ! fTpcPoints ) {
  63. cout << "-I- "<< GetName() << "::Init: No TpcPoint array!" << endl;
  64. return kERROR;
  65. }
  66. #if 0
  67. fTstPoints = (TClonesArray*) fManager->GetObject("TSTPoint");
  68. if ( ! fTstPoints ) {
  69. cout << "-I- "<< GetName() << "::Init: No TSTPoint array!" << endl;
  70. return kERROR;
  71. }
  72. #endif
  73. #ifdef GEANT
  74. fListMCtracks = (TClonesArray *)fManager->GetObject("MCTrack");
  75. #endif
  76. fNofEvents = 0;
  77. Register();
  78. return kSUCCESS;
  79. }
  80. //_________________________________________________________________
  81. void TpcLheHitsMaker::GetTstHits() {
  82. //---
  83. #ifdef TST
  84. cout << " TpcLheHitsMaker::GetTstHits() : TST Hit entries " <<
  85. fTstPoints->GetEntriesFast() << endl;
  86. #endif
  87. TpcLheHit *hit=NULL;
  88. #if 0
  89. for (int j=0; j < fTstPoints->GetEntries(); j++ ) {
  90. FairTstPoint *point = (FairTstPoint*) fTstPoints->At(j);
  91. #ifdef KTST
  92. cout << " hit: " << j << " detID " << point->GetStationNr() <<
  93. " event: " << info->GetEventNumber() << " fake: " <<
  94. info->IsFake() << "\n ";
  95. #endif
  96. hit = AddHit();
  97. hit->SetHitNumber(fNHit++);
  98. hit->SetX(point->GetX());
  99. hit->SetY(point->GetY());
  100. hit->SetZ(point->GetZ());
  101. hit->SetTrackID(point->GetTrackID());
  102. hit->SetRefIndex(j);
  103. #if 0
  104. hit->SetXerr(point->GetDx());
  105. hit->SetYerr(point->GetDy());
  106. hit->SetZerr(point->GetDz());
  107. hit->SetStation(point->GetStationNr()); //
  108. if (info->GetEventNumber() == 0 ) {
  109. //&& !info->IsFake()) {
  110. //FairMCPoint *mcpoint = (FairMCPoint *) fTpcPoints->At(point->GetRefIndex());
  111. // hit->SetTrackID(mcpoint->GetTrackID());
  112. hit->SetTrackID(info->GetTrackId());
  113. hit->SetRefIndex(point->GetRefIndex());
  114. }
  115. else {
  116. hit->SetTrackID(-1);
  117. }
  118. #endif
  119. #ifdef TST
  120. hit->Print();
  121. #endif
  122. }
  123. #endif
  124. }
  125. //_________________________________________________________________
  126. void TpcLheHitsMaker::GetTpcHits() {
  127. #ifdef TPC
  128. cout << " TpcLheHitsMaker::GetTpcHits(): Tpc points entries " <<
  129. fTpcPoints->GetEntriesFast() <<endl;
  130. #endif
  131. Int_t nTpcPoints = fTpcPoints->GetEntriesFast();
  132. for (int j = 0; j < nTpcPoints; ++j) {
  133. TpcPoint* point = (TpcPoint*) fTpcPoints->UncheckedAt(j);
  134. #ifdef TPC
  135. point->Print(" "); PR(point->GetTrackID());
  136. #endif
  137. TpcLheHit* hit = AddHit();
  138. hit->SetHitNumber(fNHit++);
  139. hit->SetX(point->GetX());
  140. hit->SetY(point->GetY());
  141. hit->SetZ(point->GetZ());
  142. hit->SetXerr(0.03);
  143. hit->SetYerr(0.03);
  144. hit->SetZerr(0.03);
  145. hit->SetTrackID(point->GetTrackID());
  146. hit->SetRefIndex(j);
  147. //AZ
  148. Double_t dRPhi = 0, dZ = 0;
  149. //gRandom->Rannor(dRPhi,dZ); //AZ
  150. hit->SetZerr(0.1); // 1 mm error in Z
  151. hit->SetZ(hit->GetZ()+dZ*hit->GetZerr()); //add error
  152. TVector3 pos;
  153. point->Position(pos);
  154. hit->SetR(pos.Pt());
  155. hit->SetRphi(pos.Phi()*hit->GetR()+dRPhi*hit->GetRphiErr()); //add error
  156. Double_t phi = hit->GetRphi() / hit->GetR();
  157. hit->SetX(hit->GetR()*TMath::Cos(phi));
  158. hit->SetY(hit->GetR()*TMath::Sin(phi));
  159. hit->SetEdep(point->GetEnergyLoss());
  160. hit->SetStep(point->GetStep());
  161. //AZ
  162. #ifdef TPC
  163. hit->Print();
  164. #endif
  165. } // for (int j=0; j < fStsHits->GetEntries(); j++ )
  166. }
  167. //_________________________________________________________________
  168. void TpcLheHitsMaker::Exec(Option_t * option) {
  169. cout << "\n\n *** Event # " << ++fNofEvents << endl;
  170. cout << " ===== TpcLheHitsMaker =====\n";
  171. Reset();
  172. // GetTstHits();
  173. GetTpcHits();
  174. cout << " Total number of hits for tracking: " <<
  175. setw(5) << fLheHits->GetEntries() << endl;
  176. #ifdef GEANT
  177. CheckTracks();
  178. #ifdef TRACK
  179. PrintTracks(0);
  180. #endif
  181. #endif
  182. #ifdef MASSPROD
  183. // fTpcPoints->Clear("C");
  184. // fTstPoints->Clear("C");
  185. #endif
  186. }
  187. //_________________________________________________________________
  188. TpcLheHit * TpcLheHitsMaker::AddHit() {
  189. // Creates a new hit in the TClonesArray.
  190. TClonesArray& hitRef = *fLheHits;
  191. Int_t size = hitRef.GetEntriesFast();
  192. return new(hitRef[size]) TpcLheHit();
  193. }
  194. //_________________________________________________________________
  195. void TpcLheHitsMaker::Register() {
  196. //---
  197. FairRootManager::Instance()->
  198. Register("LheGeantTrack","Lhe", fGeantTracks, kFALSE);
  199. // Register("LheGeantTrack","Lhe", fGeantTracks, kTRUE);
  200. FairRootManager::Instance()->
  201. Register("LheHit","Lhe", fLheHits, kFALSE);
  202. //Register("LheHit","Lhe", fLheHits, kTRUE);
  203. }
  204. //_________________________________________________________________
  205. void TpcLheHitsMaker::Finish() {
  206. // FairRootManager *fManager =FairRootManager::Instance();
  207. // fManager->Fill();
  208. cout << " TpcLheHitsMaker::Finish : Geant tracks " <<
  209. fGeantTracks->GetEntriesFast() << endl;
  210. fGeantTracks->Delete();
  211. fLheHits->Delete();
  212. }
  213. //_________________________________________________________________
  214. void TpcLheHitsMaker::Reset() {
  215. //---
  216. // cout << " TpcLheHitsMaker::Reset " << endl;
  217. fNHit = 0;
  218. fNTrack = 0;
  219. if (fLheHits->GetEntries() != 0) fLheHits->Clear("C");
  220. if (fGeantTracks->GetEntries() != 0)fGeantTracks->Clear("C");
  221. }
  222. #ifdef GEANT
  223. //_______________________________________________________________
  224. void TpcLheHitsMaker::CheckTracks() {
  225. using namespace std;
  226. Int_t parent, tr_num, tr_pdg;
  227. MpdMCTrack *gtrack = 0;
  228. //AZ Int_t nMCtracks = fListMCtracks->GetEntries();
  229. Int_t nMCtracks = fListMCtracks->GetEntriesFast();
  230. Int_t selected_tracks = 0;
  231. //AZ for (Int_t ih = 0; ih < fLheHits->GetEntries(); ih++) {
  232. for (Int_t ih = 0; ih < fLheHits->GetEntriesFast(); ih++) {
  233. TpcLheHit *hit = (TpcLheHit *)fLheHits->UncheckedAt(ih);
  234. SetTrack(hit);
  235. // hit->Print();
  236. }
  237. //#ifdef TRACK
  238. cout << " MC tracks in event: " <<fListMCtracks->GetEntries();
  239. cout << " tracks in TPC " << fNTrack << endl;
  240. //#endif
  241. for (Int_t itrack = 0; itrack < fNTrack; itrack++) {
  242. TpcLheTrack *track = (TpcLheTrack*)fGeantTracks->UncheckedAt(itrack);
  243. tr_num = track->GetTrackNumber();
  244. if(tr_num >= 0 && tr_num < nMCtracks) {
  245. tr_pdg = 0;
  246. parent = 99999;
  247. gtrack = (MpdMCTrack *) fListMCtracks->At(tr_num);
  248. tr_pdg = gtrack->GetPdgCode();
  249. track->SetPid(tr_pdg);
  250. track->SetCharge(TMath::Sign(1,tr_pdg));
  251. parent = gtrack->GetMotherId();
  252. if(parent == -1) {
  253. track->ComesFromMainVertex(kTRUE);
  254. track->SetGood(kTRUE);
  255. }
  256. TVector3 vertex;
  257. gtrack->GetStartVertex(vertex);
  258. track->SetVertex(vertex.X(), vertex.Y(), vertex.Z());
  259. TVector3 moment;
  260. gtrack->GetMomentum(moment);
  261. track->SetPx(moment.X());
  262. track->SetPy(moment.Y());
  263. track->SetPz(moment.Z());
  264. //Int_t nhit = track->GetNumberOfHits();
  265. //TObjArray *chits = track->GetRHits();
  266. // print some info about track
  267. /*
  268. cout << "\n *** track *** " << track->GetTrackNumber() <<
  269. " pdg " << track->GetPid() <<
  270. " nhits " <<nhit << " p_t " << track->GetPt() << endl;
  271. track->Print();
  272. TpcLheHit *hit = (TpcLheHit *)chits->First();
  273. hit->Print();
  274. */
  275. if (fCuts->IsGoodGeantTrack(track)) {
  276. selected_tracks++;
  277. }
  278. if (parent == 99999 ) cout << " track parent unknown";
  279. }
  280. }
  281. #if 0
  282. for (Int_t ih = 0; ih < fLheHits->GetEntries(); ih++) {
  283. TpcLheHit *hit = (TpcLheHit *)fLheHits->UncheckedAt(ih);
  284. Int_t mTrackNumber = hit->GetTrackID();
  285. for (Int_t itrack = 0; itrack < fNTrack; itrack++) {
  286. TpcLheTrack *track = (TpcLheTrack*)fGeantTracks->UncheckedAt(itrack);
  287. Int_t tr_num = track->GetTrackNumber();
  288. if (mTrackNumber == tr_num && !track->IsGood()) {
  289. hit->SetUsage(kTRUE);
  290. }
  291. }
  292. }
  293. #endif
  294. cout << "Total number of tracks in TPC: " <<
  295. setw(5) << fGeantTracks->GetEntries() << endl;
  296. cout << " Good tracks in TPC: " <<
  297. setw(5) << selected_tracks << endl;
  298. }
  299. //________________________________________________________________
  300. void TpcLheHitsMaker::SetTrack(TpcLheHit *hit) {
  301. Int_t mTrackNumber = hit->GetTrackID();
  302. Int_t mTN;
  303. #ifdef HITS
  304. cout << " fNHit = " << hit->GetHitNumber();
  305. cout << " fNTrack = " << fNTrack;
  306. cout << " mTrackNumber = " << mTrackNumber;
  307. cout << endl;
  308. #endif
  309. TpcLheTrack *track = NULL;
  310. Bool_t newtrack = kTRUE;
  311. if (fNTrack > 0) {
  312. for (Int_t i = 0; i < fGeantTracks->GetEntriesFast(); i++) {
  313. TpcLheTrack *currtrk = ((TpcLheTrack*) fGeantTracks->At(i));
  314. mTN = currtrk->GetTrackNumber();
  315. // PR( currtrk->GetTrackNumber());
  316. if (mTN == mTrackNumber) {
  317. track = currtrk;
  318. newtrack = kFALSE;
  319. break;
  320. }
  321. }
  322. }
  323. if (newtrack) track = AddTrack(mTrackNumber);
  324. track->AddHit(hit);
  325. }
  326. //________________________________________________________________
  327. TpcLheTrack *TpcLheHitsMaker::AddTrack(Int_t mTrackNumber) {
  328. // Add a new track to the list of tracks for this event.
  329. TClonesArray &tracks = *fGeantTracks;
  330. TpcLheTrack *track = new(tracks[fNTrack++]) TpcLheTrack(mTrackNumber);
  331. return track;
  332. }
  333. //_________________________________________________________________
  334. void TpcLheHitsMaker::PrintTracks(Int_t ntr) {
  335. // cout << " " << fNTrack << " tracks \n\n ";
  336. // cout << "Station X Y Z Px Py Pz \n\n";
  337. Int_t ptracks;
  338. if(ntr != 0 )
  339. ptracks = ntr;
  340. else
  341. ptracks = fGeantTracks->GetEntries();
  342. for (Int_t itrack = 0; itrack < ptracks; itrack++) {
  343. TpcLheTrack *track = (TpcLheTrack *)fGeantTracks->UncheckedAt(itrack);
  344. track->Print(); //
  345. }
  346. }
  347. #endif
  348. ClassImp(TpcLheHitsMaker)