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