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   fCurrentTrack(0),
  51   fTreeK(0),
  52   fNtrack(0),
  53   fNprimary(0),
  54   fCurrent(-1),
  55   fCurrentPrimary(-1),
  56   fHgwmk(0),
  57   fLoadPoint(0),
  58   fEventFolderName(AliConfig::GetDefaultEventFolderName())
  59 {
  60   //
  61   // Default constructor
  62   //
  63 }
  64
  65 //_______________________________________________________________________
  66 AliStack::AliStack(Int_t size, const char* evfoldname):
  67   fParticles(new TClonesArray("TParticle",1000)),
  68   fParticleMap(new TObjArray(size)),
  69   fParticleFileMap(0),
  70   fParticleBuffer(0),
  71   fCurrentTrack(0),
  72   fTreeK(0),
  73   fNtrack(0),
  74   fNprimary(0),
  75   fCurrent(-1),
  76   fCurrentPrimary(-1),
  77   fHgwmk(0),
  78   fLoadPoint(0),
  79   fEventFolderName(evfoldname)
  80 {
  81   //
  82   //  Constructor
  83   //
  84 }
  85
  86 //_______________________________________________________________________
  87 AliStack::AliStack(const AliStack& st):
  88   TVirtualMCStack(st),
  89   fParticles(0),
  90   fParticleMap(0),
  91   fParticleFileMap(0),
  92   fParticleBuffer(0),
  93   fCurrentTrack(0),
  94   fTreeK(0),
  95   fNtrack(0),
  96   fNprimary(0),
  97   fCurrent(-1),
  98   fCurrentPrimary(-1),
  99   fHgwmk(0),
 100   fLoadPoint(0),
 101   fEventFolderName()
 102 {
 103   //
 104   // Copy constructor
 105   //
 106   st.Copy(*this);
 107 }
 108
 109 //_______________________________________________________________________
 110 void AliStack::Copy(TObject&) const
 111 {
 112   AliFatal("Not implemented!");
 113 }
 114
 115 //_______________________________________________________________________
 116 AliStack::~AliStack()
 117 {
 118   //
 119   // Destructor
 120   //
 121
 122   if (fParticles) {
 123     fParticles->Delete();
 124     delete fParticles;
 125   }
 126   delete fParticleMap;
 127   //PH???  delete fTreeK;
 128 }
 129
 130 //
 131 // public methods
 132 //
 133
 134 //_____________________________________________________________________________
 135 void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
 136                         Float_t *vpos, Float_t *polar, Float_t tof,
 137                         TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is)
 138 {
 139   //
 140   // Load a track on the stack
 141   //
 142   // done     0 if the track has to be transported
 143   //          1 if not
 144   // parent   identifier of the parent track. -1 for a primary
 145   // pdg    particle code
 146   // pmom     momentum GeV/c
 147   // vpos     position
 148   // polar    polarisation
 149   // tof      time of flight in seconds
 150   // mecha    production mechanism
 151   // ntr      on output the number of the track stored
 152   //
 153
 154   //  const Float_t tlife=0;
 155
 156   //
 157   // Here we get the static mass
 158   // For MC is ok, but a more sophisticated method could be necessary
 159   // if the calculated mass is required
 160   // also, this method is potentially dangerous if the mass
 161   // used in the MC is not the same of the PDG database
 162   //
 163     TParticlePDG* pmc =  TDatabasePDG::Instance()->GetParticle(pdg);
 164     if (pmc) {
 165         Float_t mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
 166         Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+
 167                               pmom[1]*pmom[1]+pmom[2]*pmom[2]);
 168
 169 //    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",
 170 //         mass,e,fNtrack,pdg,parent,done,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS);
 171
 172
 173         PushTrack(done, parent, pdg, pmom[0], pmom[1], pmom[2], e,
 174                  vpos[0], vpos[1], vpos[2], tof, polar[0], polar[1], polar[2],
 175                  mech, ntr, weight, is);
 176     } else {
 177         AliWarning(Form("Particle type %d not defined in PDG Database !", pdg));
 178         AliWarning("Particle skipped !");
 179     }
 180 }
 181
 182 //_____________________________________________________________________________
 183 void AliStack::PushTrack(Int_t done, Int_t parent, Int_t pdg,
 184                       Double_t px, Double_t py, Double_t pz, Double_t e,
 185                       Double_t vx, Double_t vy, Double_t vz, Double_t tof,
 186                       Double_t polx, Double_t poly, Double_t polz,
 187                       TMCProcess mech, Int_t &ntr, Double_t weight, Int_t is)
 188 {
 189   //
 190   // Load a track on the stack
 191   //
 192   // done        0 if the track has to be transported
 193   //             1 if not
 194   // parent      identifier of the parent track. -1 for a primary
 195   // pdg         particle code
 196   // kS          generation status code
 197   // px, py, pz  momentum GeV/c
 198   // vx, vy, vz  position
 199   // polar       polarisation
 200   // tof         time of flight in seconds
 201   // mech        production mechanism
 202   // ntr         on output the number of the track stored
 203   //
 204   // New method interface:
 205   // arguments were changed to be in correspondence with TParticle
 206   // constructor.
 207   // Note: the energy is not calculated from the static mass but
 208   // it is passed by argument e.
 209
 210
 211   const Int_t kFirstDaughter=-1;
 212   const Int_t kLastDaughter=-1;
 213
 214
 215   TClonesArray &particles = *fParticles;
 216   TParticle* particle
 217     = new(particles[fLoadPoint++])
 218       TParticle(pdg, is, parent, -1, kFirstDaughter, kLastDaughter,
 219                 px, py, pz, e, vx, vy, vz, tof);
 220
 221   particle->SetPolarisation(polx, poly, polz);
 222   particle->SetWeight(weight);
 223   particle->SetUniqueID(mech);
 224
 225   if(!done) particle->SetBit(kDoneBit);
 226
 227   //  Declare that the daughter information is valid
 228   particle->SetBit(kDaughtersBit);
 229   //  Add the particle to the stack
 230
 231
 232   fParticleMap->AddAtAndExpand(particle, fNtrack);//CHECK!!
 233
 234   if(parent>=0) {
 235     particle = dynamic_cast<TParticle*>(fParticleMap->At(parent));
 236     if (particle) {
 237       particle->SetLastDaughter(fNtrack);
 238       if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack);
 239     }
 240     else {
 241       AliError(Form("Parent %d does not exist",parent));
 242     }
 243   }
 244   else {
 245     //
 246     // This is a primary track. Set high water mark for this event
 247     fHgwmk = fNtrack;
 248     //
 249     // Set also number if primary tracks
 250     fNprimary = fHgwmk+1;
 251     fCurrentPrimary++;
 252   }
 253   ntr = fNtrack++;
 254 }
 255
 256 //_____________________________________________________________________________
 257 TParticle*  AliStack::PopNextTrack(Int_t& itrack)
 258 {
 259   //
 260   // Returns next track from stack of particles
 261   //
 262
 263
 264   TParticle* track = GetNextParticle();
 265
 266   if (track) {
 267     itrack = fCurrent;
 268     track->SetBit(kDoneBit);
 269   }
 270   else
 271     itrack = -1;
 272
 273   fCurrentTrack = track;
 274   return track;
 275 }
 276
 277 //_____________________________________________________________________________
 278 TParticle*  AliStack::PopPrimaryForTracking(Int_t i)
 279 {
 280   //
 281   // Returns i-th primary particle if it is flagged to be tracked,
 282   // 0 otherwise
 283   //
 284
 285   TParticle* particle = Particle(i);
 286
 287   if (!particle->TestBit(kDoneBit))
 288     return particle;
 289   else
 290     return 0;
 291 }
 292
 293 //_____________________________________________________________________________
 294 void AliStack::PurifyKine()
 295 {
 296   //
 297   // Compress kinematic tree keeping only flagged particles
 298   // and renaming the particle id's in all the hits
 299   //
 300
 301   TObjArray &particles = *fParticleMap;
 302   int nkeep=fHgwmk+1, parent, i;
 303   TParticle *part, *father;
 304   TArrayI map(particles.GetLast()+1);
 305
 306   // Save in Header total number of tracks before compression
 307   // If no tracks generated return now
 308   if(fHgwmk+1 == fNtrack) return;
 309
 310   // First pass, invalid Daughter information
 311
 312   for(i=0; i<fNtrack; i++) {
 313     // Preset map, to be removed later
 314     if(i<=fHgwmk) map[i]=i ;
 315     else {
 316       map[i] = -99;
 317       if((part=dynamic_cast<TParticle*>(particles.At(i)))) {
 318 //
 319 //        Check of this track should be kept for physics reasons
 320           if (KeepPhysics(part)) KeepTrack(i);
 321 //
 322           part->ResetBit(kDaughtersBit);
 323           part->SetFirstDaughter(-1);
 324           part->SetLastDaughter(-1);
 325       }
 326     }
 327   }
 328   // Invalid daughter information for the parent of the first particle
 329   // generated. This may or may not be the current primary according to
 330   // whether decays have been recorded among the primaries
 331   part = dynamic_cast<TParticle*>(particles.At(fHgwmk+1));
 332   particles.At(part->GetFirstMother())->ResetBit(kDaughtersBit);
 333   // Second pass, build map between old and new numbering
 334   for(i=fHgwmk+1; i<fNtrack; i++) {
 335     if(particles.At(i)->TestBit(kKeepBit)) {
 336
 337       // This particle has to be kept
 338       map[i]=nkeep;
 339       // If old and new are different, have to move the pointer
 340       if(i!=nkeep) particles[nkeep]=particles.At(i);
 341       part = dynamic_cast<TParticle*>(particles.At(nkeep));
 342
 343       // as the parent is always *before*, it must be already
 344       // in place. This is what we are checking anyway!
 345       if((parent=part->GetFirstMother())>fHgwmk)
 346         if(map[parent]==-99) Fatal("PurifyKine","map[%d] = -99!\n",parent);
 347         else part->SetFirstMother(map[parent]);
 348
 349       nkeep++;
 350     }
 351   }
 352
 353   // Fix daughters information
 354   for (i=fHgwmk+1; i<nkeep; i++) {
 355     part = dynamic_cast<TParticle*>(particles.At(i));
 356     parent = part->GetFirstMother();
 357     if(parent>=0) {
 358       father = dynamic_cast<TParticle*>(particles.At(parent));
 359       if(father->TestBit(kDaughtersBit)) {
 360
 361         if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i);
 362         if(i>father->GetLastDaughter())  father->SetLastDaughter(i);
 363       } else {
 364         // Initialise daughters info for first pass
 365         father->SetFirstDaughter(i);
 366         father->SetLastDaughter(i);
 367         father->SetBit(kDaughtersBit);
 368       }
 369     }
 370   }
 371
 372   // Now loop on all registered hit lists
 373   TList* hitLists = gAlice->GetMCApp()->GetHitLists();
 374   TIter next(hitLists);
 375   TCollection *hitList;
 376
 377   while((hitList = dynamic_cast<TCollection*>(next()))) {
 378     TIter nexthit(hitList);
 379     AliHit *hit;
 380     while((hit = dynamic_cast<AliHit*>(nexthit()))) {
 381         hit->SetTrack(map[hit->GetTrack()]);
 382     }
 383   }
 384
 385   //
 386   // This for detectors which have a special mapping mechanism
 387   // for hits, such as TPC and TRD
 388   //
 389
 390    TObjArray* modules = gAlice->Modules();
 391    TIter nextmod(modules);
 392    AliModule *detector;
 393    while((detector = dynamic_cast<AliModule*>(nextmod()))) {
 394      detector->RemapTrackHitIDs(map.GetArray());
 395      detector->RemapTrackReferencesIDs(map.GetArray());
 396    }
 397    //
 398    gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray());
 399
 400    // Now the output bit, from fHgwmk to nkeep we write everything and we erase
 401    if(nkeep>fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5));
 402
 403    for (i=fHgwmk+1; i<nkeep; ++i) {
 404      fParticleBuffer = dynamic_cast<TParticle*>(particles.At(i));
 405      fParticleFileMap[i]=static_cast<Int_t>(TreeK()->GetEntries());
 406      TreeK()->Fill();
 407      particles[i]=fParticleBuffer=0;
 408     }
 409
 410    for (i=nkeep; i<fNtrack; ++i) particles[i]=0;
 411
 412    Int_t toshrink = fNtrack-fHgwmk-1;
 413    fLoadPoint-=toshrink;
 414
 415
 416    for(i=fLoadPoint; i<fLoadPoint+toshrink; ++i) fParticles->RemoveAt(i);
 417
 418    fNtrack=nkeep;
 419    fHgwmk=nkeep-1;
 420 }
 421
 422 void AliStack::ReorderKine()
 423 {
 424 //
 425 // In some transport code children might not come in a continuous sequence.
 426 // In this case the stack  has  to  be reordered in order to establish the
 427 // mother daughter relation using index ranges.
 428 //
 429   if(fHgwmk+1 == fNtrack) return;
 430
 431   //
 432   // Howmany secondaries have been produced ?
 433   Int_t nNew = fNtrack - fHgwmk - 1;
 434
 435   if (nNew > 0) {
 436       Int_t i, j;
 437       TObjArray &particles = *fParticleMap;
 438       TArrayI map1(nNew);
 439       //
 440       // Copy pointers to temporary array
 441       TParticle** tmp = new TParticle*[nNew];
 442
 443       for (i = 0; i < nNew; i++) {
 444           if (particles.At(fHgwmk + 1 + i)) {
 445 //            tmp[i] = (TParticle*) (particles.At(fHgwmk + 1 + i))->Clone();
 446               tmp[i] = (TParticle*) (particles.At(fHgwmk + 1 + i));
 447
 448 //            if (((TParticle*) (particles.At(fHgwmk + 1 + i)))->TestBit(kKeepBit))
 449 //                tmp[i]->SetBit(kKeepBit);
 450 //            if (((TParticle*) (particles.At(fHgwmk + 1 + i)))->TestBit(kDoneBit))
 451 //                tmp[i]->SetBit(kDoneBit);
 452           } else {
 453               tmp[i] = 0x0;
 454           }
 455           map1[i] = -99;
 456       }
 457
 458
 459       //
 460       // Reset  LoadPoint
 461       //
 462       fLoadPoint = fHgwmk + 1;
 463       //
 464       // Re-Push particles into stack
 465       // The outer loop is over parents, the inner over children.
 466       // -1 refers to the primary particle
 467       //
 468       for (i = -1; i < nNew - 1; i++) {
 469           Int_t ipa;
 470           TParticle* parP;
 471           if (i == -1) {
 472               ipa  = tmp[0]->GetFirstMother();
 473               parP =dynamic_cast<TParticle*>(particles.At(ipa));
 474           } else {
 475               ipa = (fHgwmk + 1 + i);
 476               // Skip deleted particles
 477               if (!tmp[i])                          continue;
 478               // Skip particles without children
 479               if (tmp[i]->GetFirstDaughter() == -1) continue;
 480               parP = tmp[i];
 481           }
 482           // Reset daughter information
 483           parP->SetFirstDaughter(-1);
 484           parP->SetLastDaughter(-1);
 485           for (j = i + 1; j < nNew; 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 }