MpdZdcTstSim.cxx 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474
  1. /*************************************************************************************
  2. *
  3. * Class MpdZdcTstSim
  4. *
  5. * Author: Elena Litvinenko
  6. * e-mail: litvin@nf.jinr.ru
  7. * Version: 8-Apr-2008
  8. *
  9. ************************************************************************************/
  10. #include "MpdZdcTstSim.h"
  11. #include "MpdZdcPoint.h"
  12. #include "TpcPoint.h"
  13. #include "FairRootManager.h"
  14. #include "MpdMCTrack.h"
  15. #include "FairRunAna.h"
  16. #include "FairRuntimeDb.h"
  17. #include "FairBaseParSet.h"
  18. using std::cout;
  19. using std::endl;
  20. MpdZdcTstSim::MpdZdcTstSim()
  21. // : FairTask ("ZdcTstSim")
  22. {
  23. fVerbose = 0;
  24. nevents = 0;
  25. fTree = 0;
  26. fH3 = 0;
  27. }
  28. MpdZdcTstSim::MpdZdcTstSim(const char *name, const char *title, Int_t verbose)
  29. : FairTask (name)
  30. {
  31. if (verbose)
  32. cout << "-Info- MpdZdcTstSim::MpdZdcTstSim begin" << endl;
  33. fVerbose = verbose;
  34. nevents = 0;
  35. fTree = 0;
  36. fTreeSummary = 0;
  37. fH3 = 0;
  38. // CreateMyTree();
  39. if (fVerbose)
  40. cout << "-Info- MpdZdcTstSim::MpdZdcTstSim end" << endl;
  41. }
  42. MpdZdcTstSim::MpdZdcTstSim(const char *name, const char *title, Int_t verbose, Int_t flag)
  43. {
  44. if (verbose)
  45. cout << "-Info- MpdZdcTstSim::MpdZdcTstSim begin" << endl;
  46. fVerbose = verbose;
  47. nevents = 0;
  48. fTree = 0;
  49. fTreeSummary = 0;
  50. // CreateMyTree();
  51. if (fVerbose)
  52. cout << "-Info- MpdZdcTstSim::MpdZdcTstSim Init starting" << endl;
  53. Init();
  54. if (fVerbose)
  55. cout << "-Info- MpdZdcTstSim::MpdZdcTstSim end" << endl;
  56. }
  57. MpdZdcTstSim::~MpdZdcTstSim()
  58. {
  59. if (fTree)
  60. fTree->Write();
  61. if (fTreeSummary)
  62. fTreeSummary->Write();
  63. if (fH3)
  64. fH3->Write();
  65. }
  66. void MpdZdcTstSim::CreateMyTree()
  67. {
  68. if (!fTree)
  69. fTree = new TNtuple ("zdcpoints","zdc points",
  70. "tpdg:tmotherid:tpx:tpy:tpz:tx:ty:tz:tmass:tenergy:tpt:tp:trapidity:tstarttime:ptrackid:peventid:px:py:pz:ppx:ppy:ppz:ptime:plength:peloss:pmcvolumeid:pcopy:pcopymother:tmotherpdg:tmotherz:tZ:tA",
  71. 1000000);
  72. if (!fTreeSummary)
  73. fTreeSummary = new TNtupleD ("zdcsummary","zdc summary",
  74. "elossall:elosszdc1:elosszdc2:tpctracks:b:ekin:ekin1:nucleons1:pions1:protons1:ekin2:nucleons2:pions2:protons2",
  75. // "elossall:elosszdc1:elosszdc2:tpctracks:b:ekin:ezdc1[4]:ezdc2[4]",
  76. 1000000);
  77. // if (fVerbose) {
  78. // cout << "-Info- MpdZdcTstSim::CreateMyTree zdcpoints:" << (fTree->GetCurrentFile()->GetName())<< endl;
  79. // cout << "-Info- MpdZdcTstSim::CreateMyTree zdcsummary:" << (fTreeSummary->GetCurrentFile()->GetName())<< endl;
  80. // }
  81. if (!fH3) {
  82. fH3 = new TH3F ("H3","Counts(ELoss,TpcTracks,b)",40,0,40,120,0,1200,32,0,16);
  83. // fH3 = new TH3F ("H3","Counts(ELoss,TpcTracks,b)",50,0,10,120,0,1200,40,0,2);
  84. // fH3->SetDirectory(0);
  85. }
  86. }
  87. InitStatus MpdZdcTstSim::Init()
  88. {
  89. if (fVerbose)
  90. cout << "-Info- MpdZdcTstSim::Init begin" << endl;
  91. FairRootManager* ioman = FairRootManager::Instance();
  92. if (! ioman) {
  93. cout << "-E- MpdZdcTstSim::Init: "
  94. << "RootManager not instantised!" << endl;
  95. return kERROR;
  96. }
  97. fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
  98. if ( ! fMCTrackArray) {
  99. cout << "-E- MpdZdcTstSim::Init: No MCTrack array!"
  100. << endl;
  101. return kERROR;
  102. }
  103. fMCZdcPointArray = (TClonesArray*) ioman->GetObject("ZdcPoint");
  104. if ( !fMCZdcPointArray) {
  105. cout << "-E- MpdZdcTstSim::Init: No ZdcPoint array!"
  106. << endl;
  107. return kERROR;
  108. }
  109. fMCTpcPointArray = (TClonesArray*) ioman->GetObject("TpcPoint");
  110. if ( !fMCTpcPointArray) {
  111. cout << "-E- MpdZdcTstSim::Init: No TpcPoint array!"
  112. << endl;
  113. // return kERROR;
  114. }
  115. fMCEventHeader = (FairMCEventHeader*) ioman->GetObject("MCEventHeader.");
  116. if ( !fMCEventHeader) {
  117. cout << "-E- MpdZdcTstSim::Init: No MCEventHeader array!"
  118. << endl;
  119. // return kERROR;
  120. }
  121. if (fVerbose)
  122. cout << "-Info- MpdZdcTstSim::Init end" << endl;
  123. CreateMyTree();
  124. fTree->SetDirectory(ioman->GetOutFile());
  125. fTreeSummary->SetDirectory(ioman->GetOutFile());
  126. if (fH3)
  127. fH3->SetDirectory(ioman->GetOutFile());
  128. ioman->Register("zdcpoints","ZDC", fTree, kTRUE);
  129. ioman->Register("zdcsummary","ZDC", fTreeSummary, kTRUE);
  130. if (fH3)
  131. ioman->Register("H3","ZDC", fH3, kTRUE);
  132. if (fVerbose) {
  133. cout << "-Info- MpdZdcTstSim::Init zdcpoints: " << (fTree->GetCurrentFile()->GetName())<< endl;
  134. cout << "-Info- MpdZdcTstSim::Init zdcsummary: " << (fTreeSummary->GetCurrentFile()->GetName())<< endl;
  135. }
  136. return kSUCCESS;
  137. }
  138. InitStatus MpdZdcTstSim::ReInit()
  139. {
  140. // // SetParContainers();
  141. // FairRunAna* ana = FairRunAna::Instance();
  142. // FairRuntimeDb* rtdb=ana->GetRuntimeDb();
  143. // MpdZdcGeoPar *fPar=(MpdZdcGeoPar*)(rtdb->getContainer("MpdZdcGeoPar"));
  144. return kSUCCESS;
  145. }
  146. void MpdZdcTstSim::Exec(Option_t * option)
  147. {
  148. if (fVerbose>3)
  149. cout << "-Info- MpdZdcTstSim::Exec begin" << endl;
  150. Double_t elossall=0,elosszdc1=0,elosszdc2=0,e=0,ezdc1[4],ezdc2[4];
  151. Int_t nPoints, nMCTracks;
  152. if (fTree) {
  153. TVector3 position, momentum, vertex;
  154. MpdZdcPoint *pPoint=NULL;
  155. MpdMCTrack *pTrack = NULL;
  156. Int_t pdg=0;
  157. Float_t fArgs[30];
  158. for (int iii=0;iii<4;iii++) {
  159. ezdc1[iii]=0; ezdc2[iii]=0;
  160. }
  161. nPoints = fMCZdcPointArray->GetEntriesFast();
  162. nMCTracks = fMCTrackArray->GetEntriesFast();
  163. if (fVerbose) {
  164. cout << " MpdZdcPoint: " << nPoints << " entries" << endl;
  165. cout << " MCTrack: " << nMCTracks << " entries" << endl;
  166. }
  167. for (Int_t i=0; i<nPoints; i++) {
  168. // if (fVerbose>2)
  169. // cout << " .";
  170. pPoint = (MpdZdcPoint*)fMCZdcPointArray->UncheckedAt(i);
  171. if (!pPoint) {
  172. if (fVerbose>2) {
  173. cout << " -";
  174. }
  175. continue;
  176. }
  177. if (pPoint->GetTrackID() == -2) {
  178. if (fVerbose>2) {
  179. // cout << i << ":" <<(pPoint->GetTrackID()) << "!";
  180. cout << "MpdZdcTstSim::Exec: MCTrack not saved: particle number:" << i << " must be too slow !" <<endl;
  181. }
  182. continue;
  183. }
  184. pTrack = (MpdMCTrack*)fMCTrackArray->UncheckedAt(pPoint->GetTrackID());
  185. if (!pTrack){
  186. if (fVerbose>2) {
  187. cout << i << ":" <<(pPoint->GetTrackID()) << "!";
  188. }
  189. continue;
  190. }
  191. try
  192. {
  193. // if ((fVerbose>2)&&(nMCTracks==174937)) {
  194. // cout << i << ":" <<(pPoint->GetTrackID()) << " track addr:" << pTrack<< "! ";
  195. // }
  196. pdg = pTrack->GetPdgCode();
  197. fArgs[0] = pdg;
  198. }
  199. catch (...)
  200. {
  201. cout << "MpdZdcTstSim::Exec: Exception 1: particle number:" << i << " track number:" <<
  202. (pPoint->GetTrackID()) << " track addr:" <<pTrack << endl;
  203. }
  204. fArgs[1] = pTrack->GetMotherId();
  205. pTrack->GetMomentum(momentum);
  206. fArgs[2] = momentum.X();
  207. fArgs[3] = momentum.Y();
  208. fArgs[4] = momentum.Z();
  209. pTrack->GetStartVertex(vertex);
  210. fArgs[5] = vertex.X();
  211. fArgs[6] = vertex.Y();
  212. fArgs[7] = vertex.Z();
  213. Int_t tZ=-1000,tA=-1000;
  214. TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle((Int_t)floor(fArgs[0]));
  215. if ( particle ) {
  216. fArgs[8] = particle->Mass();
  217. tZ=TMath::Nint(particle->Charge() / 3.);
  218. if (pdg>1e9) {
  219. tA=((pdg%10000)/10);
  220. if (!tZ)
  221. tZ=TMath::Nint((pdg-1e9)/10000);
  222. }
  223. }
  224. else {
  225. fArgs[8] = 0;
  226. if (pdg>1e9) {
  227. tZ=TMath::Nint((pdg-1e9)/10000);
  228. tA=((pdg%10000)/10);
  229. }
  230. }
  231. fArgs[9] = TMath::Sqrt(fArgs[2]*fArgs[2] +
  232. fArgs[3]*fArgs[3] + fArgs[4]*fArgs[4] + fArgs[8]*fArgs[8]) ;
  233. fArgs[10] = TMath::Sqrt(fArgs[2]*fArgs[2]+fArgs[3]*fArgs[3]);
  234. fArgs[11] = TMath::Sqrt(fArgs[2]*fArgs[2]+fArgs[3]*fArgs[3]+fArgs[4]*fArgs[4]);
  235. fArgs[12]= 0.5*TMath::Log((fArgs[9]+fArgs[4])/(fArgs[9]-fArgs[4]));
  236. fArgs[13]=1; // not implemented on the level of MpdMCTrack now
  237. //fArgs[13]=pTrack->GetStartTime();
  238. fArgs[14]= pPoint->GetTrackID();
  239. // fArgs[15]= pPoint->GetEventID();
  240. fArgs[15]= nevents;
  241. fArgs[16]= pPoint->GetX();
  242. fArgs[17]= pPoint->GetY();
  243. fArgs[18]= pPoint->GetZ();
  244. fArgs[19]= pPoint->GetPx();
  245. fArgs[20]= pPoint->GetPy();
  246. fArgs[21]= pPoint->GetPz();
  247. fArgs[22]= pPoint->GetTime();
  248. fArgs[23]= pPoint->GetLength();
  249. fArgs[24]= pPoint->GetEnergyLoss();
  250. // fArgs[25]= pPoint->GetModule();
  251. // fArgs[26]= pPoint->GetRow();
  252. // fArgs[27]= pPoint->GetCopy();
  253. fArgs[25]= pPoint->GetDetectorID();
  254. fArgs[26]= pPoint->GetCopy();
  255. fArgs[27]= pPoint->GetCopyMother();
  256. if (fArgs[1]!=-1) {
  257. pTrack = (MpdMCTrack*)fMCTrackArray->UncheckedAt(pTrack->GetMotherId());
  258. if (pTrack) {
  259. fArgs[28]= pTrack->GetPdgCode();
  260. TVector3 vertexTemp;
  261. pTrack->GetStartVertex(vertexTemp);
  262. fArgs[29]= vertexTemp.Z();
  263. }
  264. else {
  265. fArgs[28]= -2000000;
  266. fArgs[29]= -2000000;
  267. }
  268. }
  269. else {
  270. fArgs[28]= -1000000;
  271. fArgs[29]= -1000000;
  272. }
  273. //gertsen COMMENT BECAUSE OF OUT OF RANGE!!!
  274. //fArgs[30]=tZ;
  275. //fArgs[31]=tA;
  276. elossall += pPoint->GetEnergyLoss();
  277. Double_t ekin = fArgs[9] - fArgs[8];
  278. // if ((fArgs[25]<=30)&&(fArgs[1]==-1))
  279. if (fArgs[1]==-1) { // only primary particles:
  280. // if (fArgs[27]==1) // zdc1
  281. if (fArgs[18]>0) // zdc1
  282. elosszdc1 += pPoint->GetEnergyLoss();
  283. else
  284. elosszdc2 += pPoint->GetEnergyLoss();
  285. if (fArgs[26]==1) { // zdc entrances
  286. e += ekin; // kinetic energy of all primary particles at both zdc entrances
  287. // if (fArgs[27]==1) { // zdc1
  288. if (fArgs[18]>0) { // zdc1
  289. ezdc1[0] += ekin; // all
  290. switch (pdg) {
  291. case 2212: // protons
  292. ezdc1 [1] += ekin; // nucleons
  293. ezdc1 [3] += ekin; // protons
  294. break;
  295. case 2112: // neutrons
  296. ezdc1[1] += ekin; // nucleons
  297. break;
  298. case 211: // pions
  299. case -211:
  300. ezdc1 [2] += fArgs[9];
  301. break;
  302. default:
  303. break;
  304. }
  305. }
  306. else { // zdc2
  307. ezdc2[0] += ekin;
  308. switch (pdg) {
  309. case 2212: // protons
  310. ezdc2 [1] += ekin; // nucleons
  311. ezdc2 [3] += ekin; // protons
  312. break;
  313. case 2112: // neutrons
  314. ezdc2[1] += ekin; // nucleons
  315. break;
  316. case 211: // pions
  317. case -211:
  318. ezdc2 [2] += fArgs[9];
  319. break;
  320. default:
  321. break;
  322. }
  323. }
  324. } // end of zdc entrance
  325. fTree->Fill(fArgs);
  326. } // end of primary particles
  327. // fTree->Fill(fArgs);
  328. }
  329. // cout << " .finished!" << endl;
  330. nevents++;
  331. }
  332. Double_t b=0;
  333. if (fMCEventHeader) {
  334. b = ((FairMCEventHeader*)fMCEventHeader)->GetB();
  335. }
  336. if (fTreeSummary) {
  337. Int_t nTpcTracks =0;
  338. Int_t trackId;
  339. if (fMCTpcPointArray) {
  340. std::map <Int_t,Bool_t> tpc_tracks;
  341. tpc_tracks.clear();
  342. Int_t i_output=5;
  343. for (Int_t i=0;i<fMCTpcPointArray->GetEntriesFast();i++) {
  344. trackId = ((TpcPoint *)fMCTpcPointArray->UncheckedAt(i))->GetTrackID();
  345. if ((trackId>=0)&& (trackId<nMCTracks)) {
  346. if (((MpdMCTrack*)fMCTrackArray->UncheckedAt(trackId))->GetMotherId()==-1) {
  347. if (tpc_tracks.find(trackId)==tpc_tracks.end())
  348. tpc_tracks[trackId]=kTRUE;
  349. }
  350. }
  351. else {
  352. if (i_output>0) {
  353. cout << "-E- Wrong TpcPoint GetTrackID = " << trackId << " for i=" << i << " b=" << b << endl;
  354. i_output--;
  355. }
  356. }//else
  357. }//for (Int_t i=0;i<fMCTpcPointArray->GetEntriesFast();i++)
  358. nTpcTracks = tpc_tracks.size();
  359. }
  360. fTreeSummary->Fill(elossall,elosszdc1,elosszdc2,nTpcTracks,b,e,ezdc1[0],ezdc1[1],ezdc1[2],ezdc1[3],ezdc2[0],ezdc2[1],ezdc2[2],ezdc2[3]);
  361. // cout << " . " << (fTreeSummary) << endl;
  362. fH3->Fill(elossall,nTpcTracks,b);
  363. }
  364. else
  365. cout << " !";
  366. if (fVerbose>3)
  367. cout << "-Info- MpdZdcTstSim::Exec end" << endl;
  368. }
  369. void MpdZdcTstSim::Finish()
  370. {
  371. fMCTrackArray->Clear();
  372. fMCZdcPointArray->Clear();
  373. if (fMCTpcPointArray)
  374. fMCTpcPointArray->Clear();
  375. }
  376. ClassImp(MpdZdcTstSim)