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