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