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