fbf97c7e11875caff8881c2e844b0a6926602daf
[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 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 //  Particles stack class                                                    //
21 //  Implements the TMCVirtualStack of the Virtual Monte Carlo                //
22 //  Holds the particles transported during simulation                        //
23 //  Is used to compare results of reconstruction with simulation             //
24 //  Author A.Morsch                                                          //
25 //                                                                           //
26 ///////////////////////////////////////////////////////////////////////////////
27
28  
29 #include <TObjArray.h>
30 #include <TParticle.h>
31 #include <TParticlePDG.h>
32 #include <TTree.h>
33
34 #include "AliLog.h"
35 #include "AliHit.h"
36 #include "AliModule.h"
37 #include "AliRun.h"
38 #include "AliMC.h"
39 #include "AliRunLoader.h"
40 #include "AliStack.h"
41
42 ClassImp(AliStack)
43
44 //_______________________________________________________________________
45 AliStack::AliStack():
46   fParticles(0),
47   fParticleMap(0),
48   fParticleFileMap(0),
49   fParticleBuffer(0),
50   fTreeK(0),
51   fNtrack(0),
52   fNprimary(0),
53   fCurrent(-1),
54   fCurrentPrimary(-1),
55   fHgwmk(0),
56   fLoadPoint(0),
57   fEventFolderName(AliConfig::GetDefaultEventFolderName())
58 {
59   //
60   // Default constructor
61   //
62 }
63
64 //_______________________________________________________________________
65 AliStack::AliStack(Int_t size, const char* evfoldname):
66   fParticles(new TClonesArray("TParticle",1000)),
67   fParticleMap(new TObjArray(size)),
68   fParticleFileMap(0),
69   fParticleBuffer(0),
70   fTreeK(0),
71   fNtrack(0),
72   fNprimary(0),
73   fCurrent(-1),
74   fCurrentPrimary(-1),
75   fHgwmk(0),
76   fLoadPoint(0),
77   fEventFolderName(evfoldname)
78 {
79   //
80   //  Constructor
81   //
82 }
83
84 //_______________________________________________________________________
85 AliStack::AliStack(const AliStack& st):
86   TVirtualMCStack(st),
87   fParticles(0),
88   fParticleMap(0),
89   fParticleFileMap(0),
90   fParticleBuffer(0),
91   fTreeK(0),
92   fNtrack(0),
93   fNprimary(0),
94   fCurrent(-1),
95   fCurrentPrimary(-1),
96   fHgwmk(0),
97   fLoadPoint(0)
98 {
99   //
100   // Copy constructor
101   //
102   st.Copy(*this);
103 }
104
105 //_______________________________________________________________________
106 void AliStack::Copy(TObject&) const
107 {
108   AliFatal("Not implemented!");
109 }
110
111 //_______________________________________________________________________
112 AliStack::~AliStack()
113 {
114   //
115   // Destructor
116   //
117   
118   if (fParticles) {
119     fParticles->Delete();
120     delete fParticles;
121   }
122   delete fParticleMap;
123   //PH???  delete fTreeK;
124 }
125
126 //
127 // public methods
128 //
129
130 //_____________________________________________________________________________
131 void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
132                         Float_t *vpos, Float_t *polar, Float_t tof,
133                         TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
134
135   //
136   // Load a track on the stack
137   //
138   // done     0 if the track has to be transported
139   //          1 if not
140   // parent   identifier of the parent track. -1 for a primary
141   // pdg    particle code
142   // pmom     momentum GeV/c
143   // vpos     position 
144   // polar    polarisation 
145   // tof      time of flight in seconds
146   // mecha    production mechanism
147   // ntr      on output the number of the track stored
148   //
149
150   //  const Float_t tlife=0;
151   
152   //
153   // Here we get the static mass
154   // For MC is ok, but a more sophisticated method could be necessary
155   // if the calculated mass is required
156   // also, this method is potentially dangerous if the mass
157   // used in the MC is not the same of the PDG database
158   //
159     TParticlePDG* pmc =  TDatabasePDG::Instance()->GetParticle(pdg);
160     if (pmc) {
161         Float_t 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
169         PushTrack(done, parent, pdg, pmom[0], pmom[1], pmom[2], e,
170                  vpos[0], vpos[1], vpos[2], tof, polar[0], polar[1], polar[2],
171                  mech, ntr, weight, is);
172     } else {
173         AliWarning(Form("Particle type %d not defined in PDG Database !", pdg));
174         AliWarning("Particle skipped !");
175     }
176 }
177
178 //_____________________________________________________________________________
179 void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg,
180                       Double_t px, Double_t py, Double_t pz, Double_t e,
181                       Double_t vx, Double_t vy, Double_t vz, Double_t tof,
182                       Double_t polx, Double_t poly, Double_t polz,
183                       TMCProcess mech, Int_t &ntr, Double_t weight, Int_t is)
184
185   //
186   // Load a track on the stack
187   //
188   // done        0 if the track has to be transported
189   //             1 if not
190   // parent      identifier of the parent track. -1 for a primary
191   // pdg         particle code
192   // kS          generation status code
193   // px, py, pz  momentum GeV/c
194   // vx, vy, vz  position 
195   // polar       polarisation 
196   // tof         time of flight in seconds
197   // mech        production mechanism
198   // ntr         on output the number of the track stored
199   //    
200   // New method interface: 
201   // arguments were changed to be in correspondence with TParticle
202   // constructor.
203   // Note: the energy is not calculated from the static mass but
204   // it is passed by argument e.
205
206
207   const Int_t kFirstDaughter=-1;
208   const Int_t kLastDaughter=-1;
209
210
211   TClonesArray &particles = *fParticles;
212   TParticle* particle
213     = new(particles[fLoadPoint++]) 
214       TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter,
215                 px, py, pz, e, vx, vy, vz, tof);
216    
217   particle->SetPolarisation(polx, poly, polz);
218   particle->SetWeight(weight);
219   particle->SetUniqueID(mech);
220
221   if(!done) particle->SetBit(kDoneBit);
222
223   //  Declare that the daughter information is valid
224   particle->SetBit(kDaughtersBit);
225   //  Add the particle to the stack
226   
227   
228   fParticleMap->AddAtAndExpand(particle, fNtrack);//CHECK!!
229
230   if(parent>=0) {
231     particle = dynamic_cast<TParticle*>(fParticleMap->At(parent));
232     if (particle) {
233       particle->SetLastDaughter(fNtrack);
234       if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack);
235     }
236     else {
237       AliError(Form("Parent %d does not exist",parent));
238     }
239   } 
240   else { 
241     //
242     // This is a primary track. Set high water mark for this event
243     fHgwmk = fNtrack;
244     //
245     // Set also number if primary tracks
246     fNprimary = fHgwmk+1;
247     fCurrentPrimary++;
248   }
249   ntr = fNtrack++;
250 }
251
252 //_____________________________________________________________________________
253 TParticle*  AliStack::PopNextTrack(Int_t& itrack)
254 {
255   //
256   // Returns next track from stack of particles
257   //
258   
259
260   TParticle* track = GetNextParticle();
261
262   if (track) {
263     itrack = fCurrent;
264     track->SetBit(kDoneBit);
265   }
266   else
267     itrack = -1;
268   
269   fCurrentTrack = track;
270   return track;
271 }
272
273 //_____________________________________________________________________________
274 TParticle*  AliStack::PopPrimaryForTracking(Int_t i)
275 {
276   //
277   // Returns i-th primary particle if it is flagged to be tracked,
278   // 0 otherwise
279   //
280   
281   TParticle* particle = Particle(i);
282   
283   if (!particle->TestBit(kDoneBit))
284     return particle;
285   else
286     return 0;
287 }      
288
289 //_____________________________________________________________________________
290 void AliStack::PurifyKine()
291 {
292   //
293   // Compress kinematic tree keeping only flagged particles
294   // and renaming the particle id's in all the hits
295   //
296
297   TObjArray &particles = *fParticleMap;
298   int nkeep=fHgwmk+1, parent, i;
299   TParticle *part, *father;
300   TArrayI map(particles.GetLast()+1);
301
302   // Save in Header total number of tracks before compression
303   // If no tracks generated return now
304   if(fHgwmk+1 == fNtrack) return;
305
306   // First pass, invalid Daughter information
307
308   for(i=0; i<fNtrack; i++) {
309     // Preset map, to be removed later
310     if(i<=fHgwmk) map[i]=i ; 
311     else {
312       map[i] = -99;
313       if((part=dynamic_cast<TParticle*>(particles.At(i)))) {
314 //
315 //        Check of this track should be kept for physics reasons 
316           if (KeepPhysics(part)) KeepTrack(i);
317 //
318           part->ResetBit(kDaughtersBit);
319           part->SetFirstDaughter(-1);
320           part->SetLastDaughter(-1);
321       }
322     }
323   }
324   // Invalid daughter information for the parent of the first particle
325   // generated. This may or may not be the current primary according to
326   // whether decays have been recorded among the primaries
327   part = dynamic_cast<TParticle*>(particles.At(fHgwmk+1));
328   particles.At(part->GetFirstMother())->ResetBit(kDaughtersBit);
329   // Second pass, build map between old and new numbering
330   for(i=fHgwmk+1; i<fNtrack; i++) {
331     if(particles.At(i)->TestBit(kKeepBit)) {
332       
333       // This particle has to be kept
334       map[i]=nkeep;
335       // If old and new are different, have to move the pointer
336       if(i!=nkeep) particles[nkeep]=particles.At(i);
337       part = dynamic_cast<TParticle*>(particles.At(nkeep));
338       
339       // as the parent is always *before*, it must be already
340       // in place. This is what we are checking anyway!
341       if((parent=part->GetFirstMother())>fHgwmk) 
342         if(map[parent]==-99) Fatal("PurifyKine","map[%d] = -99!\n",parent);
343         else part->SetFirstMother(map[parent]);
344
345       nkeep++;
346     }
347   }
348   
349   // Fix daughters information
350   for (i=fHgwmk+1; i<nkeep; i++) {
351     part = dynamic_cast<TParticle*>(particles.At(i));
352     parent = part->GetFirstMother();
353     if(parent>=0) {
354       father = dynamic_cast<TParticle*>(particles.At(parent));
355       if(father->TestBit(kDaughtersBit)) {
356       
357         if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i);
358         if(i>father->GetLastDaughter())  father->SetLastDaughter(i);
359       } else {
360         // Initialise daughters info for first pass
361         father->SetFirstDaughter(i);
362         father->SetLastDaughter(i);
363         father->SetBit(kDaughtersBit);
364       }
365     }
366   }
367   
368   // Now loop on all registered hit lists
369   TList* hitLists = gAlice->GetMCApp()->GetHitLists();
370   TIter next(hitLists);
371   TCollection *hitList;
372   
373   while((hitList = dynamic_cast<TCollection*>(next()))) {
374     TIter nexthit(hitList);
375     AliHit *hit;
376     while((hit = dynamic_cast<AliHit*>(nexthit()))) {
377         hit->SetTrack(map[hit->GetTrack()]);
378     }
379   }
380
381   // 
382   // This for detectors which have a special mapping mechanism
383   // for hits, such as TPC and TRD
384   //
385   
386    TObjArray* modules = gAlice->Modules();
387    TIter nextmod(modules);
388    AliModule *detector;
389    while((detector = dynamic_cast<AliModule*>(nextmod()))) {
390      detector->RemapTrackHitIDs(map.GetArray());
391      detector->RemapTrackReferencesIDs(map.GetArray());
392    }
393    //
394    gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray());
395
396    // Now the output bit, from fHgwmk to nkeep we write everything and we erase
397    if(nkeep>fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5));
398
399    for (i=fHgwmk+1; i<nkeep; ++i) {
400      fParticleBuffer = dynamic_cast<TParticle*>(particles.At(i));
401      fParticleFileMap[i]=static_cast<Int_t>(TreeK()->GetEntries());
402      TreeK()->Fill();
403      particles[i]=fParticleBuffer=0;
404     }
405
406    for (i=nkeep; i<fNtrack; ++i) particles[i]=0;
407
408    Int_t toshrink = fNtrack-fHgwmk-1;
409    fLoadPoint-=toshrink;
410
411    
412    for(i=fLoadPoint; i<fLoadPoint+toshrink; ++i) fParticles->RemoveAt(i);
413
414    fNtrack=nkeep;
415    fHgwmk=nkeep-1;
416 }
417
418 void AliStack::ReorderKine()
419 {
420 //
421 // In some transport code children might not come in a continuous sequence.
422 // In this case the stack  has  to  be reordered in order to establish the 
423 // mother daughter relation using index ranges.
424 //    
425   if(fHgwmk+1 == fNtrack) return;
426
427   //
428   // Howmany secondaries have been produced ?
429   Int_t nNew = fNtrack - fHgwmk - 1;
430   
431   if (nNew > 0) {
432       Int_t i, j;
433       TObjArray &particles = *fParticleMap;
434       TArrayI map1(nNew);
435       //
436       // Copy pointers to temporary array
437       TParticle** tmp = new TParticle*[nNew];
438       
439       for (i = 0; i < nNew; i++) {
440           if (particles.At(fHgwmk + 1 + i)) {
441 //            tmp[i] = (TParticle*) (particles.At(fHgwmk + 1 + i))->Clone();
442               tmp[i] = (TParticle*) (particles.At(fHgwmk + 1 + i));
443
444 //            if (((TParticle*) (particles.At(fHgwmk + 1 + i)))->TestBit(kKeepBit))
445 //                tmp[i]->SetBit(kKeepBit);
446 //            if (((TParticle*) (particles.At(fHgwmk + 1 + i)))->TestBit(kDoneBit))
447 //                tmp[i]->SetBit(kDoneBit);
448           } else {
449               tmp[i] = 0x0;
450           }
451           map1[i] = -99;
452       }
453   
454       
455       //
456       // Reset  LoadPoint 
457       // 
458       fLoadPoint = fHgwmk + 1;
459       //
460       // Re-Push particles into stack 
461       // The outer loop is over parents, the inner over children.
462       // -1 refers to the primary particle
463       //
464       for (i = -1; i < nNew - 1; i++) {
465           Int_t ipa;
466           TParticle* parP;
467           if (i == -1) {
468               ipa  = tmp[0]->GetFirstMother();
469               parP =dynamic_cast<TParticle*>(particles.At(ipa));
470           } else {
471               ipa = (fHgwmk + 1 + i);
472               // Skip deleted particles
473               if (!tmp[i])                          continue;
474               // Skip particles without children
475               if (tmp[i]->GetFirstDaughter() == -1) continue;
476               parP = tmp[i];
477           }
478           // Reset daughter information
479           parP->SetFirstDaughter(-1);
480           parP->SetLastDaughter(-1);
481           for (j = i + 1; j < nNew; j++) {
482               // Skip deleted particles
483               if (!tmp[j])        continue;
484               // Skip particles already handled
485               if (map1[j] != -99) continue;
486               Int_t jpa = tmp[j]->GetFirstMother();
487               // Check if daughter of current parent
488               if (jpa == ipa) {
489                   particles[fLoadPoint] = tmp[j];
490                   // Re-establish daughter information
491                   parP->SetLastDaughter(fLoadPoint);
492                   if (parP->GetFirstDaughter() == -1) parP->SetFirstDaughter(fLoadPoint);
493                   // Set Mother information
494                   if (i != -1) {
495                       tmp[j]->SetFirstMother(map1[i]);
496                   } 
497                   // Build the map
498                   map1[j] = fLoadPoint;
499                   // Increase load point
500                   fLoadPoint++;
501               }
502           } // children
503       } // parents
504
505       delete[] tmp;
506
507       //
508       // Build map for remapping of hits
509       // 
510       TArrayI map(fNtrack);
511       for (i = 0; i < fNtrack; i ++) {
512           if (i <= fHgwmk) {
513               map[i] = i;
514           } else{
515               map[i] = map1[i - fHgwmk -1];
516           }
517       }
518       
519       // Now loop on all registered hit lists
520       
521       TList* hitLists = gAlice->GetMCApp()->GetHitLists();
522       TIter next(hitLists);
523       TCollection *hitList;
524       
525       while((hitList = dynamic_cast<TCollection*>(next()))) {
526           TIter nexthit(hitList);
527           AliHit *hit;
528           while((hit = dynamic_cast<AliHit*>(nexthit()))) {
529               hit->SetTrack(map[hit->GetTrack()]);
530           }
531       }
532   
533   // 
534   // This for detectors which have a special mapping mechanism
535   // for hits, such as TPC and TRD
536   //
537   
538       TObjArray* modules = gAlice->Modules();
539       TIter nextmod(modules);
540       AliModule *detector;
541       while((detector = dynamic_cast<AliModule*>(nextmod()))) {
542           detector->RemapTrackHitIDs(map.GetArray());
543           detector->RemapTrackReferencesIDs(map.GetArray());
544       }
545       gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray());
546   } // new particles poduced
547 }
548
549 Bool_t AliStack::KeepPhysics(TParticle* part)
550 {
551     //
552     // Some particles have to kept on the stack for reasons motivated
553     // by physics analysis. Decision is put here.
554     //
555     Bool_t keep = kFALSE;
556     //
557     // Keep first-generation daughter from primaries with heavy flavor 
558     //
559     Int_t parent = part->GetFirstMother();
560     if (parent >= 0 && parent <= fHgwmk) {
561         TParticle* father = dynamic_cast<TParticle*>(Particles()->At(parent));
562         Int_t kf = father->GetPdgCode();
563         kf = TMath::Abs(kf);
564         Int_t kfl = kf;
565         // meson ?
566         if  (kfl > 10) kfl/=100;
567         // baryon
568         if (kfl > 10) kfl/=10;
569         if (kfl > 10) kfl/=10;
570         if (kfl >= 4) {
571             keep = kTRUE;
572         }
573     }
574     return keep;
575 }
576
577 //_____________________________________________________________________________
578 void AliStack::FinishEvent()
579 {
580 //
581 // Write out the kinematics that was not yet filled
582 //
583   
584 // Update event header
585
586
587   if (!TreeK()) {
588 //    Fatal("FinishEvent", "No kinematics tree is defined.");
589 //    Don't panic this is a probably a lego run
590       return;
591   }  
592   
593   CleanParents();
594   if(TreeK()->GetEntries() ==0) {
595     // set the fParticleFileMap size for the first time
596     fParticleFileMap.Set(fHgwmk+1);
597   }
598
599   Bool_t allFilled = kFALSE;
600   TObject *part;
601   for(Int_t i=0; i<fHgwmk+1; ++i) 
602     if((part=fParticleMap->At(i))) {
603       fParticleBuffer = dynamic_cast<TParticle*>(part);
604       fParticleFileMap[i]= static_cast<Int_t>(TreeK()->GetEntries());
605       TreeK()->Fill();
606       //PH      (*fParticleMap)[i]=fParticleBuffer=0;      
607       fParticleBuffer=0;      
608       fParticleMap->AddAt(0,i);      
609       
610       // When all primaries were filled no particle!=0
611       // should be left => to be removed later.
612       if (allFilled) AliWarning(Form("Why != 0 part # %d?\n",i));
613     }
614     else 
615      {
616       // // printf("Why = 0 part # %d?\n",i); => We know.
617       // break;
618       // we don't break now in order to be sure there is no
619       // particle !=0 left.
620       // To be removed later and replaced with break.
621        if(!allFilled) allFilled = kTRUE;
622      } 
623
624 //_____________________________________________________________________________
625
626 void AliStack::FlagTrack(Int_t track)
627 {
628   //
629   // Flags a track and all its family tree to be kept
630   //
631   
632   TParticle *particle;
633
634   Int_t curr=track;
635   while(1) {
636     particle=dynamic_cast<TParticle*>(fParticleMap->At(curr));
637     
638     // If the particle is flagged the three from here upward is saved already
639     if(particle->TestBit(kKeepBit)) return;
640     
641     // Save this particle
642     particle->SetBit(kKeepBit);
643     
644     // Move to father if any
645     if((curr=particle->GetFirstMother())==-1) return;
646   }
647 }
648  
649 //_____________________________________________________________________________
650 void AliStack::KeepTrack(Int_t track)
651
652   //
653   // Flags a track to be kept
654   //
655   
656   fParticleMap->At(track)->SetBit(kKeepBit);
657 }
658
659 //_____________________________________________________________________________
660 void  AliStack::Reset(Int_t size) 
661 {
662   //
663   // Resets stack
664   //
665   
666   fNtrack=0;
667   fNprimary=0;
668   fHgwmk=0;
669   fLoadPoint=0;
670   fCurrent = -1;
671   fTreeK = 0x0;
672   ResetArrays(size);
673 }
674
675 //_____________________________________________________________________________
676 void  AliStack::ResetArrays(Int_t size) 
677 {
678   //
679   // Resets stack arrays
680   //
681
682   if (fParticles) 
683     fParticles->Clear();
684   else
685     fParticles = new TClonesArray("TParticle",1000);
686   if (fParticleMap) {
687     fParticleMap->Clear();
688     if (size>0) fParticleMap->Expand(size);}
689   else
690     fParticleMap = new TObjArray(size);
691 }
692
693 //_____________________________________________________________________________
694 void AliStack::SetHighWaterMark(Int_t)
695 {
696   //
697   // Set high water mark for last track in event
698   //
699
700     
701   fHgwmk = fNtrack-1;
702   fCurrentPrimary=fHgwmk;
703   // Set also number of primary tracks
704   fNprimary = fHgwmk+1;
705   fNtrack   = fHgwmk+1;      
706 }
707
708 //_____________________________________________________________________________
709 TParticle* AliStack::Particle(Int_t i)
710 {
711   //
712   // Return particle with specified ID
713   
714   //PH  if(!(*fParticleMap)[i]) {
715   if(!fParticleMap->At(i)) {
716     Int_t nentries = fParticles->GetEntriesFast();
717     // algorithmic way of getting entry index
718     // (primary particles are filled after secondaries)
719     Int_t entry = TreeKEntry(i);
720     // check whether algorithmic way and 
721     // and the fParticleFileMap[i] give the same;
722     // give the fatal error if not
723     if (entry != fParticleFileMap[i]) {
724       AliFatal(Form(
725         "!! The algorithmic way and map are different: !!\n entry: %d map: %d",
726         entry, fParticleFileMap[i])); 
727     } 
728       
729     TreeK()->GetEntry(entry);
730     new ((*fParticles)[nentries]) TParticle(*fParticleBuffer);
731     fParticleMap->AddAt((*fParticles)[nentries],i);
732   }
733   //PH  return dynamic_cast<TParticle *>((*fParticleMap)[i]);
734   return dynamic_cast<TParticle*>(fParticleMap->At(i));
735 }
736
737 //_____________________________________________________________________________
738 TParticle* AliStack::ParticleFromTreeK(Int_t id) const
739 {
740 // 
741 // return pointer to TParticle with label id
742 //
743   Int_t entry;
744   if ((entry = TreeKEntry(id)) < 0) return 0;
745   if (fTreeK->GetEntry(entry)<=0) return 0;
746   return fParticleBuffer;
747 }
748
749 //_____________________________________________________________________________
750 Int_t AliStack::TreeKEntry(Int_t id) const 
751 {
752 //
753 // return entry number in the TreeK for particle with label id
754 // return negative number if label>fNtrack
755 //
756   Int_t entry;
757   if (id<fNprimary)
758     entry = id+fNtrack-fNprimary;
759   else 
760     entry = id-fNprimary;
761   return entry;
762 }
763
764 //_____________________________________________________________________________
765 Int_t AliStack::GetCurrentParentTrackNumber() const
766 {
767   //
768   // Return number of the parent of the current track
769   //
770   
771   TParticle* current = (TParticle*)fParticleMap->At(fCurrent);
772
773   if (current) 
774     return current->GetFirstMother();
775   else {
776     AliWarning("Current track not found in the stack");
777     return -1;
778   }  
779 }
780  
781 //_____________________________________________________________________________
782 Int_t AliStack::GetPrimary(Int_t id)
783 {
784   //
785   // Return number of primary that has generated track
786   //
787   
788   int current, parent;
789   //
790   parent=id;
791   while (1) {
792     current=parent;
793     parent=Particle(current)->GetFirstMother();
794     if(parent<0) return current;
795   }
796 }
797  
798 //_____________________________________________________________________________
799 void AliStack::DumpPart (Int_t i) const
800 {
801   //
802   // Dumps particle i in the stack
803   //
804   
805   //PH  dynamic_cast<TParticle*>((*fParticleMap)[i])->Print();
806   dynamic_cast<TParticle*>(fParticleMap->At(i))->Print();
807 }
808
809 //_____________________________________________________________________________
810 void AliStack::DumpPStack ()
811 {
812   //
813   // Dumps the particle stack
814   //
815
816   Int_t i;
817
818   printf("\n\n=======================================================================\n");
819   for (i=0;i<fNtrack;i++) 
820     {
821       TParticle* particle = Particle(i);
822       if (particle) {
823         printf("-> %d ",i); particle->Print();
824         printf("--------------------------------------------------------------\n");
825       }
826       else 
827         Warning("DumpPStack", "No particle with id %d.", i); 
828     }    
829
830   printf("\n=======================================================================\n\n");
831   
832   // print  particle file map
833   printf("\nParticle file map: \n");
834   for (i=0; i<fNtrack; i++) 
835       printf("   %d th entry: %d \n",i,fParticleFileMap[i]);
836 }
837
838
839 //_____________________________________________________________________________
840 void AliStack::DumpLoadedStack() const
841 {
842   //
843   // Dumps the particle in the stack
844   // that are loaded in memory.
845   //
846
847   TObjArray &particles = *fParticleMap;
848   printf(
849          "\n\n=======================================================================\n");
850   for (Int_t i=0;i<fNtrack;i++) 
851     {
852       TParticle* particle = dynamic_cast<TParticle*>(particles[i]);
853       if (particle) {
854         printf("-> %d ",i); particle->Print();
855         printf("--------------------------------------------------------------\n");
856       }
857       else {    
858         printf("-> %d  Particle not loaded.\n",i);
859         printf("--------------------------------------------------------------\n");
860       } 
861     }
862   printf(
863          "\n=======================================================================\n\n");
864 }
865
866 //
867 // protected methods
868 //
869
870 //_____________________________________________________________________________
871 void AliStack::CleanParents()
872 {
873   //
874   // Clean particles stack
875   // Set parent/daughter relations
876   //
877   
878   TObjArray &particles = *fParticleMap;
879   TParticle *part;
880   int i;
881   for(i=0; i<fHgwmk+1; i++) {
882     part = dynamic_cast<TParticle*>(particles.At(i));
883     if(part) if(!part->TestBit(kDaughtersBit)) {
884       part->SetFirstDaughter(-1);
885       part->SetLastDaughter(-1);
886     }
887   }
888 }
889
890 //_____________________________________________________________________________
891 TParticle* AliStack::GetNextParticle()
892 {
893   //
894   // Return next particle from stack of particles
895   //
896   
897   TParticle* particle = 0;
898   
899   // search secondaries
900   //for(Int_t i=fNtrack-1; i>=0; i--) {
901   for(Int_t i=fNtrack-1; i>fHgwmk; i--) {
902       particle = dynamic_cast<TParticle*>(fParticleMap->At(i));
903       if ((particle) && (!particle->TestBit(kDoneBit))) {
904           fCurrent=i;    
905           return particle;
906       }   
907   }    
908
909   // take next primary if all secondaries were done
910   while (fCurrentPrimary>=0) {
911       fCurrent = fCurrentPrimary;    
912       particle = dynamic_cast<TParticle*>(fParticleMap->At(fCurrentPrimary--));
913       if ((particle) && (!particle->TestBit(kDoneBit))) {
914           return particle;
915       } 
916   }
917   
918   // nothing to be tracked
919   fCurrent = -1;
920  
921   
922   return particle;  
923 }
924 //__________________________________________________________________________________________
925
926 TTree* AliStack::TreeK()
927 {
928 //returns TreeK
929   if (fTreeK)
930    {
931      return fTreeK;
932    }
933   else
934    {
935      AliRunLoader *rl = AliRunLoader::GetRunLoader(fEventFolderName);
936      if (rl == 0x0)
937       {
938         AliFatal(Form("Can not get RunLoader from event folder named %s",fEventFolderName.Data()));
939         return 0x0;//pro forma
940       }
941     fTreeK = rl->TreeK();
942     if ( fTreeK )
943      {
944       ConnectTree();
945      }
946     else
947      {
948       //don't panic - could be Lego
949       AliWarning(Form("Can not get TreeK from RL. Ev. Folder is %s",fEventFolderName.Data()));
950      }
951    }
952   return fTreeK;//never reached
953 }
954 //__________________________________________________________________________________________
955
956 void AliStack::ConnectTree()
957 {
958 //
959 //  Creates branch for writing particles
960 //
961   AliDebug(1, "Connecting TreeK");
962   if (fTreeK == 0x0)
963    {
964     if (TreeK() == 0x0)
965      {
966       AliFatal("Parameter is NULL");//we don't like such a jokes
967       return;
968      }
969     return;//in this case TreeK() calls back this method (ConnectTree) 
970            //tree after setting fTreeK, the rest was already executed
971            //it is safe to return now
972    }
973
974  //  Create a branch for particles   
975   
976   AliDebug(2, Form("Tree name is %s",fTreeK->GetName()));
977    
978   if (fTreeK->GetDirectory())
979    {
980      AliDebug(2, Form("and dir is %s",fTreeK->GetDirectory()->GetName()));
981    }    
982   else
983     AliWarning("DIR IS NOT SET !!!");
984   
985   TBranch *branch=fTreeK->GetBranch(AliRunLoader::GetKineBranchName());
986   if(branch == 0x0)
987    {
988     branch = fTreeK->Branch(AliRunLoader::GetKineBranchName(), "TParticle", &fParticleBuffer, 4000);
989     AliDebug(2, "Creating Branch in Tree");
990    }  
991   else
992    {
993     AliDebug(2, "Branch Found in Tree");
994     branch->SetAddress(&fParticleBuffer);
995    }
996   if (branch->GetDirectory())
997    {
998     AliDebug(1, Form("Branch Dir Name is %s",branch->GetDirectory()->GetName()));
999    } 
1000   else
1001     AliWarning("Branch Dir is NOT SET");
1002 }
1003 //__________________________________________________________________________________________
1004
1005
1006 void AliStack::BeginEvent()
1007 {
1008 // start a new event
1009  Reset();
1010 }
1011
1012 //_____________________________________________________________________________
1013 void AliStack::FinishRun()
1014 {
1015 // Clean TreeK information
1016 }
1017 //_____________________________________________________________________________
1018
1019 Bool_t AliStack::GetEvent()
1020 {
1021 //
1022 // Get new event from TreeK
1023
1024     // Reset/Create the particle stack
1025     fTreeK = 0x0;
1026
1027     if (TreeK() == 0x0) //forces connecting
1028      {
1029       AliError("cannot find Kine Tree for current event");
1030       return kFALSE;
1031      }
1032       
1033     Int_t size = (Int_t)TreeK()->GetEntries();
1034     ResetArrays(size);
1035     return kTRUE;
1036 }
1037 //_____________________________________________________________________________
1038
1039 void AliStack::SetEventFolderName(const char* foldname)
1040 {
1041  //Sets event folder name
1042  fEventFolderName = foldname;
1043 }