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