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