TUHKMgen/UHKM/InitialState.cxx

   1 /*
   2   Ludmila Malinina  malinina@lav01.sinp.msu.ru,   SINP MSU/Moscow and JINR/Dubna
   3   Ionut Arsene  i.c.arsene@fys.uio.no,            Oslo University and ISS-Bucharest
   4   Date        : 2007/05/30
   5 */
   6
   7 #include "TVector3.h"
   8 #ifndef INITIAL_STATE
   9 #include "InitialState.h"
  10 #endif
  11 #ifndef HADRONDECAYER_INCLUDED
  12 #include "HadronDecayer.h"
  13 #endif
  14 #include <iostream>
  15 #include <fstream>
  16
  17 // aic(2008/08/08): define the fLastIndex static variable
  18 Int_t Particle::fLastIndex;
  19
  20 void InitialState::Evolve(List_t &source, List_t &secondaries, ParticleAllocator &allocator, Double_t weakDecayLimit) {
  21   LPIT_t it;
  22   LPIT_t e;
  23
  24   // Initialize the fLastIndex to -1
  25   it = source.begin();
  26   it->InitIndexing();
  27
  28   // Loop over source list (primary particles)
  29   for(it = source.begin(), e = source.end(); it != e; ++it) {
  30     // calculate the decay time of the particle
  31     Double_t decayTime = GetDecayTime(*it, weakDecayLimit);
  32     TVector3 shift(it->Mom().BoostVector());
  33     shift *= decayTime;
  34     it->SetLastInterTime(it->T() + decayTime);
  35
  36     // if the particle is unstable run the decayer
  37     if(decayTime > 0.) {
  38       it->Pos(it->Pos() += TLorentzVector(shift, 0.));
  39       it->T(it->T() + decayTime);
  40       // perform the decay procedure
  41       // the primaries are added to the list of secondaries inside Decay()
  42       Decay(secondaries, *it, allocator, fDatabase);
  43     }
  44     // else just add the primary to the list of particles
  45     else {
  46       it->SetIndex();
  47       allocator.AddParticle(*it, secondaries);
  48     }
  49   }
  50
  51   // Loop over the secondaries list and decay the cascades
  52   for(it = secondaries.begin(), e = secondaries.end(); it != e;) {
  53     // the primaries are already decayed, so just ignore it
  54     if(it->GetMother()==-1) {
  55       ++it;
  56       continue;
  57     }
  58     Double_t decayTime = GetDecayTime(*it, weakDecayLimit);
  59     it->SetLastInterTime(it->T() + decayTime);
  60     TVector3 shift(it->Mom().BoostVector());
  61     shift *= decayTime;
  62
  63     // if the particle is unstable run the decays
  64     if(decayTime > 0.) {
  65       it->Pos(it->Pos() += TLorentzVector(shift, 0.));
  66       it->T(it->T() + decayTime);
  67       // The decay products are added to the list inside the Decay() function
  68       // The mother particle is not removed from list but can be distinguished
  69       // since the daughter indexes will be initialized
  70       Decay(secondaries, *it, allocator, fDatabase);
  71       ++it;
  72     }
  73     // else just continue
  74     else
  75       ++it;
  76   }
  77 }