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