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