forked from mkrzewic/AliceO2
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathStack.cxx
More file actions
480 lines (416 loc) · 13 KB
/
Stack.cxx
File metadata and controls
480 lines (416 loc) · 13 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
/// \file Stack.cxx
/// \brief Implementation of the Stack class
/// \author M. Al-Turany - June 2014
#include "Stack.h"
#include <stddef.h> // for NULL
#include "FairDetector.h" // for FairDetector
#include "FairLogger.h" // for MESSAGE_ORIGIN, FairLogger
#include "FairMCPoint.h" // for FairMCPoint
#include "FairRootManager.h" // for FairRootManager
#include "MCTrack.h" // for MCTrack
#include "TClonesArray.h" // for TClonesArray
#include "TIterator.h" // for TIterator
#include "TLorentzVector.h" // for TLorentzVector
#include "TParticle.h" // for TParticle
#include "TRefArray.h" // for TRefArray
using std::cout;
using std::endl;
using std::pair;
using namespace AliceO2::Data;
Stack::Stack(Int_t size)
: FairGenericStack(),
mStack(),
mParticles(new TClonesArray("TParticle", size)),
mTracks(new TClonesArray("MCTrack", size)),
mStoreMap(),
mStoreIterator(),
mIndexMap(),
mIndexIterator(),
mPointsMap(),
mIndexOfCurrentTrack(-1),
mNumberOfPrimaryParticles(0),
mNumberOfEntriesInParticles(0),
mNumberOfEntriesInTracks(0),
mIndex(0),
mStoreMothers(kTRUE),
mStoreSecondaries(kTRUE),
mMinPoints(1),
mEnergyCut(0.),
mLogger(FairLogger::GetLogger())
{
}
Stack::Stack(const Stack& rhs)
: FairGenericStack(rhs),
mStack(),
mParticles(0),
mTracks(0),
mStoreMap(),
mStoreIterator(),
mIndexMap(),
mIndexIterator(),
mPointsMap(),
mIndexOfCurrentTrack(-1),
mNumberOfPrimaryParticles(0),
mNumberOfEntriesInParticles(0),
mNumberOfEntriesInTracks(0),
mIndex(0),
mStoreMothers(rhs.mStoreMothers),
mStoreSecondaries(rhs.mStoreSecondaries),
mMinPoints(rhs.mMinPoints),
mEnergyCut(rhs.mEnergyCut),
mLogger(0)
{
mParticles = new TClonesArray("TParticle", rhs.mParticles->GetSize());
mTracks = new TClonesArray("MCTrack", rhs.mTracks->GetSize());
}
Stack::~Stack()
{
if (mParticles) {
mParticles->Delete();
delete mParticles;
}
if (mTracks) {
mTracks->Delete();
delete mTracks;
}
}
Stack& Stack::operator=(const Stack& rhs)
{
// check assignment to self
if (this == &rhs) return *this;
// base class assignment
FairGenericStack::operator=(rhs);
// assignment operator
mParticles = new TClonesArray("TParticle", rhs.mParticles->GetSize());
mTracks = new TClonesArray("MCTrack", rhs.mTracks->GetSize());
mIndexOfCurrentTrack = -1;
mNumberOfPrimaryParticles = 0;
mNumberOfEntriesInParticles = 0;
mNumberOfEntriesInTracks = 0;
mIndex = 0;
mStoreMothers = rhs.mStoreMothers;
mStoreSecondaries = rhs.mStoreSecondaries;
mMinPoints = rhs.mMinPoints;
mEnergyCut = rhs.mEnergyCut;
mLogger = 0;
return *this;
}
void Stack::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)
{
PushTrack(toBeDone, parentId, pdgCode, px, py, pz, e, vx, vy, vz, time, polx, poly, polz, proc, ntr, weight, is, -1);
}
void Stack::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, Int_t secondparentID)
{
// Get TParticle array
TClonesArray& partArray = *mParticles;
// Create new TParticle and add it to the TParticle array
Int_t trackId = mNumberOfEntriesInParticles;
Int_t nPoints = 0;
Int_t daughter1Id = -1;
Int_t daughter2Id = -1;
TParticle* particle = new (partArray[mNumberOfEntriesInParticles++])
TParticle(pdgCode, trackId, parentId, nPoints, daughter1Id, daughter2Id, px, py, pz, e, vx, vy, vz, time);
particle->SetPolarisation(polx, poly, polz);
particle->SetWeight(weight);
particle->SetUniqueID(proc);
// Increment counter
if (parentId < 0) {
mNumberOfPrimaryParticles++;
}
// Set argument variable
ntr = trackId;
// Push particle on the stack if toBeDone is set
if (toBeDone == 1) {
mStack.push(particle);
}
}
TParticle* Stack::PopNextTrack(Int_t& iTrack)
{
// If end of stack: Return empty pointer
if (mStack.empty()) {
iTrack = -1;
return NULL;
}
// If not, get next particle from stack
TParticle* thisParticle = mStack.top();
mStack.pop();
if (!thisParticle) {
iTrack = 0;
return NULL;
}
mIndexOfCurrentTrack = thisParticle->GetStatusCode();
iTrack = mIndexOfCurrentTrack;
return thisParticle;
}
TParticle* Stack::PopPrimaryForTracking(Int_t iPrim)
{
// Get the iPrimth particle from the mStack TClonesArray. This
// should be a primary (if the index is correct).
// Test for index
if (iPrim < 0 || iPrim >= mNumberOfPrimaryParticles) {
if (mLogger) {
mLogger->Fatal(MESSAGE_ORIGIN, "Stack: Primary index out of range! %i ", iPrim);
}
Fatal("Stack::PopPrimaryForTracking", "Index out of range");
}
// Return the iPrim-th TParticle from the fParticle array. This should be
// a primary.
TParticle* part = (TParticle*)mParticles->At(iPrim);
if (!(part->GetMother(0) < 0)) {
if (mLogger) {
mLogger->Fatal(MESSAGE_ORIGIN, "Stack:: Not a primary track! %i ", iPrim);
}
Fatal("Stack::PopPrimaryForTracking", "Not a primary track");
}
return part;
}
TParticle* Stack::GetCurrentTrack() const
{
TParticle* currentPart = GetParticle(mIndexOfCurrentTrack);
if (!currentPart) {
if (mLogger) {
mLogger->Warning(MESSAGE_ORIGIN, "Stack: Current track not found in stack!");
}
Warning("Stack::GetCurrentTrack", "Track not found in stack");
}
return currentPart;
}
void Stack::AddParticle(TParticle* oldPart)
{
TClonesArray& array = *mParticles;
TParticle* newPart = new (array[mIndex]) TParticle(*oldPart);
newPart->SetWeight(oldPart->GetWeight());
newPart->SetUniqueID(oldPart->GetUniqueID());
mIndex++;
}
void Stack::FillTrackArray()
{
if (mLogger) {
mLogger->Debug(MESSAGE_ORIGIN, "Stack: Filling MCTrack array...");
} else {
cout << "Stack: Filling MCTrack array..." << endl;
}
// Reset index map and number of output tracks
mIndexMap.clear();
mNumberOfEntriesInTracks = 0;
// Check tracks for selection criteria
SelectTracks();
// Loop over mParticles array and copy selected tracks
for (Int_t iPart = 0; iPart < mNumberOfEntriesInParticles; iPart++) {
mStoreIterator = mStoreMap.find(iPart);
if (mStoreIterator == mStoreMap.end()) {
if (mLogger) {
mLogger->Fatal(MESSAGE_ORIGIN, "Stack: Particle %i not found in storage map! ", iPart);
}
Fatal("Stack::FillTrackArray", "Particle not found in storage map.");
}
Bool_t store = (*mStoreIterator).second;
if (store) {
MCTrack* track = new ((*mTracks)[mNumberOfEntriesInTracks]) MCTrack(GetParticle(iPart));
mIndexMap[iPart] = mNumberOfEntriesInTracks;
// Set the number of points in the detectors for this track
for (Int_t iDet = kAliIts; iDet < kSTOPHERE; iDet++) {
pair<Int_t, Int_t> a(iPart, iDet);
track->setNumberOfPoints(iDet, mPointsMap[a]);
}
mNumberOfEntriesInTracks++;
} else {
mIndexMap[iPart] = -2;
}
}
// Map index for primary mothers
mIndexMap[-1] = -1;
// Screen output
// Print(1);
}
void Stack::UpdateTrackIndex(TRefArray* detList)
{
if (mLogger) {
mLogger->Debug(MESSAGE_ORIGIN, "Stack: Updating track indices...");
} else {
cout << "Stack: Updating track indices..." << endl;
}
Int_t nColl = 0;
// First update mother ID in MCTracks
for (Int_t i = 0; i < mNumberOfEntriesInTracks; i++) {
MCTrack* track = (MCTrack*)mTracks->At(i);
Int_t iMotherOld = track->getMotherTrackId();
mIndexIterator = mIndexMap.find(iMotherOld);
if (mIndexIterator == mIndexMap.end()) {
if (mLogger) {
mLogger->Fatal(MESSAGE_ORIGIN, "Stack: Particle index %i not found in dex map! ", iMotherOld);
}
Fatal("Stack::UpdateTrackIndex", "Particle index not found in map");
}
track->SetMotherTrackId((*mIndexIterator).second);
}
if (fDetList == 0) {
// Now iterate through all active detectors
fDetIter = detList->MakeIterator();
fDetIter->Reset();
} else {
fDetIter->Reset();
}
FairDetector* det = NULL;
while ((det = (FairDetector*)fDetIter->Next())) {
// Get hit collections from detector
Int_t iColl = 0;
TClonesArray* hitArray;
while ((hitArray = det->GetCollection(iColl++))) {
nColl++;
Int_t nPoints = hitArray->GetEntriesFast();
// Update track index for all MCPoints in the collection
for (Int_t iPoint = 0; iPoint < nPoints; iPoint++) {
FairMCPoint* point = (FairMCPoint*)hitArray->At(iPoint);
Int_t iTrack = point->GetTrackID();
mIndexIterator = mIndexMap.find(iTrack);
if (mIndexIterator == mIndexMap.end()) {
if (mLogger) {
mLogger->Fatal(MESSAGE_ORIGIN, "Stack: Particle index %i not found in index map! ", iTrack);
}
Fatal("Stack::UpdateTrackIndex", "Particle index not found in map");
}
point->SetTrackID((*mIndexIterator).second);
point->SetLink(FairLink("MCTrack", (*mIndexIterator).second));
}
} // Collections of this detector
} // List of active detectors
if (mLogger) {
mLogger->Debug(MESSAGE_ORIGIN, "...stack and %i collections updated.", nColl);
} else {
cout << "...stack and " << nColl << " collections updated." << endl;
}
}
void Stack::Reset()
{
mIndex = 0;
mIndexOfCurrentTrack = -1;
mNumberOfPrimaryParticles = mNumberOfEntriesInParticles = mNumberOfEntriesInTracks = 0;
while (!mStack.empty()) {
mStack.pop();
}
mParticles->Clear();
mTracks->Clear();
mPointsMap.clear();
}
void Stack::Register()
{
FairRootManager::Instance()->Register("MCTrack", "Stack", mTracks, kTRUE);
}
void Stack::Print(Int_t iVerbose) const
{
cout << "-I- Stack: Number of primaries = " << mNumberOfPrimaryParticles << endl;
cout << " Total number of particles = " << mNumberOfEntriesInParticles << endl;
cout << " Number of tracks in output = " << mNumberOfEntriesInTracks << endl;
if (iVerbose) {
for (Int_t iTrack = 0; iTrack < mNumberOfEntriesInTracks; iTrack++) {
((MCTrack*)mTracks->At(iTrack))->Print(iTrack);
}
}
}
void Stack::AddPoint(DetectorId detId)
{
Int_t iDet = detId;
// cout << "Add point for Detektor" << iDet << endl;
pair<Int_t, Int_t> a(mIndexOfCurrentTrack, iDet);
if (mPointsMap.find(a) == mPointsMap.end()) {
mPointsMap[a] = 1;
} else {
mPointsMap[a]++;
}
}
void Stack::AddPoint(DetectorId detId, Int_t iTrack)
{
if (iTrack < 0) {
return;
}
Int_t iDet = detId;
pair<Int_t, Int_t> a(iTrack, iDet);
if (mPointsMap.find(a) == mPointsMap.end()) {
mPointsMap[a] = 1;
} else {
mPointsMap[a]++;
}
}
Int_t Stack::GetCurrentParentTrackNumber() const
{
TParticle* currentPart = GetCurrentTrack();
if (currentPart) {
return currentPart->GetFirstMother();
} else {
return -1;
}
}
TParticle* Stack::GetParticle(Int_t trackID) const
{
if (trackID < 0 || trackID >= mNumberOfEntriesInParticles) {
if (mLogger) {
mLogger->Debug(MESSAGE_ORIGIN, "Stack: Particle index %i out of range.", trackID);
}
Fatal("Stack::GetParticle", "Index out of range");
}
return (TParticle*)mParticles->At(trackID);
}
void Stack::SelectTracks()
{
// Clear storage map
mStoreMap.clear();
// Check particles in the fParticle array
for (Int_t i = 0; i < mNumberOfEntriesInParticles; i++) {
TParticle* thisPart = GetParticle(i);
Bool_t store = kTRUE;
// Get track parameters
Int_t iMother = thisPart->GetMother(0);
TLorentzVector p;
thisPart->Momentum(p);
Double_t energy = p.E();
Double_t mass = p.M();
// Double_t mass = thisPart->GetMass();
Double_t eKin = energy - mass;
// Calculate number of points
Int_t nPoints = 0;
for (Int_t iDet = kAliIts; iDet < kSTOPHERE; iDet++) {
pair<Int_t, Int_t> a(i, iDet);
if (mPointsMap.find(a) != mPointsMap.end()) {
nPoints += mPointsMap[a];
}
}
// Check for cuts (store primaries in any case)
if (iMother < 0) {
store = kTRUE;
} else {
if (!mStoreSecondaries) {
store = kFALSE;
}
if (nPoints < mMinPoints) {
store = kFALSE;
}
if (eKin < mEnergyCut) {
store = kFALSE;
}
}
// Set storage flag
mStoreMap[i] = store;
}
// If flag is set, flag recursively mothers of selected tracks
if (mStoreMothers) {
for (Int_t i = 0; i < mNumberOfEntriesInParticles; i++) {
if (mStoreMap[i]) {
Int_t iMother = GetParticle(i)->GetMother(0);
while (iMother >= 0) {
mStoreMap[iMother] = kTRUE;
iMother = GetParticle(iMother)->GetMother(0);
}
}
}
}
}
FairGenericStack* Stack::CloneStack() const
{
return new AliceO2::Data::Stack(*this);
}
ClassImp(AliceO2::Data::Stack)