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