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