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 /* $Id$ */
  17
  18 ///////////////////////////////////////////////////////////////////////////////
  19 //                                                                           //
  20 //  Particles stack class
  21 //                                                                           //
  22 ///////////////////////////////////////////////////////////////////////////////
  23
  24
  25 #include <TFile.h>
  26 #include <TFolder.h>
  27 #include <TObjArray.h>
  28 #include <TParticle.h>
  29 #include <TParticlePDG.h>
  30 #include <TTree.h>
  31
  32 #include "AliHit.h"
  33 #include "AliModule.h"
  34 #include "AliRun.h"
  35 #include "AliRunLoader.h"
  36 #include "AliStack.h"
  37
  38 ClassImp(AliStack)
  39
  40 //_______________________________________________________________________
  41 AliStack::AliStack():
  42   fParticles(0),
  43   fParticleMap(0),
  44   fParticleFileMap(0),
  45   fParticleBuffer(0),
  46   fTreeK(0),
  47   fNtrack(0),
  48   fNprimary(0),
  49   fCurrent(-1),
  50   fCurrentPrimary(-1),
  51   fHgwmk(0),
  52   fLoadPoint(0),
  53   fEventFolderName(AliConfig::fgkDefaultEventFolderName)
  54 {
  55   //
  56   // Default constructor
  57   //
  58 }
  59
  60 //_______________________________________________________________________
  61 AliStack::AliStack(Int_t size, const char* evfoldname):
  62   fParticles(new TClonesArray("TParticle",1000)),
  63   fParticleMap(new TObjArray(size)),
  64   fParticleFileMap(0),
  65   fParticleBuffer(0),
  66   fTreeK(0),
  67   fNtrack(0),
  68   fNprimary(0),
  69   fCurrent(-1),
  70   fCurrentPrimary(-1),
  71   fHgwmk(0),
  72   fLoadPoint(0),
  73   fEventFolderName(evfoldname)
  74 {
  75   //
  76   //  Constructor
  77   //
  78 }
  79
  80 //_______________________________________________________________________
  81 AliStack::AliStack(const AliStack& st):
  82   TVirtualMCStack(st),
  83   fParticles(0),
  84   fParticleMap(0),
  85   fParticleFileMap(0),
  86   fParticleBuffer(0),
  87   fTreeK(0),
  88   fNtrack(0),
  89   fNprimary(0),
  90   fCurrent(-1),
  91   fCurrentPrimary(-1),
  92   fHgwmk(0),
  93   fLoadPoint(0)
  94 {
  95   //
  96   // Copy constructor
  97   //
  98   st.Copy(*this);
  99 }
 100
 101 //_______________________________________________________________________
 102 void AliStack::Copy(AliStack&) const
 103 {
 104   Fatal("Copy","Not implemented!\n");
 105 }
 106
 107 //_______________________________________________________________________
 108 AliStack::~AliStack()
 109 {
 110   //
 111   // Destructor
 112   //
 113
 114   if (fParticles) {
 115     fParticles->Delete();
 116     delete fParticles;
 117   }
 118   delete fParticleMap;
 119   //PH???  delete fTreeK;
 120 }
 121
 122 //
 123 // public methods
 124 //
 125
 126 //_____________________________________________________________________________
 127 void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
 128                         Float_t *vpos, Float_t *polar, Float_t tof,
 129                         TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
 130 {
 131   //
 132   // Load a track on the stack
 133   //
 134   // done     0 if the track has to be transported
 135   //          1 if not
 136   // parent   identifier of the parent track. -1 for a primary
 137   // pdg    particle code
 138   // pmom     momentum GeV/c
 139   // vpos     position
 140   // polar    polarisation
 141   // tof      time of flight in seconds
 142   // mecha    production mechanism
 143   // ntr      on output the number of the track stored
 144   //
 145
 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     TParticlePDG* pmc =  TDatabasePDG::Instance()->GetParticle(pdg);
 156     if (pmc) {
 157         Float_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 158         Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+
 159                               pmom[1]*pmom[1]+pmom[2]*pmom[2]);
 160
 161 //    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",
 162 //         mass,e,fNtrack,pdg,parent,done,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS);
 163
 164
 165         PushTrack(done, parent, pdg, pmom[0], pmom[1], pmom[2], e,
 166                  vpos[0], vpos[1], vpos[2], tof, polar[0], polar[1], polar[2],
 167                  mech, ntr, weight, is);
 168     } else {
 169         Warning("PushTrack", "Particle type %d not defined in PDG Database !\n", pdg);
 170         Warning("PushTrack", "Particle skipped !\n");
 171     }
 172 }
 173
 174 //_____________________________________________________________________________
 175 void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg,
 176                       Double_t px, Double_t py, Double_t pz, Double_t e,
 177                       Double_t vx, Double_t vy, Double_t vz, Double_t tof,
 178                       Double_t polx, Double_t poly, Double_t polz,
 179                       TMCProcess mech, Int_t &ntr, Double_t weight, Int_t is)
 180 {
 181   //
 182   // Load a track on the stack
 183   //
 184   // done        0 if the track has to be transported
 185   //             1 if not
 186   // parent      identifier of the parent track. -1 for a primary
 187   // pdg         particle code
 188   // kS          generation status code
 189   // px, py, pz  momentum GeV/c
 190   // vx, vy, vz  position
 191   // polar       polarisation
 192   // tof         time of flight in seconds
 193   // mech        production mechanism
 194   // ntr         on output the number of the track stored
 195   //
 196   // New method interface:
 197   // arguments were changed to be in correspondence with TParticle
 198   // constructor.
 199   // Note: the energy is not calculated from the static mass but
 200   // it is passed by argument e.
 201
 202
 203   const Int_t kFirstDaughter=-1;
 204   const Int_t kLastDaughter=-1;
 205
 206   TClonesArray &particles = *fParticles;
 207   TParticle* particle
 208     = new(particles[fLoadPoint++])
 209       TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter,
 210                 px, py, pz, e, vx, vy, vz, tof);
 211
 212   particle->SetPolarisation(polx, poly, polz);
 213   particle->SetWeight(weight);
 214   particle->SetUniqueID(mech);
 215
 216   if(!done) particle->SetBit(kDoneBit);
 217
 218   //  Declare that the daughter information is valid
 219   particle->SetBit(kDaughtersBit);
 220   //  Add the particle to the stack
 221   fParticleMap->AddAtAndExpand(particle, fNtrack);//CHECK!!
 222
 223   if(parent>=0) {
 224     particle = dynamic_cast<TParticle*>(fParticleMap->At(parent));
 225     if (particle) {
 226       particle->SetLastDaughter(fNtrack);
 227       if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack);
 228     }
 229     else {
 230       printf("Error in AliStack::PushTrack: Parent %d does not exist\n",parent);
 231     }
 232   }
 233   else {
 234     //
 235     // This is a primary track. Set high water mark for this event
 236     fHgwmk = fNtrack;
 237     //
 238     // Set also number if primary tracks
 239     fNprimary = fHgwmk+1;
 240     fCurrentPrimary++;
 241   }
 242   ntr = fNtrack++;
 243 }
 244
 245 //_____________________________________________________________________________
 246 TParticle*  AliStack::PopNextTrack(Int_t& itrack)
 247 {
 248   //
 249   // Returns next track from stack of particles
 250   //
 251
 252
 253   TParticle* track = GetNextParticle();
 254
 255   if (track) {
 256     itrack = fCurrent;
 257     track->SetBit(kDoneBit);
 258   }
 259   else
 260     itrack = -1;
 261
 262   fCurrentTrack = track;
 263   return track;
 264 }
 265
 266 //_____________________________________________________________________________
 267 TParticle*  AliStack::PopPrimaryForTracking(Int_t i)
 268 {
 269   //
 270   // Returns i-th primary particle if it is flagged to be tracked,
 271   // 0 otherwise
 272   //
 273
 274   TParticle* particle = Particle(i);
 275
 276   if (!particle->TestBit(kDoneBit))
 277     return particle;
 278   else
 279     return 0;
 280 }
 281
 282 //_____________________________________________________________________________
 283 void AliStack::PurifyKine()
 284 {
 285   //
 286   // Compress kinematic tree keeping only flagged particles
 287   // and renaming the particle id's in all the hits
 288   //
 289
 290   TObjArray &particles = *fParticleMap;
 291   int nkeep=fHgwmk+1, parent, i;
 292   TParticle *part, *father;
 293   TArrayI map(particles.GetLast()+1);
 294
 295   // Save in Header total number of tracks before compression
 296
 297   // If no tracks generated return now
 298   if(fHgwmk+1 == fNtrack) return;
 299
 300   // First pass, invalid Daughter information
 301   for(i=0; i<fNtrack; i++) {
 302     // Preset map, to be removed later
 303     if(i<=fHgwmk) map[i]=i ;
 304     else {
 305       map[i] = -99;
 306       if((part=dynamic_cast<TParticle*>(particles.At(i)))) {
 307 //
 308 //        Check of this track should be kept for physics reasons
 309           if (KeepPhysics(part)) KeepTrack(i);
 310 //
 311           part->ResetBit(kDaughtersBit);
 312           part->SetFirstDaughter(-1);
 313           part->SetLastDaughter(-1);
 314       }
 315     }
 316   }
 317   // Invalid daughter information for the parent of the first particle
 318   // generated. This may or may not be the current primary according to
 319   // whether decays have been recorded among the primaries
 320   part = dynamic_cast<TParticle*>(particles.At(fHgwmk+1));
 321   particles.At(part->GetFirstMother())->ResetBit(kDaughtersBit);
 322   // Second pass, build map between old and new numbering
 323   for(i=fHgwmk+1; i<fNtrack; i++) {
 324     if(particles.At(i)->TestBit(kKeepBit)) {
 325
 326       // This particle has to be kept
 327       map[i]=nkeep;
 328       // If old and new are different, have to move the pointer
 329       if(i!=nkeep) particles[nkeep]=particles.At(i);
 330       part = dynamic_cast<TParticle*>(particles.At(nkeep));
 331
 332       // as the parent is always *before*, it must be already
 333       // in place. This is what we are checking anyway!
 334       if((parent=part->GetFirstMother())>fHgwmk)
 335         if(map[parent]==-99) Fatal("PurifyKine","map[%d] = -99!\n",parent);
 336         else part->SetFirstMother(map[parent]);
 337
 338       nkeep++;
 339     }
 340   }
 341
 342   // Fix daughters information
 343   for (i=fHgwmk+1; i<nkeep; i++) {
 344     part = dynamic_cast<TParticle*>(particles.At(i));
 345     parent = part->GetFirstMother();
 346     if(parent>=0) {
 347       father = dynamic_cast<TParticle*>(particles.At(parent));
 348       if(father->TestBit(kDaughtersBit)) {
 349
 350         if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i);
 351         if(i>father->GetLastDaughter())  father->SetLastDaughter(i);
 352       } else {
 353         // Initialise daughters info for first pass
 354         father->SetFirstDaughter(i);
 355         father->SetLastDaughter(i);
 356         father->SetBit(kDaughtersBit);
 357       }
 358     }
 359   }
 360
 361   // Now loop on all registered hit lists
 362   TList* hitLists = gAlice->GetHitLists();
 363   TIter next(hitLists);
 364   TCollection *hitList;
 365   while((hitList = dynamic_cast<TCollection*>(next()))) {
 366     TIter nexthit(hitList);
 367     AliHit *hit;
 368     while((hit = dynamic_cast<AliHit*>(nexthit()))) {
 369       hit->SetTrack(map[hit->GetTrack()]);
 370     }
 371   }
 372
 373   //
 374   // This for detectors which have a special mapping mechanism
 375   // for hits, such as TPC and TRD
 376   //
 377
 378    TObjArray* modules = gAlice->Modules();
 379    TIter nextmod(modules);
 380    AliModule *detector;
 381    while((detector = dynamic_cast<AliModule*>(nextmod()))) {
 382      detector->RemapTrackHitIDs(map.GetArray());
 383      detector->RemapTrackReferencesIDs(map.GetArray());
 384    }
 385    //
 386    gAlice->RemapTrackReferencesIDs(map.GetArray());
 387
 388    // Now the output bit, from fHgwmk to nkeep we write everything and we erase
 389    if(nkeep>fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5));
 390
 391    for (i=fHgwmk+1; i<nkeep; ++i) {
 392      fParticleBuffer = dynamic_cast<TParticle*>(particles.At(i));
 393      fParticleFileMap[i]=static_cast<Int_t>(TreeK()->GetEntries());
 394      TreeK()->Fill();
 395      particles[i]=fParticleBuffer=0;
 396     }
 397
 398    for (i=nkeep; i<fNtrack; ++i) particles[i]=0;
 399
 400    Int_t toshrink = fNtrack-fHgwmk-1;
 401    fLoadPoint-=toshrink;
 402    for(i=fLoadPoint; i<fLoadPoint+toshrink; ++i) fParticles->RemoveAt(i);
 403
 404    fNtrack=nkeep;
 405    fHgwmk=nkeep-1;
 406    //   delete [] map;
 407 }
 408
 409 Bool_t AliStack::KeepPhysics(TParticle* part)
 410 {
 411     //
 412     // Some particles have to kept on the stack for reasons motivated
 413     // by physics analysis. Decision is put here.
 414     //
 415     Bool_t keep = kFALSE;
 416     //
 417     // Keep first-generation daughter from primaries with heavy flavor
 418     //
 419     Int_t parent = part->GetFirstMother();
 420     if (parent >= 0 && parent <= fHgwmk) {
 421         TParticle* father = dynamic_cast<TParticle*>(Particles()->At(parent));
 422         Int_t kf = father->GetPdgCode();
 423         kf = TMath::Abs(kf);
 424         Int_t kfl = kf;
 425         // meson ?
 426         if  (kfl > 10) kfl/=100;
 427         // baryon
 428         if (kfl > 10) kfl/=10;
 429         if (kfl > 10) kfl/=10;
 430         if (kfl >= 4) {
 431             keep = kTRUE;
 432         }
 433     }
 434     return keep;
 435 }
 436
 437 //_____________________________________________________________________________
 438 void AliStack::FinishEvent()
 439 {
 440 //
 441 // Write out the kinematics that was not yet filled
 442 //
 443
 444 // Update event header
 445
 446
 447   if (!TreeK()) {
 448 //    Fatal("FinishEvent", "No kinematics tree is defined.");
 449 //    Don't panic this is a probably a lego run
 450       return;
 451   }
 452
 453   CleanParents();
 454   if(TreeK()->GetEntries() ==0) {
 455     // set the fParticleFileMap size for the first time
 456     fParticleFileMap.Set(fHgwmk+1);
 457   }
 458
 459   Bool_t allFilled = kFALSE;
 460   TObject *part;
 461   for(Int_t i=0; i<fHgwmk+1; ++i)
 462     if((part=fParticleMap->At(i))) {
 463       fParticleBuffer = dynamic_cast<TParticle*>(part);
 464       fParticleFileMap[i]= static_cast<Int_t>(TreeK()->GetEntries());
 465       TreeK()->Fill();
 466       //PH      (*fParticleMap)[i]=fParticleBuffer=0;
 467       fParticleBuffer=0;
 468       fParticleMap->AddAt(0,i);
 469
 470       // When all primaries were filled no particle!=0
 471       // should be left => to be removed later.
 472       if (allFilled) printf("Why != 0 part # %d?\n",i);
 473     }
 474     else
 475      {
 476       // // printf("Why = 0 part # %d?\n",i); => We know.
 477       // break;
 478       // we don't break now in order to be sure there is no
 479       // particle !=0 left.
 480       // To be removed later and replaced with break.
 481        if(!allFilled) allFilled = kTRUE;
 482      }
 483 }
 484 //_____________________________________________________________________________
 485
 486 void AliStack::FlagTrack(Int_t track)
 487 {
 488   //
 489   // Flags a track and all its family tree to be kept
 490   //
 491
 492   TParticle *particle;
 493
 494   Int_t curr=track;
 495   while(1) {
 496     particle=dynamic_cast<TParticle*>(fParticleMap->At(curr));
 497
 498     // If the particle is flagged the three from here upward is saved already
 499     if(particle->TestBit(kKeepBit)) return;
 500
 501     // Save this particle
 502     particle->SetBit(kKeepBit);
 503
 504     // Move to father if any
 505     if((curr=particle->GetFirstMother())==-1) return;
 506   }
 507 }
 508
 509 //_____________________________________________________________________________
 510 void AliStack::KeepTrack(Int_t track)
 511 {
 512   //
 513   // Flags a track to be kept
 514   //
 515
 516   fParticleMap->At(track)->SetBit(kKeepBit);
 517 }
 518
 519 //_____________________________________________________________________________
 520 void  AliStack::Reset(Int_t size)
 521 {
 522   //
 523   // Resets stack
 524   //
 525
 526   fNtrack=0;
 527   fNprimary=0;
 528   fHgwmk=0;
 529   fLoadPoint=0;
 530   fCurrent = -1;
 531   fTreeK = 0x0;
 532   ResetArrays(size);
 533 }
 534
 535 //_____________________________________________________________________________
 536 void  AliStack::ResetArrays(Int_t size)
 537 {
 538   //
 539   // Resets stack arrays
 540   //
 541
 542   if (fParticles)
 543     fParticles->Clear();
 544   else
 545     fParticles = new TClonesArray("TParticle",1000);
 546   if (fParticleMap) {
 547     fParticleMap->Clear();
 548     if (size>0) fParticleMap->Expand(size);}
 549   else
 550     fParticleMap = new TObjArray(size);
 551 }
 552
 553 //_____________________________________________________________________________
 554 void AliStack::SetHighWaterMark(Int_t)
 555 {
 556   //
 557   // Set high water mark for last track in event
 558   //
 559
 560   fHgwmk = fNtrack-1;
 561   fCurrentPrimary=fHgwmk;
 562
 563   // Set also number of primary tracks
 564   fNprimary = fHgwmk+1;
 565   fNtrack   = fHgwmk+1;
 566 }
 567
 568 //_____________________________________________________________________________
 569 TParticle* AliStack::Particle(Int_t i)
 570 {
 571   //
 572   // Return particle with specified ID
 573
 574   //PH  if(!(*fParticleMap)[i]) {
 575   if(!fParticleMap->At(i)) {
 576     Int_t nentries = fParticles->GetEntriesFast();
 577     // algorithmic way of getting entry index
 578     // (primary particles are filled after secondaries)
 579     Int_t entry = TreeKEntry(i);
 580     // check whether algorithmic way and
 581     // and the fParticleFileMap[i] give the same;
 582     // give the fatal error if not
 583     if (entry != fParticleFileMap[i]) {
 584       Fatal("Particle",
 585         "!! The algorithmic way and map are different: !!\n entry: %d map: %d",
 586         entry, fParticleFileMap[i]);
 587     }
 588
 589     TreeK()->GetEntry(entry);
 590     new ((*fParticles)[nentries]) TParticle(*fParticleBuffer);
 591     fParticleMap->AddAt((*fParticles)[nentries],i);
 592   }
 593   //PH  return dynamic_cast<TParticle *>((*fParticleMap)[i]);
 594   return dynamic_cast<TParticle*>(fParticleMap->At(i));
 595 }
 596
 597 //_____________________________________________________________________________
 598 TParticle* AliStack::ParticleFromTreeK(Int_t id) const
 599 {
 600 //
 601 // return pointer to TParticle with label id
 602 //
 603   Int_t entry;
 604   if ((entry = TreeKEntry(id)) < 0) return 0;
 605   if (fTreeK->GetEntry(entry)<=0) return 0;
 606   return fParticleBuffer;
 607 }
 608
 609 //_____________________________________________________________________________
 610 Int_t AliStack::TreeKEntry(Int_t id) const
 611 {
 612 //
 613 // return entry number in the TreeK for particle with label id
 614 // return negative number if label>fNtrack
 615 //
 616   Int_t entry;
 617   if (id<fNprimary)
 618     entry = id+fNtrack-fNprimary;
 619   else
 620     entry = id-fNprimary;
 621   return entry;
 622 }
 623
 624 //_____________________________________________________________________________
 625 Int_t AliStack::GetCurrentParentTrackNumber() const
 626 {
 627   //
 628   // Return number of the parent of the current track
 629   //
 630
 631   TParticle* current = (TParticle*)fParticleMap->At(fCurrent);
 632
 633   if (current)
 634     return current->GetFirstMother();
 635   else {
 636     Warning("GetCurrentParentTrackNumber", "Current track not found in the stack");
 637     return -1;
 638   }
 639 }
 640
 641 //_____________________________________________________________________________
 642 Int_t AliStack::GetPrimary(Int_t id)
 643 {
 644   //
 645   // Return number of primary that has generated track
 646   //
 647
 648   int current, parent;
 649   //
 650   parent=id;
 651   while (1) {
 652     current=parent;
 653     parent=Particle(current)->GetFirstMother();
 654     if(parent<0) return current;
 655   }
 656 }
 657
 658 //_____________________________________________________________________________
 659 void AliStack::DumpPart (Int_t i) const
 660 {
 661   //
 662   // Dumps particle i in the stack
 663   //
 664
 665   //PH  dynamic_cast<TParticle*>((*fParticleMap)[i])->Print();
 666   dynamic_cast<TParticle*>(fParticleMap->At(i))->Print();
 667 }
 668
 669 //_____________________________________________________________________________
 670 void AliStack::DumpPStack ()
 671 {
 672   //
 673   // Dumps the particle stack
 674   //
 675
 676   Int_t i;
 677
 678   printf("\n\n=======================================================================\n");
 679   for (i=0;i<fNtrack;i++)
 680     {
 681       TParticle* particle = Particle(i);
 682       if (particle) {
 683         printf("-> %d ",i); particle->Print();
 684         printf("--------------------------------------------------------------\n");
 685       }
 686       else
 687         Warning("DumpPStack", "No particle with id %d.", i);
 688     }
 689
 690   printf("\n=======================================================================\n\n");
 691
 692   // print  particle file map
 693   printf("\nParticle file map: \n");
 694   for (i=0; i<fNtrack; i++)
 695       printf("   %d th entry: %d \n",i,fParticleFileMap[i]);
 696 }
 697
 698
 699 //_____________________________________________________________________________
 700 void AliStack::DumpLoadedStack() const
 701 {
 702   //
 703   // Dumps the particle in the stack
 704   // that are loaded in memory.
 705   //
 706
 707   TObjArray &particles = *fParticleMap;
 708   printf(
 709          "\n\n=======================================================================\n");
 710   for (Int_t i=0;i<fNtrack;i++)
 711     {
 712       TParticle* particle = dynamic_cast<TParticle*>(particles[i]);
 713       if (particle) {
 714         printf("-> %d ",i); particle->Print();
 715         printf("--------------------------------------------------------------\n");
 716       }
 717       else {
 718         printf("-> %d  Particle not loaded.\n",i);
 719         printf("--------------------------------------------------------------\n");
 720       }
 721     }
 722   printf(
 723          "\n=======================================================================\n\n");
 724 }
 725
 726 //
 727 // protected methods
 728 //
 729
 730 //_____________________________________________________________________________
 731 void AliStack::CleanParents()
 732 {
 733   //
 734   // Clean particles stack
 735   // Set parent/daughter relations
 736   //
 737
 738   TObjArray &particles = *fParticleMap;
 739   TParticle *part;
 740   int i;
 741   for(i=0; i<fHgwmk+1; i++) {
 742     part = dynamic_cast<TParticle*>(particles.At(i));
 743     if(part) if(!part->TestBit(kDaughtersBit)) {
 744       part->SetFirstDaughter(-1);
 745       part->SetLastDaughter(-1);
 746     }
 747   }
 748 }
 749
 750 //_____________________________________________________________________________
 751 TParticle* AliStack::GetNextParticle()
 752 {
 753   //
 754   // Return next particle from stack of particles
 755   //
 756
 757   TParticle* particle = 0;
 758
 759   // search secondaries
 760   //for(Int_t i=fNtrack-1; i>=0; i--) {
 761   for(Int_t i=fNtrack-1; i>fHgwmk; i--) {
 762       particle = dynamic_cast<TParticle*>(fParticleMap->At(i));
 763       if ((particle) && (!particle->TestBit(kDoneBit))) {
 764           fCurrent=i;
 765           return particle;
 766       }
 767   }
 768
 769   // take next primary if all secondaries were done
 770   while (fCurrentPrimary>=0) {
 771       fCurrent = fCurrentPrimary;
 772       particle = dynamic_cast<TParticle*>(fParticleMap->At(fCurrentPrimary--));
 773       if ((particle) && (!particle->TestBit(kDoneBit))) {
 774           return particle;
 775       }
 776   }
 777
 778   // nothing to be tracked
 779   fCurrent = -1;
 780   return particle;
 781 }
 782 //__________________________________________________________________________________________
 783
 784 TTree* AliStack::TreeK()
 785 {
 786 //returns TreeK
 787   if (fTreeK)
 788    {
 789      return fTreeK;
 790    }
 791   else
 792    {
 793      AliRunLoader *rl = AliRunLoader::GetRunLoader(fEventFolderName);
 794      if (rl == 0x0)
 795       {
 796         Fatal("TreeK","Can not get RunLoader from event folder named %s",fEventFolderName.Data());
 797         return 0x0;//pro forma
 798       }
 799     fTreeK = rl->TreeK();
 800     if ( fTreeK == 0x0)
 801      {
 802       Error("TreeK","Can not get TreeK from RL. Ev. Folder is %s",fEventFolderName.Data());
 803      }
 804     else
 805      {
 806       ConnectTree();
 807      }
 808    }
 809   return fTreeK;//never reached
 810 }
 811 //__________________________________________________________________________________________
 812
 813 void AliStack::ConnectTree()
 814 {
 815 //
 816 //  Creates branch for writing particles
 817 //
 818   if (AliLoader::fgDebug) Info("ConnectTree","Connecting TreeK");
 819   if (fTreeK == 0x0)
 820    {
 821     if (TreeK() == 0x0)
 822      {
 823       Fatal("ConnectTree","Parameter is NULL");//we don't like such a jokes
 824       return;
 825      }
 826     return;//in this case TreeK() calls back this method (ConnectTree)
 827            //tree after setting fTreeK, the rest was already executed
 828            //it is safe to return now
 829    }
 830
 831  //  Create a branch for particles
 832
 833   if (AliLoader::fgDebug)
 834    Info("ConnectTree","Tree name is %s",fTreeK->GetName());
 835
 836   if (fTreeK->GetDirectory())
 837    {
 838      if (AliLoader::fgDebug)
 839       Info("ConnectTree","and dir is %s",fTreeK->GetDirectory()->GetName());
 840    }
 841   else
 842     Warning("ConnectTree","DIR IS NOT SET !!!");
 843
 844   TBranch *branch=fTreeK->GetBranch(AliRunLoader::fgkKineBranchName);
 845   if(branch == 0x0)
 846    {
 847     branch = fTreeK->Branch(AliRunLoader::fgkKineBranchName, "TParticle", &fParticleBuffer, 4000);
 848     if (AliLoader::fgDebug) Info("ConnectTree","Creating Branch in Tree");
 849    }
 850   else
 851    {
 852     if (AliLoader::fgDebug) Info("ConnectTree","Branch Found in Tree");
 853     branch->SetAddress(&fParticleBuffer);
 854    }
 855   if (branch->GetDirectory())
 856    {
 857     if (AliLoader::fgDebug)
 858       Info("ConnectTree","Branch Dir Name is %s",branch->GetDirectory()->GetName());
 859    }
 860   else
 861     Warning("ConnectTree","Branch Dir is NOT SET");
 862 }
 863 //__________________________________________________________________________________________
 864
 865
 866 void AliStack::BeginEvent()
 867 {
 868 // start a new event
 869  Reset();
 870 }
 871
 872 //_____________________________________________________________________________
 873 void AliStack::FinishRun()
 874 {
 875 // Clean TreeK information
 876 }
 877 //_____________________________________________________________________________
 878
 879 Bool_t AliStack::GetEvent()
 880 {
 881 //
 882 // Get new event from TreeK
 883
 884     // Reset/Create the particle stack
 885     fTreeK = 0x0;
 886
 887     if (TreeK() == 0x0) //forces connecting
 888      {
 889       Error("GetEvent","cannot find Kine Tree for current event\n");
 890       return kFALSE;
 891      }
 892
 893     Int_t size = (Int_t)TreeK()->GetEntries();
 894     ResetArrays(size);
 895     return kTRUE;
 896 }
 897 //_____________________________________________________________________________
 898
 899 void AliStack::SetEventFolderName(const char* foldname)
 900 {
 901  //Sets event folder name
 902  fEventFolderName = foldname;
 903 }