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