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