STEER/AliStack.cxx

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