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