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