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