MpdMCStack.cxx 7.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237
  1. // -------------------------------------------------------------------------
  2. // ----- MpdMCStack source file -----
  3. // ----- Created 09/10/08 by M. Al-Turany -----
  4. // -------------------------------------------------------------------------
  5. #include "MpdMCStack.h"
  6. #ifdef BMNROOT
  7. #include "CbmMCTrack.h"
  8. #else
  9. #include "MpdMCTrack.h"
  10. #endif
  11. #include "FairRootManager.h"
  12. #include "FairLogger.h"
  13. #include "TGeoTrack.h"
  14. #include "TClonesArray.h"
  15. #include "TEveManager.h"
  16. #include "TParticlePDG.h"
  17. #include "TDatabasePDG.h"
  18. #include <iostream>
  19. using namespace std;
  20. // ----- Default constructor -------------------------------------------
  21. MpdMCStack::MpdMCStack()
  22. {}
  23. // ----- Standard constructor ------------------------------------------
  24. MpdMCStack::MpdMCStack(const char* name, Int_t iVerbose)
  25. : FairTask(name, iVerbose),
  26. fEveTrList(new TObjArray(16))
  27. {}
  28. // ----- Destructor ----------------------------------------------------
  29. MpdMCStack::~MpdMCStack()
  30. {}
  31. // -------------------------------------------------------------------------
  32. InitStatus MpdMCStack::Init()
  33. {
  34. if (fVerbose > 1) cout<<"MpdMCStack::Init()"<<endl;
  35. FairRootManager* fManager = FairRootManager::Instance();
  36. fTrackList = (TClonesArray*) fManager->GetObject("MCTrack");
  37. if (fTrackList == 0)
  38. {
  39. LOG(ERROR)<<"MpdMCStack::Init() branch "<<GetName()<<" not found! Task will be deactivated";
  40. SetActive(kFALSE);
  41. }
  42. if (fVerbose > 2)
  43. cout<<"MpdMCStack::Init() get track list"<<fTrackList<<endl;
  44. fEventManager = MpdEventManager::Instance();
  45. if (fVerbose > 2) cout<<"MpdMCStack::Init() get instance of MpdEventManager"<<endl;
  46. MinEnergyLimit = fEventManager->GetEvtMinEnergy();
  47. MaxEnergyLimit = fEventManager->GetEvtMaxEnergy();
  48. PEnergy = 0;
  49. fPro = new FairGeanePro();
  50. gMC3 = (TGeant3*) gMC;
  51. x1[0] = 0; x1[1] = 0; x1[2] = 0;
  52. p1[0] = 0; p1[1] = 0; p1[2] = 0;
  53. x2[0] = 0; x2[1] = 0; x2[2] = 0;
  54. p2[0] = 0; p2[1] = 0; p2[2] = 0;
  55. for (Int_t i = 0; i < 15; i++)
  56. ein[i] = 0;
  57. //Float_t length[1];
  58. //length[0] = 100.0;
  59. //gMC3->Eufill(1, ein, length);
  60. fTrajFilter = FairTrajFilter::Instance();
  61. if (IsActive())
  62. return kSUCCESS;
  63. else
  64. return kERROR;
  65. }
  66. // -------------------------------------------------------------------------
  67. void MpdMCStack::Exec(Option_t* /*option*/)
  68. {
  69. if (!IsActive())
  70. return;
  71. if (fVerbose > 1) cout<<"MpdMCStack::Exec"<<endl;
  72. Reset();
  73. for (Int_t i = 0; i < fTrackList->GetEntriesFast(); i++)
  74. {
  75. if (fVerbose > 2) cout<<"MpdMCStack::Exec "<<i<<endl;
  76. #ifdef BMNROOT
  77. CbmMCTrack* tr = (CbmMCTrack*) fTrackList->At(i);
  78. #else
  79. MpdMCTrack* tr = (MpdMCTrack*) fTrackList->At(i);
  80. #endif
  81. TVector3 Ptot;
  82. tr->GetMomentum(Ptot);
  83. Int_t MotherId =tr->GetMotherId();
  84. TVector3 Vertex;
  85. tr->GetStartVertex(Vertex);
  86. Double_t time = tr->GetStartT()*1e-09;
  87. x1[0] = Vertex.x(); x1[1] = Vertex.y(); x1[2] = Vertex.z();
  88. p1[0] = Ptot.Px(); p1[1] = Ptot.Py(); p1[2] = Ptot.Pz();
  89. //TParticle* P = (TParticle*) tr->GetParticle();
  90. TParticlePDG* fParticlePDG = TDatabasePDG::Instance()->GetParticle(tr->GetPdgCode());
  91. Double_t mass = 0, ene = 0;
  92. if (fParticlePDG)
  93. mass = fParticlePDG->Mass();
  94. if (mass >= 0)
  95. ene = TMath::Sqrt(mass*mass + Ptot.Mag2());
  96. //TParticle(Int_t pdg, Int_t status, Int_t mother1, Int_t mother2, Int_t daughter1, Int_t daughter2, Double_t px, Double_t py, Double_t pz, Double_t etot, Double_t vx, Double_t vy, Double_t vz, Double_t time)
  97. TParticle* P = new TParticle(tr->GetPdgCode(), i, MotherId, -1, -1, -1, Ptot.Px(), Ptot.Py(),Ptot.Pz(),ene, Vertex.x(), Vertex.z(), Vertex.z(), time);
  98. PEnergy = P->Energy();
  99. MinEnergyLimit = TMath::Min(PEnergy, MinEnergyLimit);
  100. MaxEnergyLimit = TMath::Max(PEnergy, MaxEnergyLimit);
  101. if (fVerbose > 2)
  102. cout<<"MinEnergyLimit "<<MinEnergyLimit<<" MaxEnergyLimit "<<MaxEnergyLimit<<endl;
  103. if ((fEventManager->IsPriOnly() && (P->GetMother(0) > -1)))
  104. continue;
  105. if ((fEventManager->GetCurrentPDG() != 0) && (fEventManager->GetCurrentPDG() != tr->GetPdgCode()))
  106. continue;
  107. if (fVerbose > 2)
  108. cout<<"PEnergy "<<PEnergy<<" Min "<<fEventManager->GetMinEnergy()<<" Max "<<fEventManager->GetMaxEnergy()<<endl;
  109. if ((PEnergy < fEventManager->GetMinEnergy()) || (PEnergy > fEventManager->GetMaxEnergy()))
  110. continue;
  111. //Int_t Np = tr->GetNpoints();
  112. fTrList = GetTrGroup(P);
  113. TEveTrack* track= new TEveTrack(P, tr->GetPdgCode(), fTrPr);
  114. track->SetLineColor(fEventManager->Color(tr->GetPdgCode()));
  115. // PROPAGATION
  116. // Int_t GeantCode = TDatabasePDG::Instance()->ConvertPdgToGeant3(tr->GetPdgCode());
  117. // gMC3->Ertrak(x1, p1, x2, p2,G eantCode, "L");
  118. fPro->PropagateToLength(100.0);
  119. // fPro->SetDestinationLength(100.0);
  120. fPro->Propagate(x1, p1, x2, p2, tr->GetPdgCode());
  121. TGeoTrack* tr1 = fTrajFilter->GetCurrentTrk();
  122. const Double_t* point;
  123. Int_t Np = tr1->GetNpoints();
  124. for (Int_t n = 0; n < Np; n++)
  125. {
  126. point = tr1->GetPoint(n);
  127. track->SetPoint(n, point[0], point[1], point[2]);
  128. TEveVector pos = TEveVector(point[0], point[1], point[2]);
  129. TEvePathMark* path = new TEvePathMark();
  130. path->fV = pos ;
  131. path->fTime = point[3];
  132. if (n == 0)
  133. {
  134. TEveVector Mom = TEveVector(P->Px(), P->Py(),P->Pz());
  135. path->fP = Mom;
  136. }
  137. track->AddPathMark(*path);
  138. if (fVerbose > 3) cout<<"Path marker added "<<path<<endl;
  139. }
  140. fTrList->AddElement(track);
  141. if (fVerbose > 3) cout<<"Track added "<<track->GetName()<<endl;
  142. }
  143. //for (Int_t i = 0; i < fEveTrList->GetEntriesFast(); i++)
  144. //{
  145. // TEveTrackList* TrListIn = (TEveTrackList*) fEveTrList->At(i);
  146. // TrListIn->FindMomentumLimits(TrListIn, kFALSE);
  147. //}
  148. fEventManager->SetEvtMinEnergy(MinEnergyLimit);
  149. fEventManager->SetEvtMaxEnergy(MaxEnergyLimit);
  150. //gEve->Redraw3D(kFALSE);
  151. }
  152. // -------------------------------------------------------------------------
  153. void MpdMCStack::SetParContainers()
  154. {}
  155. // -------------------------------------------------------------------------
  156. void MpdMCStack::Finish()
  157. {}
  158. // -------------------------------------------------------------------------
  159. void MpdMCStack::Reset()
  160. {
  161. for (Int_t i = 0; i < fEveTrList->GetEntriesFast(); i++)
  162. {
  163. TEveTrackList* ele = (TEveTrackList*) fEveTrList->At(i);
  164. gEve->RemoveElement(ele, fEventManager);
  165. }
  166. fEveTrList->Clear();
  167. }
  168. TEveTrackList* MpdMCStack::GetTrGroup(TParticle* P)
  169. {
  170. fTrList = 0;
  171. for (Int_t i = 0; i < fEveTrList->GetEntriesFast(); i++)
  172. {
  173. TEveTrackList* TrListIn = (TEveTrackList*) fEveTrList->At(i);
  174. if (strcmp(TrListIn->GetName(), P->GetName()) == 0)
  175. {
  176. fTrList= TrListIn;
  177. break;
  178. }
  179. }
  180. if (fTrList == 0)
  181. {
  182. fTrPr = new TEveTrackPropagator();
  183. fTrList = new TEveTrackList(P->GetName(), fTrPr);
  184. fTrList->SetMainColor(fEventManager->Color(P->GetPdgCode()));
  185. fTrList->SetRnrLine(kTRUE);
  186. fEveTrList->Add(fTrList);
  187. gEve->AddElement(fTrList, fEventManager);
  188. }
  189. return fTrList;
  190. }
  191. ClassImp(MpdMCStack)