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