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