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