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