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