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