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