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