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 //  Implements the TMCVirtualStack of the Virtual Monte Carlo                //
  22 //  Holds the particles transported during simulation                        //
  23 //  Is used to compare results of reconstruction with simulation             //
  24 //  Author A.Morsch                                                          //
  25 //                                                                           //
  26 ///////////////////////////////////////////////////////////////////////////////
  27
  28
  29 #include <TClonesArray.h>
  30 #include <TObjArray.h>
  31 #include <TPDGCode.h>
  32 #include <TParticle.h>
  33 #include <TParticlePDG.h>
  34 #include <TTree.h>
  35
  36 #include "AliLog.h"
  37 #include "AliHit.h"
  38 #include "AliModule.h"
  39 #include "AliRun.h"
  40 #include "AliMC.h"
  41 #include "AliRunLoader.h"
  42 #include "AliStack.h"
  43
  44 ClassImp(AliStack)
  45
  46 //_______________________________________________________________________
  47 AliStack::AliStack():
  48   fParticles(0),
  49   fParticleMap(0),
  50   fParticleFileMap(0),
  51   fParticleBuffer(0),
  52   fCurrentTrack(0),
  53   fTreeK(0),
  54   fNtrack(0),
  55   fNprimary(0),
  56   fCurrent(-1),
  57   fCurrentPrimary(-1),
  58   fHgwmk(0),
  59   fLoadPoint(0),
  60   fEventFolderName(AliConfig::GetDefaultEventFolderName())
  61 {
  62   //
  63   // Default constructor
  64   //
  65 }
  66
  67 //_______________________________________________________________________
  68 AliStack::AliStack(Int_t size, const char* evfoldname):
  69   fParticles(new TClonesArray("TParticle",1000)),
  70   fParticleMap(new TObjArray(size)),
  71   fParticleFileMap(0),
  72   fParticleBuffer(0),
  73   fCurrentTrack(0),
  74   fTreeK(0),
  75   fNtrack(0),
  76   fNprimary(0),
  77   fCurrent(-1),
  78   fCurrentPrimary(-1),
  79   fHgwmk(0),
  80   fLoadPoint(0),
  81   fEventFolderName(evfoldname)
  82 {
  83   //
  84   //  Constructor
  85   //
  86 }
  87
  88 //_______________________________________________________________________
  89 AliStack::AliStack(const AliStack& st):
  90   TVirtualMCStack(st),
  91   fParticles(0),
  92   fParticleMap(0),
  93   fParticleFileMap(0),
  94   fParticleBuffer(0),
  95   fCurrentTrack(0),
  96   fTreeK(0),
  97   fNtrack(0),
  98   fNprimary(0),
  99   fCurrent(-1),
 100   fCurrentPrimary(-1),
 101   fHgwmk(0),
 102   fLoadPoint(0),
 103   fEventFolderName()
 104 {
 105   //
 106   // Copy constructor
 107   //
 108   st.Copy(*this);
 109 }
 110
 111 //_______________________________________________________________________
 112 void AliStack::Copy(TObject&) const
 113 {
 114   AliFatal("Not implemented!");
 115 }
 116
 117 //_______________________________________________________________________
 118 AliStack::~AliStack()
 119 {
 120   //
 121   // Destructor
 122   //
 123
 124   if (fParticles) {
 125     fParticles->Delete();
 126     delete fParticles;
 127   }
 128   delete fParticleMap;
 129   //PH???  delete fTreeK;
 130 }
 131
 132 //
 133 // public methods
 134 //
 135
 136 //_____________________________________________________________________________
 137 void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
 138                         Float_t *vpos, Float_t *polar, Float_t tof,
 139                         TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
 140 {
 141   //
 142   // Load a track on the stack
 143   //
 144   // done     0 if the track has to be transported
 145   //          1 if not
 146   // parent   identifier of the parent track. -1 for a primary
 147   // pdg    particle code
 148   // pmom     momentum GeV/c
 149   // vpos     position
 150   // polar    polarisation
 151   // tof      time of flight in seconds
 152   // mecha    production mechanism
 153   // ntr      on output the number of the track stored
 154   //
 155
 156   //  const Float_t tlife=0;
 157
 158   //
 159   // Here we get the static mass
 160   // For MC is ok, but a more sophisticated method could be necessary
 161   // if the calculated mass is required
 162   // also, this method is potentially dangerous if the mass
 163   // used in the MC is not the same of the PDG database
 164   //
 165     TParticlePDG* pmc =  TDatabasePDG::Instance()->GetParticle(pdg);
 166     if (pmc) {
 167         Float_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 168         Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+
 169                               pmom[1]*pmom[1]+pmom[2]*pmom[2]);
 170
 171 //    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",
 172 //         mass,e,fNtrack,pdg,parent,done,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS);
 173
 174
 175         PushTrack(done, parent, pdg, pmom[0], pmom[1], pmom[2], e,
 176                  vpos[0], vpos[1], vpos[2], tof, polar[0], polar[1], polar[2],
 177                  mech, ntr, weight, is);
 178     } else {
 179         AliWarning(Form("Particle type %d not defined in PDG Database !", pdg));
 180         AliWarning("Particle skipped !");
 181     }
 182 }
 183
 184 //_____________________________________________________________________________
 185 void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg,
 186                       Double_t px, Double_t py, Double_t pz, Double_t e,
 187                       Double_t vx, Double_t vy, Double_t vz, Double_t tof,
 188                       Double_t polx, Double_t poly, Double_t polz,
 189                       TMCProcess mech, Int_t &ntr, Double_t weight, Int_t is)
 190 {
 191   //
 192   // Load a track on the stack
 193   //
 194   // done        0 if the track has to be transported
 195   //             1 if not
 196   // parent      identifier of the parent track. -1 for a primary
 197   // pdg         particle code
 198   // kS          generation status code
 199   // px, py, pz  momentum GeV/c
 200   // vx, vy, vz  position
 201   // polar       polarisation
 202   // tof         time of flight in seconds
 203   // mech        production mechanism
 204   // ntr         on output the number of the track stored
 205   //
 206   // New method interface:
 207   // arguments were changed to be in correspondence with TParticle
 208   // constructor.
 209   // Note: the energy is not calculated from the static mass but
 210   // it is passed by argument e.
 211
 212
 213   const Int_t kFirstDaughter=-1;
 214   const Int_t kLastDaughter=-1;
 215
 216
 217   TClonesArray &particles = *fParticles;
 218   TParticle* particle
 219     = new(particles[fLoadPoint++])
 220       TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter,
 221                 px, py, pz, e, vx, vy, vz, tof);
 222
 223   particle->SetPolarisation(polx, poly, polz);
 224   particle->SetWeight(weight);
 225   particle->SetUniqueID(mech);
 226
 227   if(!done) particle->SetBit(kDoneBit);
 228
 229   //  Declare that the daughter information is valid
 230   particle->SetBit(kDaughtersBit);
 231   //  Add the particle to the stack
 232
 233
 234   fParticleMap->AddAtAndExpand(particle, fNtrack);//CHECK!!
 235
 236   if(parent>=0) {
 237     particle = dynamic_cast<TParticle*>(fParticleMap->At(parent));
 238     if (particle) {
 239       particle->SetLastDaughter(fNtrack);
 240       if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack);
 241     }
 242     else {
 243       AliError(Form("Parent %d does not exist",parent));
 244     }
 245   }
 246   else {
 247     //
 248     // This is a primary track. Set high water mark for this event
 249     fHgwmk = fNtrack;
 250     //
 251     // Set also number if primary tracks
 252     fNprimary = fHgwmk+1;
 253     fCurrentPrimary++;
 254   }
 255   ntr = fNtrack++;
 256 }
 257
 258 //_____________________________________________________________________________
 259 TParticle*  AliStack::PopNextTrack(Int_t& itrack)
 260 {
 261   //
 262   // Returns next track from stack of particles
 263   //
 264
 265
 266   TParticle* track = GetNextParticle();
 267
 268   if (track) {
 269     itrack = fCurrent;
 270     track->SetBit(kDoneBit);
 271   }
 272   else
 273     itrack = -1;
 274
 275   fCurrentTrack = track;
 276   return track;
 277 }
 278
 279 //_____________________________________________________________________________
 280 TParticle*  AliStack::PopPrimaryForTracking(Int_t i)
 281 {
 282   //
 283   // Returns i-th primary particle if it is flagged to be tracked,
 284   // 0 otherwise
 285   //
 286
 287   TParticle* particle = Particle(i);
 288
 289   if (!particle->TestBit(kDoneBit))
 290     return particle;
 291   else
 292     return 0;
 293 }
 294
 295 //_____________________________________________________________________________
 296 void AliStack::PurifyKine()
 297 {
 298   //
 299   // Compress kinematic tree keeping only flagged particles
 300   // and renaming the particle id's in all the hits
 301   //
 302
 303   TObjArray &particles = *fParticleMap;
 304   int nkeep=fHgwmk+1, parent, i;
 305   TParticle *part, *father;
 306   TArrayI map(particles.GetLast()+1);
 307
 308   // Save in Header total number of tracks before compression
 309   // If no tracks generated return now
 310   if(fHgwmk+1 == fNtrack) return;
 311
 312   // First pass, invalid Daughter information
 313    for(i=0; i<fNtrack; i++) {
 314     // Preset map, to be removed later
 315     if(i<=fHgwmk) map[i]=i ;
 316     else {
 317       map[i] = -99;
 318       if((part=dynamic_cast<TParticle*>(particles.At(i)))) {
 319 //
 320 //        Check of this track should be kept for physics reasons
 321           if (KeepPhysics(part)) KeepTrack(i);
 322 //
 323           part->ResetBit(kDaughtersBit);
 324           part->SetFirstDaughter(-1);
 325           part->SetLastDaughter(-1);
 326       }
 327     }
 328   }
 329   // Invalid daughter information for the parent of the first particle
 330   // generated. This may or may not be the current primary according to
 331   // whether decays have been recorded among the primaries
 332   part = dynamic_cast<TParticle*>(particles.At(fHgwmk+1));
 333   particles.At(part->GetFirstMother())->ResetBit(kDaughtersBit);
 334   // Second pass, build map between old and new numbering
 335   for(i=fHgwmk+1; i<fNtrack; i++) {
 336     if(particles.At(i)->TestBit(kKeepBit)) {
 337
 338       // This particle has to be kept
 339       map[i]=nkeep;
 340       // If old and new are different, have to move the pointer
 341       if(i!=nkeep) particles[nkeep]=particles.At(i);
 342       part = dynamic_cast<TParticle*>(particles.At(nkeep));
 343
 344       // as the parent is always *before*, it must be already
 345       // in place. This is what we are checking anyway!
 346       if((parent=part->GetFirstMother())>fHgwmk)
 347         if(map[parent]==-99) Fatal("PurifyKine","map[%d] = -99!\n",parent);
 348         else part->SetFirstMother(map[parent]);
 349
 350       nkeep++;
 351     }
 352   }
 353
 354   // Fix daughters information
 355   for (i=fHgwmk+1; i<nkeep; i++) {
 356     part = dynamic_cast<TParticle*>(particles.At(i));
 357     parent = part->GetFirstMother();
 358     if(parent>=0) {
 359       father = dynamic_cast<TParticle*>(particles.At(parent));
 360       if(father->TestBit(kDaughtersBit)) {
 361
 362         if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i);
 363         if(i>father->GetLastDaughter())  father->SetLastDaughter(i);
 364       } else {
 365         // Initialise daughters info for first pass
 366         father->SetFirstDaughter(i);
 367         father->SetLastDaughter(i);
 368         father->SetBit(kDaughtersBit);
 369       }
 370     }
 371   }
 372
 373
 374   // Now loop on all registered hit lists
 375   TList* hitLists = gAlice->GetMCApp()->GetHitLists();
 376   TIter next(hitLists);
 377   TCollection *hitList;
 378
 379   while((hitList = dynamic_cast<TCollection*>(next()))) {
 380     TIter nexthit(hitList);
 381     AliHit *hit;
 382
 383     while((hit = dynamic_cast<AliHit*>(nexthit()))) {
 384         hit->SetTrack(map[hit->GetTrack()]);
 385     }
 386   }
 387
 388   //
 389   // This for detectors which have a special mapping mechanism
 390   // for hits, such as TPC and TRD
 391   //
 392
 393    TObjArray* modules = gAlice->Modules();
 394    TIter nextmod(modules);
 395    AliModule *detector;
 396    while((detector = dynamic_cast<AliModule*>(nextmod()))) {
 397      detector->RemapTrackHitIDs(map.GetArray());
 398      detector->RemapTrackReferencesIDs(map.GetArray());
 399    }
 400    //
 401    gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray());
 402
 403    // Now the output bit, from fHgwmk to nkeep we write everything and we erase
 404    if(nkeep>fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5));
 405
 406    for (i=fHgwmk+1; i<nkeep; ++i) {
 407      fParticleBuffer = dynamic_cast<TParticle*>(particles.At(i));
 408      fParticleFileMap[i]=static_cast<Int_t>(TreeK()->GetEntries());
 409      TreeK()->Fill();
 410      particles[i]=fParticleBuffer=0;
 411     }
 412
 413    for (i=nkeep; i<fNtrack; ++i) particles[i]=0;
 414
 415    Int_t toshrink = fNtrack-fHgwmk-1;
 416    fLoadPoint-=toshrink;
 417
 418
 419    for(i=fLoadPoint; i<fLoadPoint+toshrink; ++i) fParticles->RemoveAt(i);
 420
 421    fNtrack=nkeep;
 422    fHgwmk=nkeep-1;
 423 }
 424
 425 void AliStack::ReorderKine()
 426 {
 427 //
 428 // In some transport code children might not come in a continuous sequence.
 429 // In this case the stack  has  to  be reordered in order to establish the
 430 // mother daughter relation using index ranges.
 431 //
 432   if(fHgwmk+1 == fNtrack) return;
 433
 434   //
 435   // Howmany secondaries have been produced ?
 436   Int_t nNew = fNtrack - fHgwmk - 1;
 437
 438   if (nNew > 0) {
 439       Int_t i, j;
 440       TObjArray &particles = *fParticleMap;
 441       TArrayI map1(nNew);
 442       //
 443       // Copy pointers to temporary array
 444       TParticle** tmp = new TParticle*[nNew];
 445
 446       for (i = 0; i < nNew; i++) {
 447           if (particles.At(fHgwmk + 1 + i)) {
 448               tmp[i] = (TParticle*) (particles.At(fHgwmk + 1 + i));
 449           } else {
 450               tmp[i] = 0x0;
 451           }
 452           map1[i] = -99;
 453       }
 454
 455
 456       //
 457       // Reset  LoadPoint
 458       //
 459       fLoadPoint = fHgwmk + 1;
 460       //
 461       // Re-Push particles into stack
 462       // The outer loop is over parents, the inner over children.
 463       // -1 refers to the primary particle
 464       //
 465       for (i = -1; i < nNew-1; i++) {
 466           Int_t ipa;
 467           TParticle* parP;
 468           if (i == -1) {
 469               ipa  = tmp[0]->GetFirstMother();
 470               parP =dynamic_cast<TParticle*>(particles.At(ipa));
 471           } else {
 472               ipa = (fHgwmk + 1 + i);
 473               // Skip deleted particles
 474               if (!tmp[i])                          continue;
 475               // Skip particles without children
 476               if (tmp[i]->GetFirstDaughter() == -1) continue;
 477               parP = tmp[i];
 478           }
 479           // Reset daughter information
 480
 481           Int_t idaumin = parP->GetFirstDaughter() - fHgwmk - 1;
 482           Int_t idaumax = parP->GetLastDaughter()  - fHgwmk - 1;
 483           parP->SetFirstDaughter(-1);
 484           parP->SetLastDaughter(-1);
 485           for (j = idaumin; j <= idaumax; j++) {
 486               // Skip deleted particles
 487               if (!tmp[j])        continue;
 488               // Skip particles already handled
 489               if (map1[j] != -99) continue;
 490               Int_t jpa = tmp[j]->GetFirstMother();
 491               // Check if daughter of current parent
 492               if (jpa == ipa) {
 493                   particles[fLoadPoint] = tmp[j];
 494                   // Re-establish daughter information
 495                   parP->SetLastDaughter(fLoadPoint);
 496                   if (parP->GetFirstDaughter() == -1) parP->SetFirstDaughter(fLoadPoint);
 497                   // Set Mother information
 498                   if (i != -1) {
 499                       tmp[j]->SetFirstMother(map1[i]);
 500                   }
 501                   // Build the map
 502                   map1[j] = fLoadPoint;
 503                   // Increase load point
 504                   fLoadPoint++;
 505               }
 506           } // children
 507       } // parents
 508
 509       delete[] tmp;
 510
 511       //
 512       // Build map for remapping of hits
 513       //
 514       TArrayI map(fNtrack);
 515       for (i = 0; i < fNtrack; i ++) {
 516           if (i <= fHgwmk) {
 517               map[i] = i;
 518           } else{
 519               map[i] = map1[i - fHgwmk -1];
 520           }
 521       }
 522
 523       // Now loop on all registered hit lists
 524
 525       TList* hitLists = gAlice->GetMCApp()->GetHitLists();
 526       TIter next(hitLists);
 527       TCollection *hitList;
 528
 529       while((hitList = dynamic_cast<TCollection*>(next()))) {
 530           TIter nexthit(hitList);
 531           AliHit *hit;
 532           while((hit = dynamic_cast<AliHit*>(nexthit()))) {
 533               hit->SetTrack(map[hit->GetTrack()]);
 534           }
 535       }
 536
 537   //
 538   // This for detectors which have a special mapping mechanism
 539   // for hits, such as TPC and TRD
 540   //
 541
 542       TObjArray* modules = gAlice->Modules();
 543       TIter nextmod(modules);
 544       AliModule *detector;
 545       while((detector = dynamic_cast<AliModule*>(nextmod()))) {
 546           detector->RemapTrackHitIDs(map.GetArray());
 547           detector->RemapTrackReferencesIDs(map.GetArray());
 548       }
 549       gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray());
 550   } // new particles poduced
 551 }
 552
 553 Bool_t AliStack::KeepPhysics(TParticle* part)
 554 {
 555     //
 556     // Some particles have to kept on the stack for reasons motivated
 557     // by physics analysis. Decision is put here.
 558     //
 559     Bool_t keep = kFALSE;
 560     //
 561     // Keep first-generation daughter from primaries with heavy flavor
 562     //
 563     Int_t parent = part->GetFirstMother();
 564     if (parent >= 0 && parent <= fHgwmk) {
 565         TParticle* father = dynamic_cast<TParticle*>(Particles()->At(parent));
 566         Int_t kf = father->GetPdgCode();
 567         kf = TMath::Abs(kf);
 568         Int_t kfl = kf;
 569         // meson ?
 570         if  (kfl > 10) kfl/=100;
 571         // baryon
 572         if (kfl > 10) kfl/=10;
 573         if (kfl > 10) kfl/=10;
 574         if (kfl >= 4) {
 575             keep = kTRUE;
 576         }
 577     }
 578     return keep;
 579 }
 580
 581 //_____________________________________________________________________________
 582 void AliStack::FinishEvent()
 583 {
 584 //
 585 // Write out the kinematics that was not yet filled
 586 //
 587
 588 // Update event header
 589
 590
 591   if (!TreeK()) {
 592 //    Fatal("FinishEvent", "No kinematics tree is defined.");
 593 //    Don't panic this is a probably a lego run
 594       return;
 595   }
 596
 597   CleanParents();
 598   if(TreeK()->GetEntries() ==0) {
 599     // set the fParticleFileMap size for the first time
 600     fParticleFileMap.Set(fHgwmk+1);
 601   }
 602
 603   Bool_t allFilled = kFALSE;
 604   TObject *part;
 605   for(Int_t i=0; i<fHgwmk+1; ++i)
 606     if((part=fParticleMap->At(i))) {
 607       fParticleBuffer = dynamic_cast<TParticle*>(part);
 608       fParticleFileMap[i]= static_cast<Int_t>(TreeK()->GetEntries());
 609       TreeK()->Fill();
 610       //PH      (*fParticleMap)[i]=fParticleBuffer=0;
 611       fParticleBuffer=0;
 612       fParticleMap->AddAt(0,i);
 613
 614       // When all primaries were filled no particle!=0
 615       // should be left => to be removed later.
 616       if (allFilled) AliWarning(Form("Why != 0 part # %d?\n",i));
 617     }
 618     else
 619      {
 620       // // printf("Why = 0 part # %d?\n",i); => We know.
 621       // break;
 622       // we don't break now in order to be sure there is no
 623       // particle !=0 left.
 624       // To be removed later and replaced with break.
 625        if(!allFilled) allFilled = kTRUE;
 626      }
 627 }
 628 //_____________________________________________________________________________
 629
 630 void AliStack::FlagTrack(Int_t track)
 631 {
 632   //
 633   // Flags a track and all its family tree to be kept
 634   //
 635
 636   TParticle *particle;
 637
 638   Int_t curr=track;
 639   while(1) {
 640     particle=dynamic_cast<TParticle*>(fParticleMap->At(curr));
 641
 642     // If the particle is flagged the three from here upward is saved already
 643     if(particle->TestBit(kKeepBit)) return;
 644
 645     // Save this particle
 646     particle->SetBit(kKeepBit);
 647
 648     // Move to father if any
 649     if((curr=particle->GetFirstMother())==-1) return;
 650   }
 651 }
 652
 653 //_____________________________________________________________________________
 654 void AliStack::KeepTrack(Int_t track)
 655 {
 656   //
 657   // Flags a track to be kept
 658   //
 659
 660   fParticleMap->At(track)->SetBit(kKeepBit);
 661 }
 662
 663 //_____________________________________________________________________________
 664 void  AliStack::Reset(Int_t size)
 665 {
 666   //
 667   // Resets stack
 668   //
 669
 670   fNtrack=0;
 671   fNprimary=0;
 672   fHgwmk=0;
 673   fLoadPoint=0;
 674   fCurrent = -1;
 675   fTreeK = 0x0;
 676   ResetArrays(size);
 677 }
 678
 679 //_____________________________________________________________________________
 680 void  AliStack::ResetArrays(Int_t size)
 681 {
 682   //
 683   // Resets stack arrays
 684   //
 685
 686   if (fParticles)
 687     fParticles->Clear();
 688   else
 689     fParticles = new TClonesArray("TParticle",1000);
 690   if (fParticleMap) {
 691     fParticleMap->Clear();
 692     if (size>0) fParticleMap->Expand(size);}
 693   else
 694     fParticleMap = new TObjArray(size);
 695 }
 696
 697 //_____________________________________________________________________________
 698 void AliStack::SetHighWaterMark(Int_t)
 699 {
 700   //
 701   // Set high water mark for last track in event
 702   //
 703
 704
 705   fHgwmk = fNtrack-1;
 706   fCurrentPrimary=fHgwmk;
 707   // Set also number of primary tracks
 708   fNprimary = fHgwmk+1;
 709   fNtrack   = fHgwmk+1;
 710 }
 711
 712 //_____________________________________________________________________________
 713 TParticle* AliStack::Particle(Int_t i)
 714 {
 715   //
 716   // Return particle with specified ID
 717
 718   //PH  if(!(*fParticleMap)[i]) {
 719   if(!fParticleMap->At(i)) {
 720     Int_t nentries = fParticles->GetEntriesFast();
 721     // algorithmic way of getting entry index
 722     // (primary particles are filled after secondaries)
 723     Int_t entry = TreeKEntry(i);
 724     // check whether algorithmic way and
 725     // and the fParticleFileMap[i] give the same;
 726     // give the fatal error if not
 727     if (entry != fParticleFileMap[i]) {
 728       AliFatal(Form(
 729         "!! The algorithmic way and map are different: !!\n entry: %d map: %d",
 730         entry, fParticleFileMap[i]));
 731     }
 732
 733     TreeK()->GetEntry(entry);
 734     new ((*fParticles)[nentries]) TParticle(*fParticleBuffer);
 735     fParticleMap->AddAt((*fParticles)[nentries],i);
 736   }
 737   //PH  return dynamic_cast<TParticle *>((*fParticleMap)[i]);
 738   return dynamic_cast<TParticle*>(fParticleMap->At(i));
 739 }
 740
 741 //_____________________________________________________________________________
 742 TParticle* AliStack::ParticleFromTreeK(Int_t id) const
 743 {
 744 //
 745 // return pointer to TParticle with label id
 746 //
 747   Int_t entry;
 748   if ((entry = TreeKEntry(id)) < 0) return 0;
 749   if (fTreeK->GetEntry(entry)<=0) return 0;
 750   return fParticleBuffer;
 751 }
 752
 753 //_____________________________________________________________________________
 754 Int_t AliStack::TreeKEntry(Int_t id) const
 755 {
 756 //
 757 // return entry number in the TreeK for particle with label id
 758 // return negative number if label>fNtrack
 759 //
 760   Int_t entry;
 761   if (id<fNprimary)
 762     entry = id+fNtrack-fNprimary;
 763   else
 764     entry = id-fNprimary;
 765   return entry;
 766 }
 767
 768 //_____________________________________________________________________________
 769 Int_t AliStack::GetCurrentParentTrackNumber() const
 770 {
 771   //
 772   // Return number of the parent of the current track
 773   //
 774
 775   TParticle* current = (TParticle*)fParticleMap->At(fCurrent);
 776
 777   if (current)
 778     return current->GetFirstMother();
 779   else {
 780     AliWarning("Current track not found in the stack");
 781     return -1;
 782   }
 783 }
 784
 785 //_____________________________________________________________________________
 786 Int_t AliStack::GetPrimary(Int_t id)
 787 {
 788   //
 789   // Return number of primary that has generated track
 790   //
 791
 792   int current, parent;
 793   //
 794   parent=id;
 795   while (1) {
 796     current=parent;
 797     parent=Particle(current)->GetFirstMother();
 798     if(parent<0) return current;
 799   }
 800 }
 801
 802 //_____________________________________________________________________________
 803 void AliStack::DumpPart (Int_t i) const
 804 {
 805   //
 806   // Dumps particle i in the stack
 807   //
 808
 809   //PH  dynamic_cast<TParticle*>((*fParticleMap)[i])->Print();
 810   dynamic_cast<TParticle*>(fParticleMap->At(i))->Print();
 811 }
 812
 813 //_____________________________________________________________________________
 814 void AliStack::DumpPStack ()
 815 {
 816   //
 817   // Dumps the particle stack
 818   //
 819
 820   Int_t i;
 821
 822   printf("\n\n=======================================================================\n");
 823   for (i=0;i<fNtrack;i++)
 824     {
 825       TParticle* particle = Particle(i);
 826       if (particle) {
 827         printf("-> %d ",i); particle->Print();
 828         printf("--------------------------------------------------------------\n");
 829       }
 830       else
 831         Warning("DumpPStack", "No particle with id %d.", i);
 832     }
 833
 834   printf("\n=======================================================================\n\n");
 835
 836   // print  particle file map
 837   printf("\nParticle file map: \n");
 838   for (i=0; i<fNtrack; i++)
 839       printf("   %d th entry: %d \n",i,fParticleFileMap[i]);
 840 }
 841
 842
 843 //_____________________________________________________________________________
 844 void AliStack::DumpLoadedStack() const
 845 {
 846   //
 847   // Dumps the particle in the stack
 848   // that are loaded in memory.
 849   //
 850
 851   TObjArray &particles = *fParticleMap;
 852   printf(
 853          "\n\n=======================================================================\n");
 854   for (Int_t i=0;i<fNtrack;i++)
 855     {
 856       TParticle* particle = dynamic_cast<TParticle*>(particles[i]);
 857       if (particle) {
 858         printf("-> %d ",i); particle->Print();
 859         printf("--------------------------------------------------------------\n");
 860       }
 861       else {
 862         printf("-> %d  Particle not loaded.\n",i);
 863         printf("--------------------------------------------------------------\n");
 864       }
 865     }
 866   printf(
 867          "\n=======================================================================\n\n");
 868 }
 869
 870 //
 871 // protected methods
 872 //
 873
 874 //_____________________________________________________________________________
 875 void AliStack::CleanParents()
 876 {
 877   //
 878   // Clean particles stack
 879   // Set parent/daughter relations
 880   //
 881
 882   TObjArray &particles = *fParticleMap;
 883   TParticle *part;
 884   int i;
 885   for(i=0; i<fHgwmk+1; i++) {
 886     part = dynamic_cast<TParticle*>(particles.At(i));
 887     if(part) if(!part->TestBit(kDaughtersBit)) {
 888       part->SetFirstDaughter(-1);
 889       part->SetLastDaughter(-1);
 890     }
 891   }
 892 }
 893
 894 //_____________________________________________________________________________
 895 TParticle* AliStack::GetNextParticle()
 896 {
 897   //
 898   // Return next particle from stack of particles
 899   //
 900
 901   TParticle* particle = 0;
 902
 903   // search secondaries
 904   //for(Int_t i=fNtrack-1; i>=0; i--) {
 905   for(Int_t i=fNtrack-1; i>fHgwmk; i--) {
 906       particle = dynamic_cast<TParticle*>(fParticleMap->At(i));
 907       if ((particle) && (!particle->TestBit(kDoneBit))) {
 908           fCurrent=i;
 909           return particle;
 910       }
 911   }
 912
 913   // take next primary if all secondaries were done
 914   while (fCurrentPrimary>=0) {
 915       fCurrent = fCurrentPrimary;
 916       particle = dynamic_cast<TParticle*>(fParticleMap->At(fCurrentPrimary--));
 917       if ((particle) && (!particle->TestBit(kDoneBit))) {
 918           return particle;
 919       }
 920   }
 921
 922   // nothing to be tracked
 923   fCurrent = -1;
 924
 925
 926   return particle;
 927 }
 928 //__________________________________________________________________________________________
 929
 930 TTree* AliStack::TreeK()
 931 {
 932 //returns TreeK
 933   if (fTreeK)
 934    {
 935      return fTreeK;
 936    }
 937   else
 938    {
 939      AliRunLoader *rl = AliRunLoader::GetRunLoader(fEventFolderName);
 940      if (rl == 0x0)
 941       {
 942         AliFatal(Form("Can not get RunLoader from event folder named %s",fEventFolderName.Data()));
 943         return 0x0;//pro forma
 944       }
 945     fTreeK = rl->TreeK();
 946     if ( fTreeK )
 947      {
 948       ConnectTree();
 949      }
 950     else
 951      {
 952       //don't panic - could be Lego
 953       AliWarning(Form("Can not get TreeK from RL. Ev. Folder is %s",fEventFolderName.Data()));
 954      }
 955    }
 956   return fTreeK;//never reached
 957 }
 958 //__________________________________________________________________________________________
 959
 960 void AliStack::ConnectTree()
 961 {
 962 //
 963 //  Creates branch for writing particles
 964 //
 965   AliDebug(1, "Connecting TreeK");
 966   if (fTreeK == 0x0)
 967    {
 968     if (TreeK() == 0x0)
 969      {
 970       AliFatal("Parameter is NULL");//we don't like such a jokes
 971       return;
 972      }
 973     return;//in this case TreeK() calls back this method (ConnectTree)
 974            //tree after setting fTreeK, the rest was already executed
 975            //it is safe to return now
 976    }
 977
 978  //  Create a branch for particles
 979
 980   AliDebug(2, Form("Tree name is %s",fTreeK->GetName()));
 981
 982   if (fTreeK->GetDirectory())
 983    {
 984      AliDebug(2, Form("and dir is %s",fTreeK->GetDirectory()->GetName()));
 985    }
 986   else
 987     AliWarning("DIR IS NOT SET !!!");
 988
 989   TBranch *branch=fTreeK->GetBranch(AliRunLoader::GetKineBranchName());
 990   if(branch == 0x0)
 991    {
 992     branch = fTreeK->Branch(AliRunLoader::GetKineBranchName(), "TParticle", &fParticleBuffer, 4000);
 993     AliDebug(2, "Creating Branch in Tree");
 994    }
 995   else
 996    {
 997     AliDebug(2, "Branch Found in Tree");
 998     branch->SetAddress(&fParticleBuffer);
 999    }
1000   if (branch->GetDirectory())
1001    {
1002     AliDebug(1, Form("Branch Dir Name is %s",branch->GetDirectory()->GetName()));
1003    }
1004   else
1005     AliWarning("Branch Dir is NOT SET");
1006 }
1007 //__________________________________________________________________________________________
1008
1009
1010 void AliStack::BeginEvent()
1011 {
1012 // start a new event
1013  Reset();
1014 }
1015
1016 //_____________________________________________________________________________
1017 void AliStack::FinishRun()
1018 {
1019 // Clean TreeK information
1020 }
1021 //_____________________________________________________________________________
1022
1023 Bool_t AliStack::GetEvent()
1024 {
1025 //
1026 // Get new event from TreeK
1027
1028     // Reset/Create the particle stack
1029     fTreeK = 0x0;
1030
1031     if (TreeK() == 0x0) //forces connecting
1032      {
1033       AliError("cannot find Kine Tree for current event");
1034       return kFALSE;
1035      }
1036
1037     Int_t size = (Int_t)TreeK()->GetEntries();
1038     ResetArrays(size);
1039     return kTRUE;
1040 }
1041 //_____________________________________________________________________________
1042
1043 void AliStack::SetEventFolderName(const char* foldname)
1044 {
1045  //Sets event folder name
1046  fEventFolderName = foldname;
1047 }
1048
1049 Bool_t AliStack::IsStable(Int_t pdg) const
1050 {
1051 //
1052 // Decide whether particle (pdg) is stable
1053 //
1054
1055     const Int_t kNstable = 14;
1056     Int_t i;
1057
1058     Int_t pdgStable[kNstable] = {
1059         kGamma,             // Photon
1060         kElectron,          // Electron
1061         kMuonPlus,          // Muon
1062         kPiPlus,            // Pion
1063         kKPlus,             // Kaon
1064         kProton,            // Proton
1065         kNeutron,           // Neutron
1066         kLambda0,           // Lambda_0
1067         kSigmaMinus,        // Sigma Minus
1068         kSigma0,            // Sigma_0
1069         kSigmaPlus,         // Sigma Plus
1070         3312,               // Xsi Minus
1071         3322,               // Xsi
1072         3334,               // Omega
1073     };
1074
1075     Bool_t isStable = kFALSE;
1076     for (i = 0; i < kNstable; i++) {
1077         if (pdg == TMath::Abs(pdgStable[i])) {
1078             isStable = kTRUE;
1079             break;
1080         }
1081     }
1082
1083     return isStable;
1084 }
1085
1086 Bool_t AliStack::IsPhysicalPrimary(Int_t index)
1087 {
1088     //
1089     // Test if a particle is a physical primary according to the following definition:
1090     // Particles produced in the collision including products of strong and
1091     // electromagnetic decay and excluding feed-down from weak decays of strange
1092     // particles.
1093     //
1094     TParticle* p = Particle(index);
1095     Int_t ist = p->GetStatusCode();
1096
1097     //
1098     // Initial state particle
1099     if (ist > 20) return kFALSE;
1100
1101     Int_t pdg = TMath::Abs(p->GetPdgCode());
1102
1103     if (!IsStable(pdg)) return kFALSE;
1104
1105     if (index < GetNprimary()) {
1106 //
1107 // Particle produced by generator
1108         return kTRUE;
1109     } else {
1110 //
1111 // Particle produced during transport
1112 //
1113 // Check if this is a heavy flavor decay product
1114         Int_t imo =  p->GetFirstMother();
1115         TParticle* pm  = Particle(imo);
1116         Int_t mpdg = TMath::Abs(pm->GetPdgCode());
1117         Int_t mfl  = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
1118         //
1119         // Light hadron
1120         if (mfl < 4) return kFALSE;
1121
1122         //
1123         // Heavy flavor hadron produced by generator
1124         if (imo <  GetNprimary()) {
1125             return kTRUE;
1126         }
1127
1128         // To be sure that heavy flavor has not been produced in a secondary interaction
1129         // Loop back to the generated mother
1130         while (imo >=  GetNprimary()) {
1131             imo = p->GetFirstMother();
1132             pm  =  Particle(imo);
1133         }
1134         mpdg = TMath::Abs(pm->GetPdgCode());
1135         mfl  = Int_t (mpdg / TMath::Power(10, Int_t(TMath::Log10(mpdg))));
1136
1137         if (mfl < 4) {
1138             return kFALSE;
1139         } else {
1140             return kTRUE;
1141         }
1142     } // produced by generator ?
1143 }
1144