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