FairRoot
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MyProjStack.cxx
Go to the documentation of this file.
1 /********************************************************************************
2  * Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH *
3  * *
4  * This software is distributed under the terms of the *
5  * GNU Lesser General Public Licence (LGPL) version 3, *
6  * copied verbatim in the file "LICENSE" *
7  ********************************************************************************/
8 
9 // -------------------------------------------------------------------------
10 // ----- MyProjStack source file -----
11 // ----- M. Al-Turany June 2014 -----
12 // -------------------------------------------------------------------------
13 
14 #include "MyProjStack.h"
15 
16 #include "FairDetector.h" // for FairDetector
17 #include "FairLink.h" // for FairLink
18 #include "FairLogger.h" // for FairLogger, MESSAGE_ORIGIN
19 #include "FairMCPoint.h" // for FairMCPoint
20 #include "FairRootManager.h" // for FairRootManager
21 #include "MyProjMCTrack.h" // for MyProjMCTrack
22 
23 #include <TClonesArray.h> // for TClonesArray
24 #include <TIterator.h> // for TIterator
25 #include <TLorentzVector.h> // for TLorentzVector
26 #include <TParticle.h> // for TParticle
27 #include <TRefArray.h> // for TRefArray
28 #include <TVirtualMC.h> // for gMC
29 #include <iosfwd> // for ostream
30 #include <iostream> // for operator<<, etc
31 #include <stddef.h> // for NULL
32 
33 using std::cout;
34 using std::endl;
35 using std::pair;
36 
37 // ----- Default constructor -------------------------------------------
38 MyProjStack::MyProjStack(Int_t size)
40  , fStack()
41  , fParticles(new TClonesArray("TParticle", size))
42  , fTracks(new TClonesArray("MyProjMCTrack", size))
43  , fStoreMap()
44  , fStoreIter()
45  , fIndexMap()
46  , fIndexIter()
47  , fPointsMap()
48  , fCurrentTrack(-1)
49  , fNPrimaries(0)
50  , fNParticles(0)
51  , fNTracks(0)
52  , fIndex(0)
53  , fStoreSecondaries(kTRUE)
54  , fMinPoints(0)
55  , fEnergyCut(0.)
56  , fStoreMothers(kTRUE)
57 {}
58 
59 // -------------------------------------------------------------------------
60 
62  : FairGenericStack(right)
63  , fStack()
64  , fParticles(new TClonesArray("TParticle", 100))
65  , fTracks(new TClonesArray("MyProjMCTrack", 100))
66  , fStoreMap()
67  , fStoreIter()
68  , fIndexMap()
69  , fIndexIter()
70  , fPointsMap()
71  , fCurrentTrack()
72  , fNPrimaries()
73  , fNParticles()
74  , fNTracks()
75  , fIndex()
76  , fStoreSecondaries(right.fStoreSecondaries)
77  , fMinPoints(right.fMinPoints)
78  , fEnergyCut(right.fEnergyCut)
79  , fStoreMothers(right.fStoreMothers)
80 {}
81 
82 // ----- Destructor ----------------------------------------------------
84 {
85  if (fParticles) {
86  fParticles->Delete();
87  delete fParticles;
88  }
89  if (fTracks) {
90  fTracks->Delete();
91  delete fTracks;
92  }
93 }
94 // -------------------------------------------------------------------------
95 
96 void MyProjStack::PushTrack(Int_t toBeDone,
97  Int_t parentId,
98  Int_t pdgCode,
99  Double_t px,
100  Double_t py,
101  Double_t pz,
102  Double_t e,
103  Double_t vx,
104  Double_t vy,
105  Double_t vz,
106  Double_t time,
107  Double_t polx,
108  Double_t poly,
109  Double_t polz,
110  TMCProcess proc,
111  Int_t& ntr,
112  Double_t weight,
113  Int_t is)
114 {
115 
116  PushTrack(
117  toBeDone, parentId, pdgCode, px, py, pz, e, vx, vy, vz, time, polx, poly, polz, proc, ntr, weight, is, -1);
118 }
119 
120 // ----- Virtual public method PushTrack -------------------------------
121 void MyProjStack::PushTrack(Int_t toBeDone,
122  Int_t parentId,
123  Int_t pdgCode,
124  Double_t px,
125  Double_t py,
126  Double_t pz,
127  Double_t e,
128  Double_t vx,
129  Double_t vy,
130  Double_t vz,
131  Double_t time,
132  Double_t polx,
133  Double_t poly,
134  Double_t polz,
135  TMCProcess proc,
136  Int_t& ntr,
137  Double_t weight,
138  Int_t is,
139  Int_t secondparentID)
140 {
141 
142  // --> Get TParticle array
143  TClonesArray& partArray = *fParticles;
144 
145  // --> Create new TParticle and add it to the TParticle array
146  Int_t trackId = fNParticles;
147  Int_t nPoints = 0;
148  Int_t daughter1Id = -1;
149  Int_t daughter2Id = -1;
150  TParticle* particle = new (partArray[fNParticles++])
151  TParticle(pdgCode, trackId, parentId, nPoints, daughter1Id, daughter2Id, px, py, pz, e, vx, vy, vz, time);
152  particle->SetPolarisation(polx, poly, polz);
153  particle->SetWeight(weight);
154  particle->SetUniqueID(proc);
155 
156  // --> Increment counter
157  if (parentId < 0) {
158  fNPrimaries++;
159  }
160 
161  // --> Set argument variable
162  ntr = trackId;
163 
164  // --> Push particle on the stack if toBeDone is set
165  if (toBeDone == 1) {
166  fStack.push(particle);
167  }
168 }
169 // -------------------------------------------------------------------------
170 
171 // ----- Virtual method PopNextTrack -----------------------------------
172 TParticle* MyProjStack::PopNextTrack(Int_t& iTrack)
173 {
174 
175  // If end of stack: Return empty pointer
176  if (fStack.empty()) {
177  iTrack = -1;
178  return NULL;
179  }
180 
181  // If not, get next particle from stack
182  TParticle* thisParticle = fStack.top();
183  fStack.pop();
184 
185  if (!thisParticle) {
186  iTrack = 0;
187  return NULL;
188  }
189 
190  fCurrentTrack = thisParticle->GetStatusCode();
191  iTrack = fCurrentTrack;
192 
193  return thisParticle;
194 }
195 // -------------------------------------------------------------------------
196 
197 // ----- Virtual method PopPrimaryForTracking --------------------------
198 TParticle* MyProjStack::PopPrimaryForTracking(Int_t iPrim)
199 {
200 
201  // Get the iPrimth particle from the fStack TClonesArray. This
202  // should be a primary (if the index is correct).
203 
204  // Test for index
205  if (iPrim < 0 || iPrim >= fNPrimaries) {
206  LOG(fatal) << "MyProjStack: Primary index out of range! " << iPrim;
207  }
208 
209  // Return the iPrim-th TParticle from the fParticle array. This should be
210  // a primary.
211  TParticle* part = (TParticle*)fParticles->At(iPrim);
212  if (!(part->GetMother(0) < 0)) {
213  LOG(fatal) << "MyProjStack:: Not a primary track! " << iPrim;
214  }
215 
216  return part;
217 }
218 // -------------------------------------------------------------------------
219 
220 // ----- Virtual public method GetCurrentTrack -------------------------
221 TParticle* MyProjStack::GetCurrentTrack() const
222 {
223  TParticle* currentPart = GetParticle(fCurrentTrack);
224  if (!currentPart) {
225  LOG(warning) << "MyProjStack: Current track not found in stack!";
226  }
227  return currentPart;
228 }
229 // -------------------------------------------------------------------------
230 
231 // ----- Public method AddParticle -------------------------------------
232 void MyProjStack::AddParticle(TParticle* oldPart)
233 {
234  TClonesArray& array = *fParticles;
235  TParticle* newPart = new (array[fIndex]) TParticle(*oldPart);
236  newPart->SetWeight(oldPart->GetWeight());
237  newPart->SetUniqueID(oldPart->GetUniqueID());
238  fIndex++;
239 }
240 // -------------------------------------------------------------------------
241 
242 // ----- Public method FillTrackArray ----------------------------------
244 {
245 
246  LOG(debug) << "MyProjStack: Filling MCTrack array...";
247 
248  // --> Reset index map and number of output tracks
249  fIndexMap.clear();
250  fNTracks = 0;
251 
252  // --> Check tracks for selection criteria
253  // SelectTracks();
254 
255  // --> Loop over fParticles array and copy selected tracks
256  for (Int_t iPart = 0; iPart < fNParticles; iPart++) {
257 
258  /*fStoreIter = fStoreMap.find(iPart);
259  if (fStoreIter == fStoreMap.end() ) {
260  LOG(fatal) << "MyProjStack: Particle " << iPart << " not found in storage map! ";
261  }
262  Bool_t store = (*fStoreIter).second;*/
263 
264  Bool_t store = kTRUE;
265 
266  if (store) {
267  MyProjMCTrack* track = new ((*fTracks)[fNTracks]) MyProjMCTrack(GetParticle(iPart));
268  fIndexMap[iPart] = fNTracks;
269  // --> Set the number of points in the detectors for this track
270  for (Int_t iDet = kNewDetector; iDet < kSTOPHERE; iDet++) {
271  pair<Int_t, Int_t> a(iPart, iDet);
272  track->SetNPoints(iDet, fPointsMap[a]);
273  }
274  fNTracks++;
275  } else {
276  fIndexMap[iPart] = -2;
277  }
278  }
279 
280  // --> Map index for primary mothers
281  fIndexMap[-1] = -1;
282 
283  // --> Screen output
284  Print(1);
285 }
286 // -------------------------------------------------------------------------
287 
288 // ----- Public method UpdateTrackIndex --------------------------------
289 void MyProjStack::UpdateTrackIndex(TRefArray* detList)
290 {
291  /*
292  LOG(debug) << "MyProjStack: Updating track indizes...";
293  Int_t nColl = 0;
294 
295  // First update mother ID in MCTracks
296  for (Int_t i=0; i<fNTracks; i++) {
297  MyProjMCTrack* track = (MyProjMCTrack*)fTracks->At(i);
298  Int_t iMotherOld = track->GetMotherId();
299  fIndexIter = fIndexMap.find(iMotherOld);
300  if (fIndexIter == fIndexMap.end()) {
301  LOG(fatal) << "MyProjStack: Particle index " << iMotherOld << " not found in dex map! ";
302  }
303  track->SetMotherId( (*fIndexIter).second );
304  }
305 
306 
307  if(fDetList==0) {
308  // Now iterate through all active detectors
309  fDetIter = detList->MakeIterator();
310  fDetIter->Reset();
311  } else {
312  fDetIter->Reset();
313  }
314 
315  FairDetector* det = NULL;
316  while( (det = (FairDetector*)fDetIter->Next() ) ) {
317 
318 
319  // --> Get hit collections from detector
320  Int_t iColl = 0;
321  TClonesArray* hitArray;
322  while ( (hitArray = det->GetCollection(iColl++)) ) {
323  nColl++;
324  Int_t nPoints = hitArray->GetEntriesFast();
325 
326  // --> Update track index for all MCPoints in the collection
327  for (Int_t iPoint=0; iPoint<nPoints; iPoint++) {
328  FairMCPoint* point = (FairMCPoint*)hitArray->At(iPoint);
329  Int_t iTrack = point->GetTrackID();
330 
331  fIndexIter = fIndexMap.find(iTrack);
332  if (fIndexIter == fIndexMap.end()) {
333  LOG(fatal) << "MyProjStack: Particle index " << iTrack << " not found in index map! ";
334  }
335  point->SetTrackID((*fIndexIter).second);
336  point->SetLink(FairLink("MCTrack", (*fIndexIter).second));
337  }
338 
339  } // Collections of this detector
340  } // List of active detectors
341  LOG(debug) << "...stack and " << nColl << " collections updated.";
342  */
343 }
344 // -------------------------------------------------------------------------
345 
346 // ----- Public method Reset -------------------------------------------
347 void MyProjStack::Reset()
348 {
349  fIndex = 0;
350  fCurrentTrack = -1;
351  fNPrimaries = fNParticles = fNTracks = 0;
352  while (!fStack.empty()) {
353  fStack.pop();
354  }
355  fParticles->Clear();
356  fTracks->Clear();
357  fPointsMap.clear();
358 }
359 // -------------------------------------------------------------------------
360 
361 // ----- Public method Register ----------------------------------------
363 {
364  if (gMC && (!gMC->IsMT())) {
365  FairRootManager::Instance()->Register("MCTrack", "Stack", fTracks, kTRUE);
366  } else {
367  FairRootManager::Instance()->RegisterAny("MCTrack", fTracks, kTRUE);
368  }
369 }
370 // -------------------------------------------------------------------------
371 
372 // ----- Public method Print --------------------------------------------
373 void MyProjStack::Print(Int_t iVerbose) const
374 {
375  cout << "-I- MyProjStack: Number of primaries = " << fNPrimaries << endl;
376  cout << " Total number of particles = " << fNParticles << endl;
377  cout << " Number of tracks in output = " << fNTracks << endl;
378  if (iVerbose) {
379  for (Int_t iTrack = 0; iTrack < fNTracks; iTrack++) {
380  ((MyProjMCTrack*)fTracks->At(iTrack))->Print(iTrack);
381  }
382  }
383 }
384 // -------------------------------------------------------------------------
385 
386 // ----- Public method AddPoint (for current track) --------------------
388 {
389  Int_t iDet = detId;
390  // cout << "Add point for Detektor" << iDet << endl;
391  pair<Int_t, Int_t> a(fCurrentTrack, iDet);
392  if (fPointsMap.find(a) == fPointsMap.end()) {
393  fPointsMap[a] = 1;
394  } else {
395  fPointsMap[a]++;
396  }
397 }
398 // -------------------------------------------------------------------------
399 
400 // ----- Public method AddPoint (for arbitrary track) -------------------
401 void MyProjStack::AddPoint(DetectorId detId, Int_t iTrack)
402 {
403  if (iTrack < 0) {
404  return;
405  }
406  Int_t iDet = detId;
407  pair<Int_t, Int_t> a(iTrack, iDet);
408  if (fPointsMap.find(a) == fPointsMap.end()) {
409  fPointsMap[a] = 1;
410  } else {
411  fPointsMap[a]++;
412  }
413 }
414 // -------------------------------------------------------------------------
415 
416 // ----- Virtual method GetCurrentParentTrackNumber --------------------
418 {
419  TParticle* currentPart = GetCurrentTrack();
420  if (currentPart) {
421  return currentPart->GetFirstMother();
422  } else {
423  return -1;
424  }
425 }
426 // -------------------------------------------------------------------------
427 
428 // ----- Public method GetParticle -------------------------------------
429 TParticle* MyProjStack::GetParticle(Int_t trackID) const
430 {
431  if (trackID < 0 || trackID >= fNParticles) {
432  LOG(fatal) << "MyProjStack: Particle index " << trackID << " out of range.";
433  }
434  return (TParticle*)fParticles->At(trackID);
435 }
436 // -------------------------------------------------------------------------
437 
438 // ----- Private method SelectTracks -----------------------------------
439 void MyProjStack::SelectTracks()
440 {
441 
442  // --> Clear storage map
443  fStoreMap.clear();
444 
445  // --> Check particles in the fParticle array
446  for (Int_t i = 0; i < fNParticles; i++) {
447 
448  TParticle* thisPart = GetParticle(i);
449  Bool_t store = kTRUE;
450 
451  // --> Get track parameters
452  Int_t iMother = thisPart->GetMother(0);
453  TLorentzVector p;
454  thisPart->Momentum(p);
455  Double_t energy = p.E();
456  Double_t mass = p.M();
457  // Double_t mass = thisPart->GetMass();
458  Double_t eKin = energy - mass;
459 
460  // --> Calculate number of points
461  Int_t nPoints = 0;
462  for (Int_t iDet = kNewDetector; iDet < kSTOPHERE; iDet++) {
463  pair<Int_t, Int_t> a(i, iDet);
464  if (fPointsMap.find(a) != fPointsMap.end()) {
465  nPoints += fPointsMap[a];
466  }
467  }
468 
469  // --> Check for cuts (store primaries in any case)
470  if (iMother < 0) {
471  store = kTRUE;
472  } else {
473  if (!fStoreSecondaries) {
474  store = kFALSE;
475  }
476  if (nPoints < fMinPoints) {
477  store = kFALSE;
478  }
479  if (eKin < fEnergyCut) {
480  store = kFALSE;
481  }
482  }
483 
484  // --> Set storage flag
485  fStoreMap[i] = store;
486  }
487 
488  // --> If flag is set, flag recursively mothers of selected tracks
489  if (fStoreMothers) {
490  for (Int_t i = 0; i < fNParticles; i++) {
491  if (fStoreMap[i]) {
492  Int_t iMother = GetParticle(i)->GetMother(0);
493  while (iMother >= 0) {
494  fStoreMap[iMother] = kTRUE;
495  iMother = GetParticle(iMother)->GetMother(0);
496  }
497  }
498  }
499  }
500 }
501 // -------------------------------------------------------------------------
502 
DetectorId
virtual void FillTrackArray()
virtual TParticle * PopPrimaryForTracking(Int_t iPrim)
static FairRootManager * Instance()
virtual TParticle * PopNextTrack(Int_t &iTrack)
ClassImp(FairEventBuilder)
TParticle * GetParticle(Int_t trackId) const
virtual void Register()
virtual TParticle * GetCurrentTrack() const
virtual void Print(Int_t iVerbose=0) const
virtual void AddParticle(TParticle *part)
virtual void PushTrack(Int_t toBeDone, Int_t parentID, Int_t pdgCode, Double_t px, Double_t py, Double_t pz, Double_t e, Double_t vx, Double_t vy, Double_t vz, Double_t time, Double_t polx, Double_t poly, Double_t polz, TMCProcess proc, Int_t &ntr, Double_t weight, Int_t is)
Definition: MyProjStack.cxx:96
void SetNPoints(Int_t iDet, Int_t np)
virtual ~MyProjStack()
Definition: MyProjStack.cxx:83
void Register(const char *name, const char *Foldername, TNamed *obj, Bool_t toFile)
virtual void UpdateTrackIndex(TRefArray *detArray=0)
virtual Int_t GetCurrentParentTrackNumber() const
void AddPoint(DetectorId iDet)
virtual void Reset()
void RegisterAny(const char *name, T *&obj, Bool_t toFile)
MyProjStack(Int_t size=100)
Definition: MyProjStack.cxx:38