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
 373   // Now loop on all registered hit lists
 374   TList* hitLists = gAlice->GetMCApp()->GetHitLists();
 375   TIter next(hitLists);
 376   TCollection *hitList;
 377
 378   while((hitList = dynamic_cast<TCollection*>(next()))) {
 379     TIter nexthit(hitList);
 380     AliHit *hit;
 381     while((hit = dynamic_cast<AliHit*>(nexthit()))) {
 382         hit->SetTrack(map[hit->GetTrack()]);
 383     }
 384   }
 385
 386   //
 387   // This for detectors which have a special mapping mechanism
 388   // for hits, such as TPC and TRD
 389   //
 390
 391    TObjArray* modules = gAlice->Modules();
 392    TIter nextmod(modules);
 393    AliModule *detector;
 394    while((detector = dynamic_cast<AliModule*>(nextmod()))) {
 395      detector->RemapTrackHitIDs(map.GetArray());
 396      detector->RemapTrackReferencesIDs(map.GetArray());
 397    }
 398    //
 399    gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray());
 400
 401    // Now the output bit, from fHgwmk to nkeep we write everything and we erase
 402    if(nkeep>fParticleFileMap.GetSize()) fParticleFileMap.Set(Int_t (nkeep*1.5));
 403
 404    for (i=fHgwmk+1; i<nkeep; ++i) {
 405      fParticleBuffer = dynamic_cast<TParticle*>(particles.At(i));
 406      fParticleFileMap[i]=static_cast<Int_t>(TreeK()->GetEntries());
 407      TreeK()->Fill();
 408      particles[i]=fParticleBuffer=0;
 409     }
 410
 411    for (i=nkeep; i<fNtrack; ++i) particles[i]=0;
 412
 413    Int_t toshrink = fNtrack-fHgwmk-1;
 414    fLoadPoint-=toshrink;
 415
 416
 417    for(i=fLoadPoint; i<fLoadPoint+toshrink; ++i) fParticles->RemoveAt(i);
 418
 419    fNtrack=nkeep;
 420    fHgwmk=nkeep-1;
 421 }
 422
 423 void AliStack::ReorderKine()
 424 {
 425 //
 426 // In some transport code children might not come in a continuous sequence.
 427 // In this case the stack  has  to  be reordered in order to establish the
 428 // mother daughter relation using index ranges.
 429 //
 430   if(fHgwmk+1 == fNtrack) return;
 431
 432   //
 433   // Howmany secondaries have been produced ?
 434   Int_t nNew = fNtrack - fHgwmk - 1;
 435
 436   if (nNew > 0) {
 437       Int_t i, j;
 438       TObjArray &particles = *fParticleMap;
 439       TArrayI map1(nNew);
 440       //
 441       // Copy pointers to temporary array
 442       TParticle** tmp = new TParticle*[nNew];
 443
 444       for (i = 0; i < nNew; i++) {
 445           if (particles.At(fHgwmk + 1 + i)) {
 446               tmp[i] = (TParticle*) (particles.At(fHgwmk + 1 + i));
 447           } else {
 448               tmp[i] = 0x0;
 449           }
 450           map1[i] = -99;
 451       }
 452
 453
 454       //
 455       // Reset  LoadPoint
 456       //
 457       fLoadPoint = fHgwmk + 1;
 458       //
 459       // Re-Push particles into stack
 460       // The outer loop is over parents, the inner over children.
 461       // -1 refers to the primary particle
 462       //
 463       for (i = -1; i < nNew-1; i++) {
 464           Int_t ipa;
 465           TParticle* parP;
 466           if (i == -1) {
 467               ipa  = tmp[0]->GetFirstMother();
 468               parP =dynamic_cast<TParticle*>(particles.At(ipa));
 469           } else {
 470               ipa = (fHgwmk + 1 + i);
 471               // Skip deleted particles
 472               if (!tmp[i])                          continue;
 473               // Skip particles without children
 474               if (tmp[i]->GetFirstDaughter() == -1) continue;
 475               parP = tmp[i];
 476           }
 477           // Reset daughter information
 478
 479           Int_t idaumin = parP->GetFirstDaughter() - fHgwmk - 1;
 480           Int_t idaumax = parP->GetLastDaughter()  - fHgwmk - 1;
 481           parP->SetFirstDaughter(-1);
 482           parP->SetLastDaughter(-1);
 483           for (j = idaumin; j <= idaumax; j++) {
 484               // Skip deleted particles
 485               if (!tmp[j])        continue;
 486               // Skip particles already handled
 487               if (map1[j] != -99) continue;
 488               Int_t jpa = tmp[j]->GetFirstMother();
 489               // Check if daughter of current parent
 490               if (jpa == ipa) {
 491                   particles[fLoadPoint] = tmp[j];
 492                   // Re-establish daughter information
 493                   parP->SetLastDaughter(fLoadPoint);
 494                   if (parP->GetFirstDaughter() == -1) parP->SetFirstDaughter(fLoadPoint);
 495                   // Set Mother information
 496                   if (i != -1) {
 497                       tmp[j]->SetFirstMother(map1[i]);
 498                   }
 499                   // Build the map
 500                   map1[j] = fLoadPoint;
 501                   // Increase load point
 502                   fLoadPoint++;
 503               }
 504           } // children
 505       } // parents
 506
 507       delete[] tmp;
 508
 509       //
 510       // Build map for remapping of hits
 511       //
 512       TArrayI map(fNtrack);
 513       for (i = 0; i < fNtrack; i ++) {
 514           if (i <= fHgwmk) {
 515               map[i] = i;
 516           } else{
 517               map[i] = map1[i - fHgwmk -1];
 518           }
 519       }
 520
 521       // Now loop on all registered hit lists
 522
 523       TList* hitLists = gAlice->GetMCApp()->GetHitLists();
 524       TIter next(hitLists);
 525       TCollection *hitList;
 526
 527       while((hitList = dynamic_cast<TCollection*>(next()))) {
 528           TIter nexthit(hitList);
 529           AliHit *hit;
 530           while((hit = dynamic_cast<AliHit*>(nexthit()))) {
 531               hit->SetTrack(map[hit->GetTrack()]);
 532           }
 533       }
 534
 535   //
 536   // This for detectors which have a special mapping mechanism
 537   // for hits, such as TPC and TRD
 538   //
 539
 540       TObjArray* modules = gAlice->Modules();
 541       TIter nextmod(modules);
 542       AliModule *detector;
 543       while((detector = dynamic_cast<AliModule*>(nextmod()))) {
 544           detector->RemapTrackHitIDs(map.GetArray());
 545           detector->RemapTrackReferencesIDs(map.GetArray());
 546       }
 547       gAlice->GetMCApp()->RemapTrackReferencesIDs(map.GetArray());
 548   } // new particles poduced
 549 }
 550
 551 Bool_t AliStack::KeepPhysics(TParticle* part)
 552 {
 553     //
 554     // Some particles have to kept on the stack for reasons motivated
 555     // by physics analysis. Decision is put here.
 556     //
 557     Bool_t keep = kFALSE;
 558     //
 559     // Keep first-generation daughter from primaries with heavy flavor
 560     //
 561     Int_t parent = part->GetFirstMother();
 562     if (parent >= 0 && parent <= fHgwmk) {
 563         TParticle* father = dynamic_cast<TParticle*>(Particles()->At(parent));
 564         Int_t kf = father->GetPdgCode();
 565         kf = TMath::Abs(kf);
 566         Int_t kfl = kf;
 567         // meson ?
 568         if  (kfl > 10) kfl/=100;
 569         // baryon
 570         if (kfl > 10) kfl/=10;
 571         if (kfl > 10) kfl/=10;
 572         if (kfl >= 4) {
 573             keep = kTRUE;
 574         }
 575     }
 576     return keep;
 577 }
 578
 579 //_____________________________________________________________________________
 580 void AliStack::FinishEvent()
 581 {
 582 //
 583 // Write out the kinematics that was not yet filled
 584 //
 585
 586 // Update event header
 587
 588
 589   if (!TreeK()) {
 590 //    Fatal("FinishEvent", "No kinematics tree is defined.");
 591 //    Don't panic this is a probably a lego run
 592       return;
 593   }
 594
 595   CleanParents();
 596   if(TreeK()->GetEntries() ==0) {
 597     // set the fParticleFileMap size for the first time
 598     fParticleFileMap.Set(fHgwmk+1);
 599   }
 600
 601   Bool_t allFilled = kFALSE;
 602   TObject *part;
 603   for(Int_t i=0; i<fHgwmk+1; ++i)
 604     if((part=fParticleMap->At(i))) {
 605       fParticleBuffer = dynamic_cast<TParticle*>(part);
 606       fParticleFileMap[i]= static_cast<Int_t>(TreeK()->GetEntries());
 607       TreeK()->Fill();
 608       //PH      (*fParticleMap)[i]=fParticleBuffer=0;
 609       fParticleBuffer=0;
 610       fParticleMap->AddAt(0,i);
 611
 612       // When all primaries were filled no particle!=0
 613       // should be left => to be removed later.
 614       if (allFilled) AliWarning(Form("Why != 0 part # %d?\n",i));
 615     }
 616     else
 617      {
 618       // // printf("Why = 0 part # %d?\n",i); => We know.
 619       // break;
 620       // we don't break now in order to be sure there is no
 621       // particle !=0 left.
 622       // To be removed later and replaced with break.
 623        if(!allFilled) allFilled = kTRUE;
 624      }
 625 }
 626 //_____________________________________________________________________________
 627
 628 void AliStack::FlagTrack(Int_t track)
 629 {
 630   //
 631   // Flags a track and all its family tree to be kept
 632   //
 633
 634   TParticle *particle;
 635
 636   Int_t curr=track;
 637   while(1) {
 638     particle=dynamic_cast<TParticle*>(fParticleMap->At(curr));
 639
 640     // If the particle is flagged the three from here upward is saved already
 641     if(particle->TestBit(kKeepBit)) return;
 642
 643     // Save this particle
 644     particle->SetBit(kKeepBit);
 645
 646     // Move to father if any
 647     if((curr=particle->GetFirstMother())==-1) return;
 648   }
 649 }
 650
 651 //_____________________________________________________________________________
 652 void AliStack::KeepTrack(Int_t track)
 653 {
 654   //
 655   // Flags a track to be kept
 656   //
 657
 658   fParticleMap->At(track)->SetBit(kKeepBit);
 659 }
 660
 661 //_____________________________________________________________________________
 662 void  AliStack::Reset(Int_t size)
 663 {
 664   //
 665   // Resets stack
 666   //
 667
 668   fNtrack=0;
 669   fNprimary=0;
 670   fHgwmk=0;
 671   fLoadPoint=0;
 672   fCurrent = -1;
 673   fTreeK = 0x0;
 674   ResetArrays(size);
 675 }
 676
 677 //_____________________________________________________________________________
 678 void  AliStack::ResetArrays(Int_t size)
 679 {
 680   //
 681   // Resets stack arrays
 682   //
 683
 684   if (fParticles)
 685     fParticles->Clear();
 686   else
 687     fParticles = new TClonesArray("TParticle",1000);
 688   if (fParticleMap) {
 689     fParticleMap->Clear();
 690     if (size>0) fParticleMap->Expand(size);}
 691   else
 692     fParticleMap = new TObjArray(size);
 693 }
 694
 695 //_____________________________________________________________________________
 696 void AliStack::SetHighWaterMark(Int_t)
 697 {
 698   //
 699   // Set high water mark for last track in event
 700   //
 701
 702
 703   fHgwmk = fNtrack-1;
 704   fCurrentPrimary=fHgwmk;
 705   // Set also number of primary tracks
 706   fNprimary = fHgwmk+1;
 707   fNtrack   = fHgwmk+1;
 708 }
 709
 710 //_____________________________________________________________________________
 711 TParticle* AliStack::Particle(Int_t i)
 712 {
 713   //
 714   // Return particle with specified ID
 715
 716   //PH  if(!(*fParticleMap)[i]) {
 717   if(!fParticleMap->At(i)) {
 718     Int_t nentries = fParticles->GetEntriesFast();
 719     // algorithmic way of getting entry index
 720     // (primary particles are filled after secondaries)
 721     Int_t entry = TreeKEntry(i);
 722     // check whether algorithmic way and
 723     // and the fParticleFileMap[i] give the same;
 724     // give the fatal error if not
 725     if (entry != fParticleFileMap[i]) {
 726       AliFatal(Form(
 727         "!! The algorithmic way and map are different: !!\n entry: %d map: %d",
 728         entry, fParticleFileMap[i]));
 729     }
 730
 731     TreeK()->GetEntry(entry);
 732     new ((*fParticles)[nentries]) TParticle(*fParticleBuffer);
 733     fParticleMap->AddAt((*fParticles)[nentries],i);
 734   }
 735   //PH  return dynamic_cast<TParticle *>((*fParticleMap)[i]);
 736   return dynamic_cast<TParticle*>(fParticleMap->At(i));
 737 }
 738
 739 //_____________________________________________________________________________
 740 TParticle* AliStack::ParticleFromTreeK(Int_t id) const
 741 {
 742 //
 743 // return pointer to TParticle with label id
 744 //
 745   Int_t entry;
 746   if ((entry = TreeKEntry(id)) < 0) return 0;
 747   if (fTreeK->GetEntry(entry)<=0) return 0;
 748   return fParticleBuffer;
 749 }
 750
 751 //_____________________________________________________________________________
 752 Int_t AliStack::TreeKEntry(Int_t id) const
 753 {
 754 //
 755 // return entry number in the TreeK for particle with label id
 756 // return negative number if label>fNtrack
 757 //
 758   Int_t entry;
 759   if (id<fNprimary)
 760     entry = id+fNtrack-fNprimary;
 761   else
 762     entry = id-fNprimary;
 763   return entry;
 764 }
 765
 766 //_____________________________________________________________________________
 767 Int_t AliStack::GetCurrentParentTrackNumber() const
 768 {
 769   //
 770   // Return number of the parent of the current track
 771   //
 772
 773   TParticle* current = (TParticle*)fParticleMap->At(fCurrent);
 774
 775   if (current)
 776     return current->GetFirstMother();
 777   else {
 778     AliWarning("Current track not found in the stack");
 779     return -1;
 780   }
 781 }
 782
 783 //_____________________________________________________________________________
 784 Int_t AliStack::GetPrimary(Int_t id)
 785 {
 786   //
 787   // Return number of primary that has generated track
 788   //
 789
 790   int current, parent;
 791   //
 792   parent=id;
 793   while (1) {
 794     current=parent;
 795     parent=Particle(current)->GetFirstMother();
 796     if(parent<0) return current;
 797   }
 798 }
 799
 800 //_____________________________________________________________________________
 801 void AliStack::DumpPart (Int_t i) const
 802 {
 803   //
 804   // Dumps particle i in the stack
 805   //
 806
 807   //PH  dynamic_cast<TParticle*>((*fParticleMap)[i])->Print();
 808   dynamic_cast<TParticle*>(fParticleMap->At(i))->Print();
 809 }
 810
 811 //_____________________________________________________________________________
 812 void AliStack::DumpPStack ()
 813 {
 814   //
 815   // Dumps the particle stack
 816   //
 817
 818   Int_t i;
 819
 820   printf("\n\n=======================================================================\n");
 821   for (i=0;i<fNtrack;i++)
 822     {
 823       TParticle* particle = Particle(i);
 824       if (particle) {
 825         printf("-> %d ",i); particle->Print();
 826         printf("--------------------------------------------------------------\n");
 827       }
 828       else
 829         Warning("DumpPStack", "No particle with id %d.", i);
 830     }
 831
 832   printf("\n=======================================================================\n\n");
 833
 834   // print  particle file map
 835   printf("\nParticle file map: \n");
 836   for (i=0; i<fNtrack; i++)
 837       printf("   %d th entry: %d \n",i,fParticleFileMap[i]);
 838 }
 839
 840
 841 //_____________________________________________________________________________
 842 void AliStack::DumpLoadedStack() const
 843 {
 844   //
 845   // Dumps the particle in the stack
 846   // that are loaded in memory.
 847   //
 848
 849   TObjArray &particles = *fParticleMap;
 850   printf(
 851          "\n\n=======================================================================\n");
 852   for (Int_t i=0;i<fNtrack;i++)
 853     {
 854       TParticle* particle = dynamic_cast<TParticle*>(particles[i]);
 855       if (particle) {
 856         printf("-> %d ",i); particle->Print();
 857         printf("--------------------------------------------------------------\n");
 858       }
 859       else {
 860         printf("-> %d  Particle not loaded.\n",i);
 861         printf("--------------------------------------------------------------\n");
 862       }
 863     }
 864   printf(
 865          "\n=======================================================================\n\n");
 866 }
 867
 868 //
 869 // protected methods
 870 //
 871
 872 //_____________________________________________________________________________
 873 void AliStack::CleanParents()
 874 {
 875   //
 876   // Clean particles stack
 877   // Set parent/daughter relations
 878   //
 879
 880   TObjArray &particles = *fParticleMap;
 881   TParticle *part;
 882   int i;
 883   for(i=0; i<fHgwmk+1; i++) {
 884     part = dynamic_cast<TParticle*>(particles.At(i));
 885     if(part) if(!part->TestBit(kDaughtersBit)) {
 886       part->SetFirstDaughter(-1);
 887       part->SetLastDaughter(-1);
 888     }
 889   }
 890 }
 891
 892 //_____________________________________________________________________________
 893 TParticle* AliStack::GetNextParticle()
 894 {
 895   //
 896   // Return next particle from stack of particles
 897   //
 898
 899   TParticle* particle = 0;
 900
 901   // search secondaries
 902   //for(Int_t i=fNtrack-1; i>=0; i--) {
 903   for(Int_t i=fNtrack-1; i>fHgwmk; i--) {
 904       particle = dynamic_cast<TParticle*>(fParticleMap->At(i));
 905       if ((particle) && (!particle->TestBit(kDoneBit))) {
 906           fCurrent=i;
 907           return particle;
 908       }
 909   }
 910
 911   // take next primary if all secondaries were done
 912   while (fCurrentPrimary>=0) {
 913       fCurrent = fCurrentPrimary;
 914       particle = dynamic_cast<TParticle*>(fParticleMap->At(fCurrentPrimary--));
 915       if ((particle) && (!particle->TestBit(kDoneBit))) {
 916           return particle;
 917       }
 918   }
 919
 920   // nothing to be tracked
 921   fCurrent = -1;
 922
 923
 924   return particle;
 925 }
 926 //__________________________________________________________________________________________
 927
 928 TTree* AliStack::TreeK()
 929 {
 930 //returns TreeK
 931   if (fTreeK)
 932    {
 933      return fTreeK;
 934    }
 935   else
 936    {
 937      AliRunLoader *rl = AliRunLoader::GetRunLoader(fEventFolderName);
 938      if (rl == 0x0)
 939       {
 940         AliFatal(Form("Can not get RunLoader from event folder named %s",fEventFolderName.Data()));
 941         return 0x0;//pro forma
 942       }
 943     fTreeK = rl->TreeK();
 944     if ( fTreeK )
 945      {
 946       ConnectTree();
 947      }
 948     else
 949      {
 950       //don't panic - could be Lego
 951       AliWarning(Form("Can not get TreeK from RL. Ev. Folder is %s",fEventFolderName.Data()));
 952      }
 953    }
 954   return fTreeK;//never reached
 955 }
 956 //__________________________________________________________________________________________
 957
 958 void AliStack::ConnectTree()
 959 {
 960 //
 961 //  Creates branch for writing particles
 962 //
 963   AliDebug(1, "Connecting TreeK");
 964   if (fTreeK == 0x0)
 965    {
 966     if (TreeK() == 0x0)
 967      {
 968       AliFatal("Parameter is NULL");//we don't like such a jokes
 969       return;
 970      }
 971     return;//in this case TreeK() calls back this method (ConnectTree)
 972            //tree after setting fTreeK, the rest was already executed
 973            //it is safe to return now
 974    }
 975
 976  //  Create a branch for particles
 977
 978   AliDebug(2, Form("Tree name is %s",fTreeK->GetName()));
 979
 980   if (fTreeK->GetDirectory())
 981    {
 982      AliDebug(2, Form("and dir is %s",fTreeK->GetDirectory()->GetName()));
 983    }
 984   else
 985     AliWarning("DIR IS NOT SET !!!");
 986
 987   TBranch *branch=fTreeK->GetBranch(AliRunLoader::GetKineBranchName());
 988   if(branch == 0x0)
 989    {
 990     branch = fTreeK->Branch(AliRunLoader::GetKineBranchName(), "TParticle", &fParticleBuffer, 4000);
 991     AliDebug(2, "Creating Branch in Tree");
 992    }
 993   else
 994    {
 995     AliDebug(2, "Branch Found in Tree");
 996     branch->SetAddress(&fParticleBuffer);
 997    }
 998   if (branch->GetDirectory())
 999    {
1000     AliDebug(1, Form("Branch Dir Name is %s",branch->GetDirectory()->GetName()));
1001    }
1002   else
1003     AliWarning("Branch Dir is NOT SET");
1004 }
1005 //__________________________________________________________________________________________
1006
1007
1008 void AliStack::BeginEvent()
1009 {
1010 // start a new event
1011  Reset();
1012 }
1013
1014 //_____________________________________________________________________________
1015 void AliStack::FinishRun()
1016 {
1017 // Clean TreeK information
1018 }
1019 //_____________________________________________________________________________
1020
1021 Bool_t AliStack::GetEvent()
1022 {
1023 //
1024 // Get new event from TreeK
1025
1026     // Reset/Create the particle stack
1027     fTreeK = 0x0;
1028
1029     if (TreeK() == 0x0) //forces connecting
1030      {
1031       AliError("cannot find Kine Tree for current event");
1032       return kFALSE;
1033      }
1034
1035     Int_t size = (Int_t)TreeK()->GetEntries();
1036     ResetArrays(size);
1037     return kTRUE;
1038 }
1039 //_____________________________________________________________________________
1040
1041 void AliStack::SetEventFolderName(const char* foldname)
1042 {
1043  //Sets event folder name
1044  fEventFolderName = foldname;
1045 }