EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FairMCStack.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FairMCStack.cxx
1 // -------------------------------------------------------------------------
2 // ----- FairMCStack source file -----
3 // ----- Created 09/10/08 by M. Al-Turany -----
4 // -------------------------------------------------------------------------
5 
6 #include "FairMCStack.h"
7 #include "TEveTrack.h"
8 #include <iostream>
9 #include "TEveTrackPropagator.h"
10 #include "TGeoTrack.h"
11 #include "FairRootManager.h"
12 #include "TClonesArray.h"
13 #include "TObjArray.h"
14 #include "TEveManager.h"
15 #include "FairEventManager.h"
16 #include "FairMCTrack.h"
17 #include "TGeant3.h"
18 #include "FairGeanePro.h"
19 #include "FairTrajFilter.h"
20 
21 using std::cout;
22 using std::endl;
23 
24 // ----- Default constructor -------------------------------------------
25 FairMCStack::FairMCStack()
26 {
27 
28 }
29 // -------------------------------------------------------------------------
30 
31 
32 // ----- Standard constructor ------------------------------------------
33 FairMCStack::FairMCStack(const char* name, Int_t iVerbose)
34  : FairTask(name, iVerbose),
35  fEveTrList( new TObjArray(16))
36 {
37 
38 }
39 // -------------------------------------------------------------------------
40 InitStatus FairMCStack::Init()
41 {
42  if(fVerbose>1) { cout<< "FairMCStack::Init()" << endl; }
43  FairRootManager* fManager = FairRootManager::Instance();
44  fTrackList = (TClonesArray*)fManager->GetObject("MCTrack");
45  if(fTrackList==0) {
46  cout << "FairMCPointDraw::Init() branch " << GetName() << " Not found! Task will be deactivated "<< endl;
47  SetActive(kFALSE);
48  }
49  if(fVerbose>2) { cout<< "FairMCStack::Init() get track list" << fTrackList<< endl; }
50  if(fVerbose>2) { cout<< "FairMCStack::Init() create propagator" << endl; }
51  fEventManager =FairEventManager::Instance();
52  if(fVerbose>2) { cout<< "FairMCStack::Init() get instance of FairEventManager " << endl; }
53  fEvent = "Current Event";
54  MinEnergyLimit=fEventManager->GetEvtMinEnergy();
55  MaxEnergyLimit=fEventManager->GetEvtMaxEnergy();
56  PEnergy=0;
57  fPro = new FairGeanePro();
58  gMC3 = (TGeant3*) gMC;
59  x1[0]=0;
60  x1[1]=0;
61  x1[2]=0;
62  p1[0]=0;
63  p1[1]=0;
64  p1[2]=0;
65 
66 
67  x2[0]=0;
68  x2[1]=0;
69  x2[2]=0;
70  p2[0]=0;
71  p2[1]=0;
72  p2[2]=0;
73 
74 
75  for(Int_t i=0; i<15; i++) { ein[i]=0; }
76  Float_t length[1];
77  length[0]=100.0;
78  //gMC3->Eufill(1, ein, length);
79 
80 
81  fTrajFilter = FairTrajFilter::Instance();
82 
83  if(IsActive()) { return kSUCCESS; }
84  else { return kERROR; }
85 
86 }
87 // -------------------------------------------------------------------------
88 void FairMCStack::Exec(Option_t* option)
89 {
90 
91  if (IsActive()) {
92 
93  if(fVerbose>1) { cout << " FairMCStack::Exec "<< endl; }
94  FairMCTrack* tr;
95  const Double_t* point;
96 
97  Reset();
98 
99  for (Int_t i=0; i<fTrackList->GetEntriesFast(); i++) {
100  if(fVerbose>2) { cout << "FairMCStack::Exec "<< i << endl; }
101  tr=(FairMCTrack*)fTrackList->At(i);
102 
103  TVector3 Ptot = tr->GetMomentum();
104  Int_t MotherId =tr->GetMotherID();
105  TVector3 Vertex =tr->GetStartVertex();
106 
107  Double_t time= tr->GetStartTime()*1e-09;
108 
109  x1[0]=Vertex.x();
110  x1[1]=Vertex.y();
111  x1[2]=Vertex.z();
112  p1[0]=Ptot.Px();
113  p1[1]=Ptot.Py();
114  p1[2]=Ptot.Pz();
115 
116 
117 
118  //TParticle *P=(TParticle *)tr->GetParticle();
119  Double_t mass=0.0;
120  Double_t ene=0.0;
121  TParticlePDG* fParticlePDG = TDatabasePDG::Instance()->GetParticle(tr->GetPdgCode());
122 
123  if (fParticlePDG) {
124  mass = fParticlePDG->Mass();
125  }
126  if ( mass >= 0 ) {
127  ene = TMath::Sqrt(mass*mass + Ptot.Mag2());
128  }
129  //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)
130 
131  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);
132 
133  PEnergy=P->Energy();
134  MinEnergyLimit=TMath::Min(PEnergy,MinEnergyLimit) ;
135  MaxEnergyLimit=TMath::Max(PEnergy,MaxEnergyLimit) ;
136  if(fVerbose>2) { cout << "MinEnergyLimit " << MinEnergyLimit << " MaxEnergyLimit " << MaxEnergyLimit << endl; }
137  if (fEventManager->IsPriOnly() && P->GetMother(0)>-1) { continue; }
138  if(fEventManager->GetCurrentPDG()!=0 && fEventManager->GetCurrentPDG()!= tr->GetPdgCode()) { continue; }
139  if(fVerbose>2) { cout << "PEnergy " << PEnergy << " Min " << fEventManager->GetMinEnergy() << " Max " << fEventManager->GetMaxEnergy() <<endl; }
140  if( (PEnergy<fEventManager->GetMinEnergy()) || (PEnergy >fEventManager->GetMaxEnergy())) { continue; }
141 
142  //here we have to propagate
143 
144  //Int_t Np=tr->GetNpoints();
145 
146  fTrList= GetTrGroup(P);
147  TEveTrack* track= new TEveTrack(P, tr->GetPdgCode(), fTrPr);
148  track->SetLineColor(fEventManager->Color(tr->GetPdgCode()));
149 
150  // Int_t GeantCode=TDatabasePDG::Instance()->ConvertPdgToGeant3(tr->GetPdgCode());
151 
152  // gMC3->Ertrak(x1,p1,x2,p2,GeantCode,"L");
153  fPro->PropagateToLength(100.0);
154  fPro->Propagate(x1, p1, x2, p2,tr->GetPdgCode());
155  TGeoTrack* tr1= fTrajFilter->GetCurrentTrk();
156 
157  Int_t Np=tr1->GetNpoints();
158 
159  for (Int_t n=0; n<Np; n++) {
160  point=tr1->GetPoint(n);
161  track->SetPoint(n,point[0],point[1],point[2]);
162  TEveVector pos= TEveVector(point[0], point[1],point[2]);
163  TEvePathMark* path = new TEvePathMark();
164  path->fV=pos ;
165  path->fTime= point[3];
166  if(n==0) {
167  TEveVector Mom= TEveVector(P->Px(), P->Py(),P->Pz());
168  path->fP=Mom;
169  }
170  if(fVerbose>3) { cout << "Path marker added " << path << endl; }
171 
172  track->AddPathMark(*path);
173 
174  if(fVerbose>3) { cout << "Path marker added " << path << endl; }
175  }
176  fTrList->AddElement(track);
177  if(fVerbose>3) { cout << "track added " << track->GetName() << endl; }
178 
179 
180 // cout << "CurrentTrack " << fTrajFilter->GetCurrentTrk() << endl;
181 
182  /* for (Int_t n=0; n<Np; n++){
183  point=tr->GetPoint(n);
184  track->SetPoint(n,point[0],point[1],point[2]);
185  TEveVector pos= TEveVector(point[0], point[1],point[2]);
186  TEvePathMark *path = new TEvePathMark();
187  path->fV=pos ;
188  path->fTime= point[3];
189  if(n==0){
190  TEveVector Mom= TEveVector(P->Px(), P->Py(),P->Pz());
191  path->fP=Mom;
192  }
193  if(fVerbose>3) cout << "Path marker added " << path << endl;
194  track->AddPathMark(*path);
195  if(fVerbose>3) cout << "Path marker added " << path << endl;
196  }
197 
198  fTrList->AddElement(track);
199  if(fVerbose>3)cout << "track added " << track->GetName() << endl;
200  */
201  }
202  for (Int_t i=0; i<fEveTrList->GetEntriesFast(); i++) {
203  TEveTrackList* TrListIn=( TEveTrackList*) fEveTrList->At(i);
204  TrListIn->FindMomentumLimits(TrListIn, kFALSE);
205  }
206  fEventManager->SetEvtMaxEnergy(MaxEnergyLimit);
207  fEventManager->SetEvtMinEnergy(MinEnergyLimit);
208  gEve->Redraw3D(kFALSE);
209  }
210 }
211 // ----- Destructor ----------------------------------------------------
212 FairMCStack::~FairMCStack()
213 {
214 }
215 // -------------------------------------------------------------------------
216 void FairMCStack::SetParContainers()
217 {
218 
219 }
220 
221 // -------------------------------------------------------------------------
222 void FairMCStack::Finish()
223 {
224 
225 }
226 // -------------------------------------------------------------------------
227 void FairMCStack::Reset()
228 {
229  for (Int_t i=0; i<fEveTrList->GetEntriesFast(); i++) {
230  TEveTrackList* ele=( TEveTrackList*) fEveTrList->At(i);
231  gEve->RemoveElement(ele,fEventManager);
232  }
233  fEveTrList->Clear();
234 }
235 
236 TEveTrackList* FairMCStack::GetTrGroup(TParticle* P)
237 {
238 
239  fTrList=0;
240  for (Int_t i=0; i<fEveTrList->GetEntriesFast(); i++) {
241  TEveTrackList* TrListIn=( TEveTrackList*) fEveTrList->At(i);
242  if ( strcmp(TrListIn->GetName(),P->GetName())==0 ) {
243  fTrList= TrListIn;
244  break;
245  }
246  }
247  if(fTrList ==0) {
248  fTrPr=new TEveTrackPropagator();
249  fTrList= new TEveTrackList(P->GetName(),fTrPr );
250  fTrList->SetMainColor(fEventManager->Color(P->GetPdgCode()));
251  fEveTrList->Add(fTrList);
252  gEve->AddElement( fTrList ,fEventManager );
253  fTrList->SetRnrLine(kTRUE);
254  }
255  return fTrList;
256 }
257 
258 ClassImp(FairMCStack)
259 
260