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