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