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