]>
Commit | Line | Data |
---|---|---|
9e1a0ddb | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
16 | /* | |
17 | $Log$ | |
ae9676e4 | 18 | Revision 1.21 2002/05/28 14:24:57 hristov |
19 | Correct warning messages | |
20 | ||
3d8c4988 | 21 | Revision 1.20 2002/04/30 11:47:30 morsch |
22 | KeepPhysics method called by PurifyKine added (N. Carrer, A.M.) | |
23 | ||
e44c5340 | 24 | Revision 1.19 2002/03/12 11:06:03 morsch |
25 | Add particle status code to argument list of SetTrack(..). | |
26 | ||
47c8bcbe | 27 | Revision 1.18 2002/02/20 16:14:41 hristov |
28 | fParticleBuffer points to object which doesn't belong to AliStack, do not delete it (J.Chudoba) | |
29 | ||
ba9938e6 | 30 | Revision 1.17 2001/12/05 08:51:56 hristov |
31 | The default constructor now creates no objects (thanks to r.Brun). Some corrections required by the previous changes. | |
32 | ||
be3380a3 | 33 | Revision 1.16 2001/11/20 09:27:55 hristov |
34 | Possibility to investigate a primary of not yet loaded particle (I.Hrivnacova) | |
35 | ||
f5f55563 | 36 | Revision 1.15 2001/09/04 15:10:37 hristov |
37 | Additional protection is included to avoid some problems using Hijing | |
38 | ||
1305ca2f | 39 | Revision 1.14 2001/08/30 09:44:06 hristov |
40 | VertexSource_t added to avoid the warnings | |
41 | ||
da2af24c | 42 | Revision 1.13 2001/08/29 13:31:42 morsch |
43 | Protection against (fTreeK == 0) in destructor. | |
44 | ||
cfb2cc9c | 45 | Revision 1.12 2001/07/27 13:03:13 hristov |
46 | Default Branch split level set to 99 | |
47 | ||
d0f40f23 | 48 | Revision 1.11 2001/07/27 12:34:20 jchudoba |
49 | remove the dummy argument in GetEvent method | |
50 | ||
cec14bc5 | 51 | Revision 1.10 2001/07/20 10:13:54 morsch |
52 | In Particle(Int_t) use GetEntriesFast to speed up the procedure. | |
53 | ||
911a3148 | 54 | Revision 1.9 2001/07/03 08:10:57 hristov |
55 | J.Chudoba's changes merged correctly with the HEAD | |
56 | ||
d1750cff | 57 | Revision 1.6 2001/05/31 06:59:06 fca |
58 | Clean setting and deleting of fParticleBuffer | |
59 | ||
17f4a3d9 | 60 | Revision 1.5 2001/05/30 12:18:46 hristov |
61 | Loop variables declared once | |
62 | ||
3ab6f951 | 63 | Revision 1.4 2001/05/25 07:25:20 hristov |
64 | AliStack destructor corrected (I.Hrivnacova) | |
65 | ||
81694e6a | 66 | Revision 1.3 2001/05/22 14:33:16 hristov |
67 | Minor changes | |
68 | ||
ed54bf31 | 69 | Revision 1.2 2001/05/17 05:49:39 fca |
70 | Reset pointers to daughters | |
71 | ||
76969085 | 72 | Revision 1.1 2001/05/16 14:57:22 alibrary |
73 | New files for folders and Stack | |
74 | ||
9e1a0ddb | 75 | */ |
76 | ||
77 | /////////////////////////////////////////////////////////////////////////////// | |
78 | // // | |
79 | // Particles stack class | |
80 | // // | |
81 | /////////////////////////////////////////////////////////////////////////////// | |
82 | ||
83 | #include <iostream.h> | |
84 | ||
85 | #include <TObjArray.h> | |
86 | #include <TParticle.h> | |
87 | #include <TTree.h> | |
88 | #include <TFile.h> | |
89 | #include <TFolder.h> | |
90 | #include <TROOT.h> | |
91 | ||
92 | #include "AliStack.h" | |
93 | #include "AliRun.h" | |
94 | #include "AliModule.h" | |
95 | #include "AliHit.h" | |
9e1a0ddb | 96 | |
97 | ClassImp(AliStack) | |
98 | ||
99 | //_____________________________________________________________________________ | |
100 | AliStack::AliStack(Int_t size) | |
101 | { | |
102 | // | |
103 | // Constructor | |
104 | // | |
105 | ||
106 | // Create the particles arrays | |
107 | fParticles = new TClonesArray("TParticle",1000); | |
108 | fParticleMap = new TObjArray(size); | |
17f4a3d9 | 109 | fParticleBuffer = 0; |
ccf7a81f | 110 | fNtrack = 0; |
111 | fNprimary = 0; | |
9e1a0ddb | 112 | fCurrent = -1; |
113 | fCurrentPrimary = -1; | |
ccf7a81f | 114 | fTreeK = 0; |
9e1a0ddb | 115 | } |
116 | ||
117 | ||
118 | //_____________________________________________________________________________ | |
119 | AliStack::AliStack() | |
120 | { | |
121 | // | |
122 | // Default constructor | |
123 | // | |
124 | ||
125 | // Create the particles arrays | |
be3380a3 | 126 | fParticles = 0; |
127 | fParticleMap = 0; | |
17f4a3d9 | 128 | fParticleBuffer = 0; |
ccf7a81f | 129 | fNtrack = 0; |
9e1a0ddb | 130 | fCurrent = -1; |
ccf7a81f | 131 | fNprimary = 0; |
9e1a0ddb | 132 | fCurrentPrimary = -1; |
ccf7a81f | 133 | fTreeK = 0; |
9e1a0ddb | 134 | } |
135 | ||
136 | ||
137 | //_____________________________________________________________________________ | |
138 | AliStack::~AliStack() | |
139 | { | |
140 | // | |
141 | // Destructor | |
142 | // | |
143 | ||
144 | if (fParticles) { | |
145 | fParticles->Delete(); | |
146 | delete fParticles; | |
147 | } | |
148 | delete fParticleMap; | |
cfb2cc9c | 149 | if (fTreeK) delete fTreeK; |
9e1a0ddb | 150 | } |
151 | ||
152 | // | |
153 | // public methods | |
154 | // | |
155 | ||
156 | //_____________________________________________________________________________ | |
157 | void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, | |
158 | Float_t *vpos, Float_t *polar, Float_t tof, | |
47c8bcbe | 159 | AliMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) |
9e1a0ddb | 160 | { |
161 | // | |
162 | // Load a track on the stack | |
163 | // | |
164 | // done 0 if the track has to be transported | |
165 | // 1 if not | |
166 | // parent identifier of the parent track. -1 for a primary | |
167 | // pdg particle code | |
168 | // pmom momentum GeV/c | |
169 | // vpos position | |
170 | // polar polarisation | |
171 | // tof time of flight in seconds | |
172 | // mecha production mechanism | |
173 | // ntr on output the number of the track stored | |
174 | // | |
175 | ||
176 | Float_t mass; | |
177 | const Int_t kfirstdaughter=-1; | |
178 | const Int_t klastdaughter=-1; | |
9e1a0ddb | 179 | // const Float_t tlife=0; |
180 | ||
181 | // | |
182 | // Here we get the static mass | |
183 | // For MC is ok, but a more sophisticated method could be necessary | |
184 | // if the calculated mass is required | |
185 | // also, this method is potentially dangerous if the mass | |
186 | // used in the MC is not the same of the PDG database | |
187 | // | |
188 | mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass(); | |
189 | Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+ | |
190 | pmom[1]*pmom[1]+pmom[2]*pmom[2]); | |
191 | ||
ccf7a81f | 192 | // printf("Loading mass %f ene %f No %d ip %d parent %d done %d pos %f %f %f mom %f %f %f kS %d m \n", |
193 | // mass,e,fNtrack,pdg,parent,done,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS); | |
194 | ||
9e1a0ddb | 195 | TClonesArray &particles = *fParticles; |
196 | TParticle* particle | |
197 | = new(particles[fLoadPoint++]) | |
47c8bcbe | 198 | TParticle(pdg, is, parent, -1, kfirstdaughter, klastdaughter, |
9e1a0ddb | 199 | pmom[0], pmom[1], pmom[2], e, vpos[0], vpos[1], vpos[2], tof); |
200 | particle->SetPolarisation(TVector3(polar[0],polar[1],polar[2])); | |
201 | particle->SetWeight(weight); | |
202 | particle->SetUniqueID(mech); | |
203 | if(!done) particle->SetBit(kDoneBit); | |
204 | ||
205 | ||
206 | // Declare that the daughter information is valid | |
207 | particle->SetBit(kDaughtersBit); | |
208 | // Add the particle to the stack | |
209 | fParticleMap->AddAtAndExpand(particle, fNtrack); | |
210 | ||
211 | if(parent>=0) { | |
212 | particle = (TParticle*) fParticleMap->At(parent); | |
1305ca2f | 213 | if (particle) { |
214 | particle->SetLastDaughter(fNtrack); | |
215 | if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); | |
216 | } | |
217 | else { | |
218 | printf("Error in AliStack::SetTrack: Parent %d does not exist\n",parent); | |
219 | } | |
9e1a0ddb | 220 | } |
221 | else { | |
222 | // | |
223 | // This is a primary track. Set high water mark for this event | |
224 | fHgwmk = fNtrack; | |
225 | // | |
226 | // Set also number if primary tracks | |
227 | fNprimary = fHgwmk+1; | |
228 | fCurrentPrimary++; | |
229 | } | |
230 | ntr = fNtrack++; | |
231 | } | |
232 | ||
233 | //_____________________________________________________________________________ | |
234 | void AliStack::SetTrack(Int_t done, Int_t parent, Int_t pdg, | |
235 | Double_t px, Double_t py, Double_t pz, Double_t e, | |
236 | Double_t vx, Double_t vy, Double_t vz, Double_t tof, | |
237 | Double_t polx, Double_t poly, Double_t polz, | |
47c8bcbe | 238 | AliMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) |
9e1a0ddb | 239 | { |
240 | // | |
241 | // Load a track on the stack | |
242 | // | |
243 | // done 0 if the track has to be transported | |
244 | // 1 if not | |
245 | // parent identifier of the parent track. -1 for a primary | |
246 | // pdg particle code | |
247 | // kS generation status code | |
248 | // px, py, pz momentum GeV/c | |
249 | // vx, vy, vz position | |
250 | // polar polarisation | |
251 | // tof time of flight in seconds | |
252 | // mech production mechanism | |
253 | // ntr on output the number of the track stored | |
254 | // | |
255 | // New method interface: | |
256 | // arguments were changed to be in correspondence with TParticle | |
257 | // constructor. | |
258 | // Note: the energy is not calculated from the static mass but | |
259 | // it is passed by argument e. | |
260 | ||
261 | ||
9e1a0ddb | 262 | const Int_t kFirstDaughter=-1; |
263 | const Int_t kLastDaughter=-1; | |
264 | ||
265 | TClonesArray &particles = *fParticles; | |
266 | TParticle* particle | |
267 | = new(particles[fLoadPoint++]) | |
47c8bcbe | 268 | TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter, |
9e1a0ddb | 269 | px, py, pz, e, vx, vy, vz, tof); |
270 | ||
271 | particle->SetPolarisation(polx, poly, polz); | |
272 | particle->SetWeight(weight); | |
273 | particle->SetUniqueID(mech); | |
274 | ||
275 | if(!done) particle->SetBit(kDoneBit); | |
276 | ||
277 | // Declare that the daughter information is valid | |
278 | particle->SetBit(kDaughtersBit); | |
279 | // Add the particle to the stack | |
280 | fParticleMap->AddAtAndExpand(particle, fNtrack);//CHECK!! | |
281 | ||
282 | if(parent>=0) { | |
283 | particle = (TParticle*) fParticleMap->At(parent); | |
1305ca2f | 284 | if (particle) { |
285 | particle->SetLastDaughter(fNtrack); | |
286 | if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack); | |
287 | } | |
288 | else { | |
289 | printf("Error in AliStack::SetTrack: Parent %d does not exist\n",parent); | |
290 | } | |
9e1a0ddb | 291 | } |
292 | else { | |
293 | // | |
294 | // This is a primary track. Set high water mark for this event | |
295 | fHgwmk = fNtrack; | |
296 | // | |
297 | // Set also number if primary tracks | |
298 | fNprimary = fHgwmk+1; | |
299 | fCurrentPrimary++; | |
300 | } | |
301 | ntr = fNtrack++; | |
302 | } | |
303 | ||
304 | //_____________________________________________________________________________ | |
305 | void AliStack::GetNextTrack(Int_t &mtrack, Int_t &ipart, Float_t *pmom, | |
306 | Float_t &e, Float_t *vpos, Float_t *polar, | |
307 | Float_t &tof) | |
308 | { | |
309 | // | |
310 | // Return next track from stack of particles | |
311 | // | |
312 | ||
313 | ||
314 | TParticle* track = GetNextParticle(); | |
ccf7a81f | 315 | // cout << "GetNextTrack():" << fCurrent << fNprimary << endl; |
9e1a0ddb | 316 | |
317 | if(track) { | |
318 | mtrack=fCurrent; | |
319 | ipart=track->GetPdgCode(); | |
320 | pmom[0]=track->Px(); | |
321 | pmom[1]=track->Py(); | |
322 | pmom[2]=track->Pz(); | |
323 | e =track->Energy(); | |
324 | vpos[0]=track->Vx(); | |
325 | vpos[1]=track->Vy(); | |
326 | vpos[2]=track->Vz(); | |
327 | TVector3 pol; | |
328 | track->GetPolarisation(pol); | |
329 | polar[0]=pol.X(); | |
330 | polar[1]=pol.Y(); | |
331 | polar[2]=pol.Z(); | |
332 | tof=track->T(); | |
333 | track->SetBit(kDoneBit); | |
ccf7a81f | 334 | // cout << "Filled params" << endl; |
9e1a0ddb | 335 | } |
336 | else | |
337 | mtrack=-1; | |
338 | ||
339 | // | |
340 | // stop and start timer when we start a primary track | |
341 | Int_t nprimaries = fNprimary; | |
342 | if (fCurrent >= nprimaries) return; | |
343 | if (fCurrent < nprimaries-1) { | |
344 | fTimer.Stop(); | |
345 | track=(TParticle*) fParticleMap->At(fCurrent+1); | |
346 | // track->SetProcessTime(fTimer.CpuTime()); | |
347 | } | |
348 | fTimer.Start(); | |
349 | } | |
350 | ||
351 | ||
352 | //_____________________________________________________________________________ | |
353 | void AliStack::PurifyKine() | |
354 | { | |
355 | // | |
356 | // Compress kinematic tree keeping only flagged particles | |
357 | // and renaming the particle id's in all the hits | |
358 | // | |
359 | ||
360 | TObjArray &particles = *fParticleMap; | |
361 | int nkeep=fHgwmk+1, parent, i; | |
362 | TParticle *part, *father; | |
363 | TArrayI map(particles.GetLast()+1); | |
364 | ||
365 | // Save in Header total number of tracks before compression | |
366 | ||
367 | // If no tracks generated return now | |
368 | if(fHgwmk+1 == fNtrack) return; | |
369 | ||
9e1a0ddb | 370 | // First pass, invalid Daughter information |
371 | for(i=0; i<fNtrack; i++) { | |
372 | // Preset map, to be removed later | |
373 | if(i<=fHgwmk) map[i]=i ; | |
374 | else { | |
375 | map[i] = -99; | |
e44c5340 | 376 | if((part=(TParticle*) particles.At(i))) { |
377 | // | |
378 | // Check of this track should be kept for physics reasons | |
379 | if (KeepPhysics(part)) KeepTrack(i); | |
380 | // | |
76969085 | 381 | part->ResetBit(kDaughtersBit); |
382 | part->SetFirstDaughter(-1); | |
383 | part->SetLastDaughter(-1); | |
384 | } | |
9e1a0ddb | 385 | } |
386 | } | |
387 | // Invalid daughter information for the parent of the first particle | |
388 | // generated. This may or may not be the current primary according to | |
389 | // whether decays have been recorded among the primaries | |
390 | part = (TParticle *)particles.At(fHgwmk+1); | |
391 | particles.At(part->GetFirstMother())->ResetBit(kDaughtersBit); | |
392 | // Second pass, build map between old and new numbering | |
393 | for(i=fHgwmk+1; i<fNtrack; i++) { | |
394 | if(particles.At(i)->TestBit(kKeepBit)) { | |
395 | ||
396 | // This particle has to be kept | |
397 | map[i]=nkeep; | |
398 | // If old and new are different, have to move the pointer | |
399 | if(i!=nkeep) particles[nkeep]=particles.At(i); | |
400 | part = (TParticle*) particles.At(nkeep); | |
401 | ||
402 | // as the parent is always *before*, it must be already | |
403 | // in place. This is what we are checking anyway! | |
404 | if((parent=part->GetFirstMother())>fHgwmk) | |
405 | if(map[parent]==-99) Fatal("PurifyKine","map[%d] = -99!\n",parent); | |
406 | else part->SetFirstMother(map[parent]); | |
407 | ||
408 | nkeep++; | |
409 | } | |
410 | } | |
411 | ||
412 | // Fix daughters information | |
413 | for (i=fHgwmk+1; i<nkeep; i++) { | |
414 | part = (TParticle *)particles.At(i); | |
415 | parent = part->GetFirstMother(); | |
416 | if(parent>=0) { | |
417 | father = (TParticle *)particles.At(parent); | |
418 | if(father->TestBit(kDaughtersBit)) { | |
419 | ||
420 | if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i); | |
421 | if(i>father->GetLastDaughter()) father->SetLastDaughter(i); | |
422 | } else { | |
423 | // Initialise daughters info for first pass | |
424 | father->SetFirstDaughter(i); | |
425 | father->SetLastDaughter(i); | |
426 | father->SetBit(kDaughtersBit); | |
427 | } | |
428 | } | |
429 | } | |
430 | ||
431 | // Now loop on all registered hit lists | |
432 | TList* hitLists = gAlice->GetHitLists(); | |
433 | TIter next(hitLists); | |
434 | TCollection *hitList; | |
435 | while((hitList = (TCollection*)next())) { | |
436 | TIter nexthit(hitList); | |
437 | AliHit *hit; | |
438 | while((hit = (AliHit*)nexthit())) { | |
439 | hit->SetTrack(map[hit->GetTrack()]); | |
440 | } | |
441 | } | |
442 | ||
443 | // | |
444 | // This for detectors which have a special mapping mechanism | |
445 | // for hits, such as TPC and TRD | |
446 | // | |
447 | ||
448 | TObjArray* modules = gAlice->Modules(); | |
449 | TIter nextmod(modules); | |
450 | AliModule *detector; | |
451 | while((detector = (AliModule*)nextmod())) { | |
452 | detector->RemapTrackHitIDs(map.GetArray()); | |
ae9676e4 | 453 | detector->RemapTrackReferencesIDs(map.GetArray()); |
9e1a0ddb | 454 | } |
455 | ||
456 | // Now the output bit, from fHgwmk to nkeep we write everything and we erase | |
457 | if(nkeep>fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5)); | |
458 | ||
459 | for (i=fHgwmk+1; i<nkeep; ++i) { | |
460 | fParticleBuffer = (TParticle*) particles.At(i); | |
461 | fParticleFileMap[i]=(Int_t) fTreeK->GetEntries(); | |
462 | fTreeK->Fill(); | |
17f4a3d9 | 463 | particles[i]=fParticleBuffer=0; |
9e1a0ddb | 464 | } |
465 | ||
466 | for (i=nkeep; i<fNtrack; ++i) particles[i]=0; | |
467 | ||
76969085 | 468 | Int_t toshrink = fNtrack-fHgwmk-1; |
9e1a0ddb | 469 | fLoadPoint-=toshrink; |
470 | for(i=fLoadPoint; i<fLoadPoint+toshrink; ++i) fParticles->RemoveAt(i); | |
471 | ||
472 | fNtrack=nkeep; | |
473 | fHgwmk=nkeep-1; | |
474 | // delete [] map; | |
475 | } | |
476 | ||
e44c5340 | 477 | Bool_t AliStack::KeepPhysics(TParticle* part) |
478 | { | |
479 | // | |
480 | // Some particles have to kept on the stack for reasons motivated | |
481 | // by physics analysis. Decision is put here. | |
482 | // | |
483 | Bool_t keep = kFALSE; | |
484 | // | |
485 | // Keep first-generation daughter from primaries with heavy flavor | |
486 | // | |
487 | Int_t parent = part->GetFirstMother(); | |
488 | if (parent >= 0 && parent <= fHgwmk) { | |
489 | TParticle* father = (TParticle*) Particles()->At(parent); | |
490 | Int_t kf = father->GetPdgCode(); | |
491 | kf = TMath::Abs(kf); | |
492 | Int_t kfl = kf; | |
493 | // meson ? | |
494 | if (kfl > 10) kfl/=100; | |
495 | // baryon | |
496 | if (kfl > 10) kfl/=10; | |
497 | if (kfl > 10) kfl/=10; | |
498 | if (kfl >= 4) { | |
499 | keep = kTRUE; | |
500 | } | |
501 | } | |
502 | return keep; | |
503 | } | |
504 | ||
9e1a0ddb | 505 | //_____________________________________________________________________________ |
506 | void AliStack::FinishEvent() | |
507 | { | |
508 | // | |
509 | // Write out the kinematics that was not yet filled | |
510 | // | |
511 | ||
512 | // Update event header | |
513 | ||
514 | ||
515 | if (!fTreeK) { | |
516 | // Fatal("FinishEvent", "No kinematics tree is defined."); | |
517 | // Don't panic this is a probably a lego run | |
518 | return; | |
519 | ||
520 | } | |
521 | ||
522 | CleanParents(); | |
523 | if(fTreeK->GetEntries() ==0) { | |
524 | // set the fParticleFileMap size for the first time | |
525 | fParticleFileMap.Set(fHgwmk+1); | |
526 | } | |
527 | ||
528 | Bool_t allFilled = kFALSE; | |
529 | TObject *part; | |
530 | for(Int_t i=0; i<fHgwmk+1; ++i) | |
531 | if((part=fParticleMap->At(i))) { | |
532 | fParticleBuffer = (TParticle*) part; | |
533 | fParticleFileMap[i]= (Int_t) fTreeK->GetEntries(); | |
534 | fTreeK->Fill(); | |
da2af24c | 535 | //PH (*fParticleMap)[i]=fParticleBuffer=0; |
536 | fParticleBuffer=0; | |
537 | fParticleMap->AddAt(0,i); | |
9e1a0ddb | 538 | |
539 | // When all primaries were filled no particle!=0 | |
540 | // should be left => to be removed later. | |
541 | if (allFilled) printf("Why != 0 part # %d?\n",i); | |
542 | } | |
543 | else { | |
544 | // // printf("Why = 0 part # %d?\n",i); => We know. | |
545 | // break; | |
546 | // we don't break now in order to be sure there is no | |
547 | // particle !=0 left. | |
548 | // To be removed later and replaced with break. | |
549 | if(!allFilled) allFilled = kTRUE; | |
550 | } | |
ccf7a81f | 551 | // cout << "Nof particles: " << fNtrack << endl; |
9e1a0ddb | 552 | //Reset(); |
553 | } | |
554 | ||
555 | //_____________________________________________________________________________ | |
556 | void AliStack::FlagTrack(Int_t track) | |
557 | { | |
558 | // | |
559 | // Flags a track and all its family tree to be kept | |
560 | // | |
561 | ||
562 | TParticle *particle; | |
563 | ||
564 | Int_t curr=track; | |
565 | while(1) { | |
566 | particle=(TParticle*)fParticleMap->At(curr); | |
567 | ||
568 | // If the particle is flagged the three from here upward is saved already | |
569 | if(particle->TestBit(kKeepBit)) return; | |
570 | ||
571 | // Save this particle | |
572 | particle->SetBit(kKeepBit); | |
573 | ||
574 | // Move to father if any | |
575 | if((curr=particle->GetFirstMother())==-1) return; | |
576 | } | |
577 | } | |
578 | ||
579 | //_____________________________________________________________________________ | |
ed54bf31 | 580 | void AliStack::KeepTrack(Int_t track) |
9e1a0ddb | 581 | { |
582 | // | |
583 | // Flags a track to be kept | |
584 | // | |
585 | ||
586 | fParticleMap->At(track)->SetBit(kKeepBit); | |
587 | } | |
588 | ||
589 | //_____________________________________________________________________________ | |
590 | void AliStack::Reset(Int_t size) | |
591 | { | |
592 | // | |
593 | // Resets stack | |
594 | // | |
595 | ||
596 | fNtrack=0; | |
597 | fNprimary=0; | |
598 | fHgwmk=0; | |
599 | fLoadPoint=0; | |
600 | fCurrent = -1; | |
601 | ResetArrays(size); | |
602 | } | |
603 | ||
604 | //_____________________________________________________________________________ | |
605 | void AliStack::ResetArrays(Int_t size) | |
606 | { | |
607 | // | |
608 | // Resets stack arrays | |
609 | // | |
610 | ||
be3380a3 | 611 | if (fParticles) |
612 | fParticles->Clear(); | |
613 | else | |
614 | fParticles = new TClonesArray("TParticle",1000); | |
615 | if (fParticleMap) { | |
616 | fParticleMap->Clear(); | |
617 | if (size>0) fParticleMap->Expand(size);} | |
618 | else | |
619 | fParticleMap = new TObjArray(size); | |
9e1a0ddb | 620 | } |
621 | ||
622 | //_____________________________________________________________________________ | |
ed54bf31 | 623 | void AliStack::SetHighWaterMark(Int_t nt) |
9e1a0ddb | 624 | { |
625 | // | |
626 | // Set high water mark for last track in event | |
627 | // | |
628 | ||
629 | fHgwmk = fNtrack-1; | |
630 | fCurrentPrimary=fHgwmk; | |
631 | ||
632 | // Set also number of primary tracks | |
633 | fNprimary = fHgwmk+1; | |
634 | fNtrack = fHgwmk+1; | |
635 | } | |
636 | ||
637 | //_____________________________________________________________________________ | |
638 | TParticle* AliStack::Particle(Int_t i) | |
639 | { | |
640 | // | |
641 | // Return particle with specified ID | |
642 | ||
da2af24c | 643 | //PH if(!(*fParticleMap)[i]) { |
644 | if(!fParticleMap->At(i)) { | |
911a3148 | 645 | Int_t nentries = fParticles->GetEntriesFast(); |
9e1a0ddb | 646 | // algorithmic way of getting entry index |
647 | // (primary particles are filled after secondaries) | |
648 | Int_t entry; | |
649 | if (i<fNprimary) | |
650 | entry = i+fNtrack-fNprimary; | |
651 | else | |
652 | entry = i-fNprimary; | |
653 | // check whether algorithmic way and | |
654 | // and the fParticleFileMap[i] give the same; | |
655 | // give the fatal error if not | |
656 | if (entry != fParticleFileMap[i]) { | |
657 | Fatal("Particle", | |
658 | "!! The algorithmic way and map are different: !!\n entry: %d map: %d", | |
659 | entry, fParticleFileMap[i]); | |
660 | } | |
661 | ||
662 | fTreeK->GetEntry(entry); | |
663 | new ((*fParticles)[nentries]) TParticle(*fParticleBuffer); | |
664 | fParticleMap->AddAt((*fParticles)[nentries],i); | |
665 | } | |
da2af24c | 666 | //PH return (TParticle *) (*fParticleMap)[i]; |
667 | return (TParticle *) fParticleMap->At(i); | |
9e1a0ddb | 668 | } |
669 | ||
670 | //_____________________________________________________________________________ | |
f5f55563 | 671 | Int_t AliStack::GetPrimary(Int_t id) |
9e1a0ddb | 672 | { |
673 | // | |
674 | // Return number of primary that has generated track | |
675 | // | |
676 | ||
677 | int current, parent; | |
9e1a0ddb | 678 | // |
679 | parent=id; | |
680 | while (1) { | |
681 | current=parent; | |
f5f55563 | 682 | parent=Particle(current)->GetFirstMother(); |
9e1a0ddb | 683 | if(parent<0) return current; |
684 | } | |
685 | } | |
686 | ||
687 | //_____________________________________________________________________________ | |
688 | void AliStack::DumpPart (Int_t i) const | |
689 | { | |
690 | // | |
691 | // Dumps particle i in the stack | |
692 | // | |
693 | ||
da2af24c | 694 | //PH ((TParticle*) (*fParticleMap)[i])->Print(); |
695 | ((TParticle*) fParticleMap->At(i))->Print(); | |
9e1a0ddb | 696 | } |
697 | ||
698 | //_____________________________________________________________________________ | |
699 | void AliStack::DumpPStack () | |
700 | { | |
701 | // | |
702 | // Dumps the particle stack | |
703 | // | |
704 | ||
3ab6f951 | 705 | Int_t i; |
706 | ||
9e1a0ddb | 707 | printf( |
708 | "\n\n=======================================================================\n"); | |
3ab6f951 | 709 | for (i=0;i<fNtrack;i++) |
9e1a0ddb | 710 | { |
711 | TParticle* particle = Particle(i); | |
712 | if (particle) { | |
713 | printf("-> %d ",i); particle->Print(); | |
714 | printf("--------------------------------------------------------------\n"); | |
715 | } | |
716 | else | |
717 | Warning("DumpPStack", "No particle with id %d.", i); | |
718 | } | |
719 | ||
720 | printf( | |
721 | "\n=======================================================================\n\n"); | |
722 | ||
723 | // print particle file map | |
724 | printf("\nParticle file map: \n"); | |
3ab6f951 | 725 | for (i=0; i<fNtrack; i++) |
9e1a0ddb | 726 | printf(" %d th entry: %d \n",i,fParticleFileMap[i]); |
727 | } | |
728 | ||
729 | ||
730 | //_____________________________________________________________________________ | |
731 | void AliStack::DumpLoadedStack() const | |
732 | { | |
733 | // | |
734 | // Dumps the particle in the stack | |
735 | // that are loaded in memory. | |
736 | // | |
737 | ||
738 | TObjArray &particles = *fParticleMap; | |
739 | printf( | |
740 | "\n\n=======================================================================\n"); | |
741 | for (Int_t i=0;i<fNtrack;i++) | |
742 | { | |
743 | TParticle* particle = (TParticle*) particles[i]; | |
744 | if (particle) { | |
745 | printf("-> %d ",i); particle->Print(); | |
746 | printf("--------------------------------------------------------------\n"); | |
747 | } | |
748 | else { | |
749 | printf("-> %d Particle not loaded.\n",i); | |
750 | printf("--------------------------------------------------------------\n"); | |
751 | } | |
752 | } | |
753 | printf( | |
754 | "\n=======================================================================\n\n"); | |
755 | } | |
756 | ||
757 | // | |
758 | // protected methods | |
759 | // | |
760 | ||
761 | //_____________________________________________________________________________ | |
762 | void AliStack::CleanParents() | |
763 | { | |
764 | // | |
765 | // Clean particles stack | |
766 | // Set parent/daughter relations | |
767 | // | |
768 | ||
769 | TObjArray &particles = *fParticleMap; | |
770 | TParticle *part; | |
771 | int i; | |
772 | for(i=0; i<fHgwmk+1; i++) { | |
773 | part = (TParticle *)particles.At(i); | |
774 | if(part) if(!part->TestBit(kDaughtersBit)) { | |
775 | part->SetFirstDaughter(-1); | |
776 | part->SetLastDaughter(-1); | |
777 | } | |
778 | } | |
779 | } | |
780 | ||
781 | //_____________________________________________________________________________ | |
782 | TParticle* AliStack::GetNextParticle() | |
783 | { | |
784 | // | |
785 | // Return next particle from stack of particles | |
786 | // | |
787 | ||
788 | TParticle* particle = 0; | |
789 | ||
790 | // search secondaries | |
791 | //for(Int_t i=fNtrack-1; i>=0; i--) { | |
792 | for(Int_t i=fNtrack-1; i>fHgwmk; i--) { | |
793 | particle = (TParticle*) fParticleMap->At(i); | |
794 | if ((particle) && (!particle->TestBit(kDoneBit))) { | |
795 | fCurrent=i; | |
796 | //cout << "GetNextParticle() - secondary " | |
797 | // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; | |
798 | return particle; | |
799 | } | |
800 | } | |
801 | ||
802 | // take next primary if all secondaries were done | |
803 | while (fCurrentPrimary>=0) { | |
804 | fCurrent = fCurrentPrimary; | |
805 | particle = (TParticle*) fParticleMap->At(fCurrentPrimary--); | |
806 | if ((particle) && (!particle->TestBit(kDoneBit))) { | |
807 | //cout << "GetNextParticle() - primary " | |
808 | // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; | |
809 | return particle; | |
810 | } | |
811 | } | |
812 | ||
813 | // nothing to be tracked | |
814 | fCurrent = -1; | |
815 | //cout << "GetNextParticle() - none " | |
816 | // << fNtrack << " " << fHgwmk << " " << fCurrent << endl; | |
817 | return particle; | |
818 | } | |
819 | ||
820 | //__________________________________________________________________________________________ | |
821 | void AliStack::MakeTree(Int_t event, const char *file) | |
822 | { | |
823 | // | |
824 | // Make Kine tree and creates branch for writing particles | |
825 | // | |
826 | TBranch *branch=0; | |
827 | // Make Kinematics Tree | |
828 | char hname[30]; | |
9e1a0ddb | 829 | if (!fTreeK) { |
830 | sprintf(hname,"TreeK%d",event); | |
831 | fTreeK = new TTree(hname,"Kinematics"); | |
832 | // Create a branch for particles | |
d0f40f23 | 833 | branch = fTreeK->Branch("Particles", "TParticle", &fParticleBuffer, 4000); |
9e1a0ddb | 834 | fTreeK->Write(0,TObject::kOverwrite); |
835 | } | |
836 | } | |
837 | ||
81694e6a | 838 | //_____________________________________________________________________________ |
9e1a0ddb | 839 | void AliStack::BeginEvent(Int_t event) |
840 | { | |
841 | // start a new event | |
842 | // | |
843 | // | |
844 | fNprimary = 0; | |
845 | fNtrack = 0; | |
846 | ||
847 | char hname[30]; | |
848 | if(fTreeK) { | |
849 | fTreeK->Reset(); | |
850 | sprintf(hname,"TreeK%d",event); | |
851 | fTreeK->SetName(hname); | |
852 | } | |
853 | } | |
854 | ||
81694e6a | 855 | //_____________________________________________________________________________ |
9e1a0ddb | 856 | void AliStack::FinishRun() |
857 | { | |
858 | // Clean TreeK information | |
859 | if (fTreeK) { | |
860 | delete fTreeK; fTreeK = 0; | |
861 | } | |
862 | } | |
863 | ||
cec14bc5 | 864 | Bool_t AliStack::GetEvent(Int_t event) |
9e1a0ddb | 865 | { |
866 | // | |
867 | // Get new event from TreeK | |
ccf7a81f | 868 | |
9e1a0ddb | 869 | // Reset/Create the particle stack |
870 | if (fTreeK) delete fTreeK; | |
871 | ||
872 | // Get Kine Tree from file | |
873 | char treeName[20]; | |
874 | sprintf(treeName,"TreeK%d",event); | |
875 | fTreeK = (TTree*)gDirectory->Get(treeName); | |
876 | ||
ccf7a81f | 877 | if (fTreeK) |
878 | fTreeK->SetBranchAddress("Particles", &fParticleBuffer); | |
879 | else { | |
3d8c4988 | 880 | // Error("GetEvent","cannot find Kine Tree for event:%d\n",event); |
881 | Warning("GetEvent","cannot find Kine Tree for event:%d\n",event); | |
ccf7a81f | 882 | return kFALSE; |
883 | } | |
884 | // printf("\n primaries %d", fNprimary); | |
885 | // printf("\n tracks %d", fNtrack); | |
886 | ||
9e1a0ddb | 887 | Int_t size = (Int_t)fTreeK->GetEntries(); |
888 | ResetArrays(size); | |
ccf7a81f | 889 | return kTRUE; |
9e1a0ddb | 890 | } |