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