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