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