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