MpdNDetAnalysis.cxx 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603
  1. /*
  2. First simple version of Nuetron Detector (NDet) analysis
  3. May 2009
  4. */
  5. #include "MpdNDetAnalysis.h"
  6. #include "MpdNDetPoint.h"
  7. #include "MpdNDetPointLite.h"
  8. //... #include "MpdNDetHitFastMC.h"
  9. //... #include "MpdNDetRecParticle.h"
  10. #include "FairRootManager.h"
  11. #include "FairMCTrack.h"
  12. #include "TParticle.h"
  13. #include "TH1F.h"
  14. #include "TH2F.h"
  15. #include "TVector3.h"
  16. #include "TLorentzVector.h"
  17. #include <iostream>
  18. #include <stdlib.h>
  19. using namespace std;
  20. //------------------------------------------------------------------------------
  21. MpdNDetAnalysis::MpdNDetAnalysis() :FairTask()
  22. {
  23. }
  24. //------------------------------------------------------------------------------
  25. MpdNDetAnalysis::MpdNDetAnalysis(const char *name, const char *title)
  26. :FairTask(name)
  27. {
  28. SetTitle(title);
  29. fEvent = 0;
  30. fDebug = "";
  31. Double_t pi = TMath::Pi();
  32. cout<<"Start ANALYSIS-------------------->"<<endl;
  33. fhNevents = new TH1F("hNevents" ,"Number of events",1, 0.5,1.5);
  34. // Neutron tracks histograms
  35. //Primary
  36. fhNeut0Pabs = new TH1F("hNeut0Pabs" ,"Momentum of prim. neutrons",100, 0.,10.);
  37. fhNeut0Pt = new TH1F("hNeut0Pt" ,"p_{T} of prim. neutrons" ,100, 0.,2.);
  38. fhNeut0Eta = new TH1F("hNeut0Eta" ,"#eta of prim. neutrons" , 50,-2., 2.);
  39. fhNeut0Theta= new TH1F("hNeut0Theta","#theta of prim. neutrons" ,100, 0., pi);
  40. fhNeut0Phi = new TH1F("hNeut0Phi" ,"#phi of prim. neutrons" ,100,-pi, pi);
  41. fhNeut0PtEta= new TH2F("hNeut0PtEta","p_{T},#eta of prim. neutrons",100, 0.,2., 50,-2., 2.);
  42. fhNeut0Rad = new TH1F("hNeut0Rad","r=#sqrt{x^{2}+y^{2}}",200,0.,200.);
  43. fhNeut0RadZ = new TH2F("hNeut0RadZ","n_{prim} r=#sqrt{x^{2}+y^{2}} vs Z",200,-500.,500.,200,0.,200.);
  44. //Secondary
  45. fhNeut1Pabs = new TH1F("hNeut1Pabs" ,"Momentum of sec. neutrons",100, 0.,10.);
  46. fhNeut1Pt = new TH1F("hNeut1Pt" ,"p_{T} of sec. neutrons" ,100, 0.,2.);
  47. fhNeut1Eta = new TH1F("hNeut1Eta" ,"#eta of sec. neutrons" , 50,-2., 2.);
  48. fhNeut1Theta= new TH1F("hNeut1Theta","#theta of sec. neutrons" ,100, 0., pi);
  49. fhNeut1Phi = new TH1F("hNeut1Phi" ,"#phi of sec. neutrons" ,100,-pi, pi);
  50. fhNeut1PtEta= new TH2F("hNeut1PtEta","p_{T},#eta of sec. neutrons",100, 0.,2., 50,-2., 2.);
  51. fhNeut1Rad = new TH1F("hNeut1Rad","r=#sqrt{x^{2}+y^{2}}",200,0.,200.);
  52. fhNeut1RadZ = new TH2F("hNeut1RadZ","n_{seco} r=#sqrt{x^{2}+y^{2}} vs Z",200,-500.,500.,200,0.,200.);
  53. fhNeut1RadE = new TH2F("hNeut1RadE","n_{seco} r=#sqrt{x^{2}+y^{2}} vs Ekin",200,0.,2.,200,0.,200.);
  54. // Not Neutoron tracks histograms
  55. //Protons prim
  56. fhProt0Pabs = new TH1F("hProt0Pabs" ,"Momentum of prim. protons",100, 0.,10.);
  57. fhProt0Pt = new TH1F("hProt0Pt" ,"p_{T} of prim. protons" ,100, 0.,2.);
  58. fhProt0Eta = new TH1F("hProt0Eta" ,"#eta of prim. protons" , 50,-2., 2.);
  59. fhProt0Theta= new TH1F("hProt0Theta","#theta of prim. protons" ,100, 0., pi);
  60. fhProt0Phi = new TH1F("hProt0Phi" ,"#phi of prim. protons" ,100,-pi, pi);
  61. fhProt0PtEta= new TH2F("hProt0PtEta","p_{T},#eta of prim. protons",100, 0.,2., 50,-2., 2.);
  62. fhProt0Rad = new TH1F("hProt0Rad","r=#sqrt{x^{2}+y^{2}}",200,0.,200.);
  63. //Protons sec
  64. fhProt1Pabs = new TH1F("hProt1Pabs" ,"Momentum of sec. protons",100, 0.,10.);
  65. fhProt1Pt = new TH1F("hProt1Pt" ,"p_{T} of sec. protons" ,100, 0.,2.);
  66. fhProt1Eta = new TH1F("hProt1Eta" ,"#eta of sec. protons" , 50,-2., 2.);
  67. fhProt1Theta= new TH1F("hProt1Theta","#theta of sec. protons" ,100, 0., pi);
  68. fhProt1Phi = new TH1F("hProt1Phi" ,"#phi of sec. protons" ,100,-pi, pi);
  69. fhProt1PtEta= new TH2F("hProt1PtEta","p_{T},#eta of sec. protons",100, 0.,2., 50,-2., 2.);
  70. fhProt1Rad = new TH1F("hProt1Rad","r=#sqrt{x^{2}+y^{2}}",200,0.,200.);
  71. // MC points histograms
  72. fhMcpPabs = new TH1F("hMcpPabs" ,"Momentum of detected particles",120, 0.,30.);
  73. fhMcpPt = new TH1F("hMcpPt" ,"p_{T} of detected particles" ,100, 0.,10.);
  74. fhMcpEta = new TH1F("hMcpEta" ,"#eta of detected particles" , 80,-2., 6.);
  75. fhMcpTheta= new TH1F("hMcpTheta","#theta of detected particles" ,100, 0., pi);
  76. fhMcpPhi = new TH1F("hMcpPhi" ,"#phi of detected particles" ,100,-pi, pi);
  77. fhMcpPtEta= new TH2F("hMcpPtEta","p_{T},#eta of detected particles",
  78. 100, 0.,10., 80,-2., 6.);
  79. fhMcpXY = new TH2F("hMcpXY","(X,Y) coordinates of MC points",
  80. 300,-600,600,300,-600,600);
  81. //MC pionts PDG code
  82. fhMcpPdgNdet = new TH1D("hMcpPdgNdet" ,"PDD code on the enter of Ndet",10001, -5000.,5000.);
  83. // Hits histograms
  84. fhHitXY = new TH2F("hHitXY","NDET hit (x,y)",100,-600.,600.,100,-500.,500.);
  85. fhHitE = new TH1F("hHitE" ,"NDET hit energy",120,0.,30.);
  86. // Reconstructed particles histograms
  87. fhRecPabs = new TH1F("hRecPabs","Reconstructed |p|" ,120, 0.,30.);
  88. fhRecPt = new TH1F("hRecPt" ,"Reconstructed p_{T}",100, 0.,10.);
  89. fhRecYrap = new TH1F("hRecYrap","Reconstructed y" , 80,-2., 6.);
  90. fhRecPhi = new TH1F("hRecPhi" ,"Reconstructed #phi" ,100,-pi, pi);
  91. // Acceptance
  92. fhPtYrap = new TH2F("PtYrap","pt vs y of primary particles", 100,0.,6.,100,0.,5.);
  93. fhPtYrapSignal = new TH2F("PtYrapSignal","pt vs y of signal particles", 100,0.,6.,100,0.,5.);
  94. fhPtYrapAccept = new TH2F("PtYrapAccept","NDET p_{T}-y acceptance",100,0.,6.,100,0.,5.);
  95. }
  96. //------------------------------------------------------------------------------
  97. void MpdNDetAnalysis::WriteOutput()
  98. {
  99. cout<<"Write histograms to MY FILE neutron.root"<<endl;
  100. ftest = new TFile("neutron.root","recreate");
  101. if(ftest) {
  102. fhNevents->Write();
  103. fhMcpPdgNdet->Write();
  104. fhNeut0Pabs ->Write();
  105. fhNeut0Pt ->Write();
  106. fhNeut0Eta ->Write();
  107. fhNeut0Theta->Write();
  108. fhNeut0Phi ->Write();
  109. fhNeut0PtEta->Write();
  110. fhNeut0Rad->Write();
  111. fhNeut0RadZ->Write();
  112. fhNeut1Pabs ->Write();
  113. fhNeut1Pt ->Write();
  114. fhNeut1Eta ->Write();
  115. fhNeut1Theta->Write();
  116. fhNeut1Phi ->Write();
  117. fhNeut1PtEta->Write();
  118. fhNeut1Rad->Write();
  119. fhNeut1RadZ->Write();
  120. fhNeut1RadE->Write();
  121. fhProt0Pabs ->Write();
  122. fhProt0Pt ->Write();
  123. fhProt0Eta ->Write();
  124. fhProt0Theta->Write();
  125. fhProt0Phi ->Write();
  126. fhProt0PtEta->Write();
  127. fhProt0Rad->Write();
  128. fhProt1Pabs ->Write();
  129. fhProt1Pt ->Write();
  130. fhProt1Eta ->Write();
  131. fhProt1Theta->Write();
  132. fhProt1Phi ->Write();
  133. fhProt1PtEta->Write();
  134. fhProt1Rad->Write();
  135. ftest->Close();
  136. //cout<<"finish writting..."<<endl;
  137. }
  138. }
  139. //------------------------------------------------------------------------------
  140. MpdNDetAnalysis::~MpdNDetAnalysis()
  141. {
  142. cout<<"Destruct the run<----------------"<<endl;
  143. }
  144. //------------------------------------------------------------------------------
  145. InitStatus MpdNDetAnalysis::Init()
  146. {
  147. // Activate branches with NDET objects
  148. FairRootManager *manager= FairRootManager::Instance();
  149. // all tracks
  150. // fMCtrack = (TClonesArray*)manager->ActivateBranch("MCTrack");
  151. fMCtrack = (TClonesArray*)manager->GetObject("MCTrack"); // EL
  152. //NDET MC points inside NDET
  153. // fNDetPointLite = (TClonesArray*)manager->ActivateBranch("NDetPointLite");
  154. fNDetPointLite = (TClonesArray*)manager->GetObject("NDetPointLite"); // EL
  155. //NDET MC points on entrance to NDET
  156. // fNDetPoint = (TClonesArray*)manager->ActivateBranch("NDetPoint");
  157. fNDetPoint = (TClonesArray*)manager->GetObject("NDetPoint"); // EL
  158. //NDET hits
  159. /*
  160. fListNDEThits = (TClonesArray*)manager->ActivateBranch("EcalHitFastMC");
  161. */
  162. //NDET reconstructed points
  163. /*
  164. fListNDETrp = (TClonesArray*)manager->ActivateBranch("EcalRecParticle");
  165. */
  166. return kSUCCESS;
  167. }
  168. //------------------------------------------------------------------------------
  169. void MpdNDetAnalysis::Exec(Option_t* option)
  170. {
  171. // make NDET analysis
  172. printf("NDET analysis: event %d\n",fEvent);
  173. OneEventAnalysis();
  174. //cout<<"option = "<<*option<<endl;
  175. //if(*option==-1)
  176. Finish();
  177. }
  178. //------------------------------------------------------------------------------
  179. void MpdNDetAnalysis::OneEventAnalysis()
  180. {
  181. //ReadMCTrack();
  182. Acceptance();
  183. /* AnalyzePrimary(); */
  184. /* AnalyzeMCPointsEdep(); */
  185. /* AnalyzeMCPointsWall(); */
  186. /* AnalyzeHits(); */
  187. /* AnalyzeRecParticles(); */
  188. fEvent++;fhNevents->Fill(1.);
  189. }
  190. //------------------------------------------------------------------------------
  191. void MpdNDetAnalysis::ReadMCTrack()
  192. {
  193. //Test of MC tracks and points ...
  194. FairMCTrack *tr;
  195. FairMCTrack *mtr;
  196. MpdNDetPoint *wp;
  197. MpdNDetPointLite *lp;
  198. TVector3 mom3;
  199. Int_t n=fMCtrack->GetEntries();
  200. cout<<"Number of MCTracks: "<<n<<endl;
  201. Int_t lite = fNDetPointLite->GetEntries();
  202. cout<<"Number of NDETPointLite (inside volume): "<<lite<<endl;
  203. Int_t point = fNDetPoint->GetEntries();
  204. cout<<"Number of NDETPoint (on front wall): "<<point<<endl;
  205. Double_t Mn=0.939565560;
  206. //loop over front hits(points)
  207. for(Int_t i=0;i<point;i++) {
  208. wp = (MpdNDetPoint*)fNDetPoint->At(i);
  209. if(wp->TrackID()>-1&&wp->TrackID()<n+1)
  210. tr = (FairMCTrack*)fMCtrack->At(wp->TrackID());
  211. else {cout<<"Wrong TrackID at MpdNDetPoint"<<endl; exit(-1);}
  212. Int_t MoPdg = -111111111;
  213. // if(tr->GetMotherID()>-1) {
  214. if(tr->GetMotherId()>-1) {
  215. // mtr=(FairMCTrack*)fMCtrack->At(tr->GetMotherID());
  216. mtr=(FairMCTrack*)fMCtrack->At(tr->GetMotherId());
  217. MoPdg=mtr->GetPdgCode();
  218. }
  219. if(tr->GetPdgCode() == 2212)//neutron
  220. cout<<"TID "<<wp->TrackID()
  221. <<"\t C: "<<wp->X()<<" "<<wp->Y()<<" "<<wp->Z()
  222. <<"\t P: "<<wp->Px()<<" "<<wp->Py()<<" "<<wp->Pz()
  223. <<"\t Ekin: "<<wp->ELoss()
  224. //<<" E "<<TMath::Sqrt(wp->Px()*wp->Px()+wp->Py()*wp->Py()+wp->Pz()*wp->Pz()+Mn*Mn)-Mn
  225. <<"\t T: "<<wp->Time()
  226. <<"\t L: "<<wp->Length()
  227. <<"\t Di: "<<wp->DetectorID()
  228. <<"\t Ekin: "<<wp->ELoss()
  229. <<"\t PDG: "<<tr->GetPdgCode()
  230. <<" MoPGD "<<MoPdg
  231. <<endl;
  232. }
  233. cout<<"===================================="<<endl;
  234. cout<<"===================================="<<endl;
  235. cout<<"===================================="<<endl;
  236. cout<<"===================================="<<endl;
  237. //loop over hits(points) inside ndet
  238. for(Int_t i=0;i<lite;i++) {
  239. lp = (MpdNDetPointLite*)fNDetPointLite->At(i);
  240. if(lp->GetTrackID()>-1&&lp->GetTrackID()<n+1)
  241. tr = (FairMCTrack*)fMCtrack->At(lp->GetTrackID());
  242. else {cout<<"Wrong TrackID at MpdNDetPointLite"<<endl; exit(-1);}
  243. Int_t MoPdg = -111111111;
  244. // if(tr->GetMotherID()>-1) {
  245. if(tr->GetMotherId()>-1) {
  246. // mtr=(FairMCTrack*)fMCtrack->At(tr->GetMotherID());
  247. mtr=(FairMCTrack*)fMCtrack->At(tr->GetMotherId());
  248. MoPdg=mtr->GetPdgCode();
  249. }
  250. if(tr->GetPdgCode() == 12212)//neutron
  251. cout<<"TID "<<lp->GetTrackID()
  252. <<"\t T "<<lp->GetTime()
  253. <<"\t Eloss "<<lp->GetEnergyLoss()
  254. <<"\t PDG: "<<tr->GetPdgCode()
  255. <<"\t MoPGD "<<MoPdg
  256. <<"\t Det "<<lp->GetDetectorID()
  257. <<endl;
  258. }
  259. }
  260. //------------------------------------------------------------------------------
  261. void MpdNDetAnalysis::Acceptance()
  262. {
  263. FairMCTrack *tr;//track
  264. FairMCTrack *mtr;//mother track
  265. MpdNDetPoint *wp;//MC pionts on the enter of ndet
  266. MpdNDetPointLite *lp;
  267. TLorentzVector momentum;
  268. TLorentzVector coordinate;
  269. TVector3 mvtx;
  270. TVector3 vtx;
  271. TVector3 tmom;
  272. Int_t n=fMCtrack->GetEntries();
  273. cout<<"Number of MCTracks: "<<n<<endl;
  274. Int_t lite = fNDetPointLite->GetEntries();
  275. cout<<"Number of NDETPointLite (inside volume): "<<lite<<endl;
  276. Int_t point = fNDetPoint->GetEntries();
  277. cout<<"Number of NDETPoint (on front wall): "<<point<<endl;
  278. Double_t Mn=0.939565560;
  279. //loop over front hits(points)
  280. for(Int_t i=0;i<point;i++) {
  281. wp = (MpdNDetPoint*)fNDetPoint->At(i);
  282. if(wp->TrackID()>-1&&wp->TrackID()<n+1)
  283. tr = (FairMCTrack*)fMCtrack->At(wp->TrackID());
  284. else {cout<<"Wrong TrackID at MpdNDetPoint"<<endl; exit(-1);}
  285. // vtx=tr->GetStartVertex();
  286. tr->GetStartVertex(vtx);
  287. // tmom=tr->GetMomentum();
  288. tr->GetMomentum(tmom);
  289. Int_t MoPdg = -111111111;
  290. // if(tr->GetMotherID()>-1) {
  291. if(tr->GetMotherId()>-1) {
  292. // mtr=(FairMCTrack*)fMCtrack->At(tr->GetMotherID());
  293. mtr=(FairMCTrack*)fMCtrack->At(tr->GetMotherId());
  294. MoPdg=mtr->GetPdgCode();
  295. // mvtx=mtr->GetStartVertex();
  296. mtr->GetStartVertex(mvtx);
  297. //if(tr->GetPdgCode() == 2112) cout<<TMath::Sqrt(mvtx.X()*mvtx.X()+mvtx.Y()*mvtx.Y())<<endl;
  298. }
  299. momentum.SetPxPyPzE(wp->Px(),wp->Py(),wp->Pz(),wp->ELoss());
  300. coordinate.SetXYZT(wp->X(),wp->Y(),wp->Z(),wp->Time());
  301. fhMcpPdgNdet->Fill(tr->GetPdgCode());
  302. Double_t rad=TMath::Sqrt(wp->X()*wp->X()+wp->Y()*wp->Y());
  303. Double_t mrad=TMath::Sqrt(vtx.X()*vtx.X()+vtx.Y()*vtx.Y());
  304. //TMath::Sqrt(mvtx.X()*mvtx.X()+mvtx.Y()*mvtx.Y());
  305. //Coming in particles
  306. if(
  307. momentum.E() > Mn+0.001 //1 MeV
  308. //&& momentum.Px()*coordinate.X()>0
  309. //&& momentum.Py()*coordinate.Y()>0
  310. && tmom.Px()*vtx.X()>=0
  311. && tmom.Py()*vtx.Y()>=0
  312. && rad<155.01
  313. && mrad<155. //cut secondary neutrons which creates in ndet
  314. ) {
  315. //neutrons
  316. if(tr->GetPdgCode() == 2112) {
  317. // if(tr->GetMotherID()==-1) {//primary
  318. if(tr->GetMotherId()==-1) {//primary
  319. fhNeut0Pabs ->Fill(momentum.P());
  320. fhNeut0Pt ->Fill(momentum.Pt());
  321. fhNeut0Eta ->Fill(momentum.Eta());
  322. fhNeut0Theta->Fill(momentum.Theta());
  323. fhNeut0Phi ->Fill(momentum.Phi());
  324. fhNeut0PtEta->Fill(momentum.Pt(),momentum.Eta());
  325. fhNeut0Rad->Fill(rad);
  326. fhNeut0RadZ->Fill(vtx.Z(),mrad);
  327. } else {
  328. //cout<<momentum.Px()<<" "<<tmom.Px()<<endl;
  329. fhNeut1Pabs ->Fill(momentum.P());
  330. fhNeut1Pt ->Fill(momentum.Pt());
  331. fhNeut1Eta ->Fill(momentum.Eta());
  332. fhNeut1Theta->Fill(momentum.Theta());
  333. fhNeut1Phi ->Fill(momentum.Phi());
  334. fhNeut1PtEta->Fill(momentum.Pt(),momentum.Eta());
  335. fhNeut1Rad->Fill(mrad);
  336. fhNeut1RadZ->Fill(vtx.Z(),mrad);
  337. //fhNeut1RadE->Fill(momentum.E()-Mn,mrad);
  338. fhNeut1RadE->Fill(momentum.P(),mrad);
  339. }
  340. }
  341. if(tr->GetPdgCode() == 2212) {//protons
  342. // if(tr->GetMotherID()==-1) {//primary
  343. if(tr->GetMotherId()==-1) {//primary
  344. fhProt0Pabs ->Fill(momentum.P());
  345. fhProt0Pt ->Fill(momentum.Pt());
  346. fhProt0Eta ->Fill(momentum.Eta());
  347. fhProt0Theta->Fill(momentum.Theta());
  348. fhProt0Phi ->Fill(momentum.Phi());
  349. fhProt0PtEta->Fill(momentum.Pt(),momentum.Eta());
  350. fhProt0Rad->Fill(rad);
  351. } else {
  352. fhProt1Pabs ->Fill(momentum.P());
  353. fhProt1Pt ->Fill(momentum.Pt());
  354. fhProt1Eta ->Fill(momentum.Eta());
  355. fhProt1Theta->Fill(momentum.Theta());
  356. fhProt1Phi ->Fill(momentum.Phi());
  357. fhProt1PtEta->Fill(momentum.Pt(),momentum.Eta());
  358. fhProt1Rad->Fill(mrad);
  359. }
  360. }
  361. }//coming in ...
  362. if(tr->GetPdgCode() == 12212)//neutron
  363. cout<<"TID "<<wp->TrackID()
  364. <<"\t C: "<<wp->X()<<" "<<wp->Y()<<" "<<wp->Z()
  365. <<"\t P: "<<wp->Px()<<" "<<wp->Py()<<" "<<wp->Pz()
  366. <<"\t Ekin: "<<wp->ELoss()
  367. //<<" E "<<TMath::Sqrt(wp->Px()*wp->Px()+wp->Py()*wp->Py()+wp->Pz()*wp->Pz()+Mn*Mn)-Mn
  368. <<"\t T: "<<wp->Time()
  369. <<"\t L: "<<wp->Length()
  370. <<"\t Di: "<<wp->DetectorID()
  371. <<"\t Etot: "<<wp->ELoss()
  372. <<"\t PDG: "<<tr->GetPdgCode()
  373. <<" MoPGD "<<MoPdg
  374. <<endl;
  375. }
  376. }
  377. //------------------------------------------------------------------------------
  378. void MpdNDetAnalysis::AnalyzePrimary()
  379. {
  380. // Neutary track analysis
  381. FairMCTrack *primary;
  382. fNprimary = 0;
  383. TVector3 mom3;
  384. TLorentzVector momentum;
  385. Float_t energyTotal = 0;
  386. cout<<"Number of MCTracks: "<<fMCtrack->GetEntries()<<endl;
  387. for (Int_t iNeutary=0; iNeutary<fMCtrack->GetEntries(); iNeutary++) {
  388. primary = (FairMCTrack*)fMCtrack->At(iNeutary);
  389. //if (primary->GetMotherID() != -1) continue;
  390. fNprimary++;
  391. // mom3=primary->GetMomentum();
  392. primary->GetMomentum(mom3);
  393. /*
  394. cout<<"# "<<iNeutary
  395. <<" ID: "<<primary->GetPdgCode()
  396. <<"MID: "<<primary->GetMotherID()
  397. <<" Px: "<<mom3.Px()
  398. <<endl;
  399. */
  400. /*
  401. if (strstr(fDebug,"prim"))
  402. printf("track %d: id=%d, p=(%f,%f,%f,%f) GeV/c, pt,y,theta=%f,%f,%f\n",
  403. iNeutary,
  404. primary->GetPdgCode(),
  405. momentum.Px(),
  406. momentum.Py(),
  407. momentum.Pz(),
  408. momentum.E(),
  409. momentum.Pt(),
  410. momentum.Eta(),
  411. momentum.Theta());
  412. fhNeutPabs ->Fill(momentum.P());
  413. fhNeutPt ->Fill(momentum.Pt());
  414. fhNeutEta ->Fill(momentum.Eta());
  415. fhNeutTheta->Fill(momentum.Theta());
  416. fhNeutPhi ->Fill(momentum.Phi());
  417. fhNeutPtEta->Fill(momentum.Pt(),momentum.Eta());
  418. energyTotal+=momentum.P();
  419. */
  420. }
  421. printf("Total energy of primaries: %F GeV\n",energyTotal);
  422. if (strstr(fDebug,"prim"))
  423. printf("\tNumber of primaries: %d\n",fNprimary);
  424. }
  425. //------------------------------------------------------------------------------
  426. void MpdNDetAnalysis::AnalyzeMCPointsEdep()
  427. {
  428. // Monte-Carlo points inside NDET
  429. if (fNDetPointLite == 0) return;
  430. if (strstr(fDebug,"mcp"))
  431. printf("\tNumber of NDET MC points in NDET: %d\n",fNDetPointLite->GetEntries());
  432. MpdNDetPointLite* mcPoint=NULL;
  433. TVector3 xyz;
  434. TLorentzVector momentum;
  435. Float_t energyTotal = 0;
  436. for (Int_t iPoint=0; iPoint<fNDetPointLite->GetEntries(); iPoint++) {
  437. mcPoint = (MpdNDetPointLite*)fNDetPointLite->At(iPoint); // MC point
  438. energyTotal+=mcPoint->GetEnergyLoss();
  439. }
  440. printf("Total energy in MC points: %F GeV\n",energyTotal);
  441. }
  442. //------------------------------------------------------------------------------
  443. void MpdNDetAnalysis::AnalyzeMCPointsWall()
  444. {
  445. // Monte-Carlo points on entrance to NDET
  446. if (fNDetPoint == 0) return;
  447. if (strstr(fDebug,"mcp"))
  448. printf("\tNumber of NDET MC points on NDET: %d\n",fNDetPoint->GetEntries());
  449. MpdNDetPoint* mcPoint=NULL;
  450. TVector3 xyz;
  451. TVector3 mom3;
  452. TLorentzVector momentum;
  453. for (Int_t iPoint=0; iPoint<fNDetPoint->GetEntries(); iPoint++) {
  454. mcPoint = (MpdNDetPoint*)fNDetPoint->At(iPoint); // MC point
  455. FairMCTrack* track = (FairMCTrack*)fMCtrack->At(mcPoint->GetTrackID());
  456. if (track->GetPdgCode() != 22 ) continue; // not a photon!
  457. mcPoint->Position(xyz);
  458. fhMcpXY->Fill(xyz.X(),xyz.Y());
  459. // mom3=track->GetMomentum();
  460. track->GetMomentum(mom3);
  461. if (strstr(fDebug,"mcp"))
  462. printf("MC point %d: id=%d, p=(%f,%f,%f,%f) GeV/c, pt,y,theta=%f,%f,%f\n",
  463. iPoint,
  464. track->GetPdgCode(),
  465. momentum.Px(),
  466. momentum.Py(),
  467. momentum.Pz(),
  468. momentum.Energy(),
  469. momentum.Pt(),
  470. momentum.Eta(),
  471. momentum.Theta());
  472. fhMcpPabs ->Fill(momentum.P());
  473. fhMcpPt ->Fill(momentum.Pt());
  474. fhMcpEta ->Fill(momentum.Eta());
  475. fhMcpTheta->Fill(momentum.Theta());
  476. fhMcpPhi ->Fill(momentum.Phi());
  477. fhMcpPtEta->Fill(momentum.Pt(),momentum.Eta());
  478. }
  479. }
  480. //------------------------------------------------------------------------------
  481. //... void MpdNDetAnalysis::AnalyzeHits()
  482. //... {
  483. //... // NDET hit analysis for Fast MC
  484. //...
  485. //... if (fListNDEThits == 0) return;
  486. //... if (strstr(fDebug,"hit"))
  487. //... printf("\tNumber of NDET hits: %d\n",fListNDEThits->GetEntries());
  488. //...
  489. //... MpdNDetHitFastMC* hit=NULL;
  490. //...
  491. //... for (Int_t iHit=0; iHit<fListNDEThits->GetEntries(); iHit++) {
  492. //... hit = (MpdNDetHitFastMC*)fListNDEThits->At(iHit);
  493. //... fhHitXY->Fill(hit->GetX(),hit->GetY());
  494. //... fhHitE ->Fill(hit->GetAmplitude());
  495. //... }
  496. //... }
  497. //------------------------------------------------------------------------------
  498. //... void MpdNDetAnalysis::AnalyzeRecParticles()
  499. //... {
  500. //... // NDET reconstructed particle analysis
  501. //...
  502. //... if (fListNDETrp == 0) return;
  503. //... if (strstr(fDebug,"rp"))
  504. //... printf("\tNumber of NDET reconstructed particles: %d\n",fListNDETrp->GetEntries());
  505. //...
  506. //... MpdNDetRecParticle* recpart=NULL;
  507. //...
  508. //... for (Int_t iRP=0; iRP<fListNDETrp->GetEntries(); iRP++) {
  509. //... recpart = (MpdNDetRecParticle*)fListNDETrp->At(iRP);
  510. //... TLorentzVector momentum = recpart->GetMomentum();
  511. //... fhRecPabs->Fill(momentum.P());
  512. //... fhRecPt ->Fill(momentum.Pt());
  513. //... fhRecYrap->Fill(momentum.Rapidity());
  514. //... fhRecPhi ->Fill(momentum.Phi());
  515. //... }
  516. //... }
  517. //------------------------------------------------------------------------------
  518. void MpdNDetAnalysis::Finish()
  519. {
  520. WriteOutput();
  521. }
  522. //------------------------------------------------------------------------------
  523. ClassImp(MpdNDetAnalysis)