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