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 #include <iostream>
   7 #include <fstream>
   8 using namespace std;
   9
  10 #include <TVector3.h>
  11 #include "HadronDecayer.h"
  12 #include "InitialState.h"
  13
  14 void InitialState::Evolve(List_t &secondaries, ParticleAllocator &allocator, Double_t weakDecayLimit) {
  15   // Particle indexes are set for primaries already from InitialStateHydjet::Initialize()
  16
  17   // particle list iterators
  18   LPIT_t it;
  19   LPIT_t e;
  20
  21   // Particle indexes are set for primaries already from InitialStateHydjet::Initialize()
  22
  23   // Decay loop
  24   // Note that the decay products are always added at the end of list so the unstable products are
  25   // decayed when the iterator reaches them (used for cascade decays)
  26
  27   it = secondaries.begin();
  28
  29   for(it = secondaries.begin(), e = secondaries.end(); it != e; it++) {
  30     // if the decay procedure was applied already skip ... (e.g. particles from pythia history information)
  31     if(it->GetDecayed()) {
  32       continue;
  33     }
  34
  35     // generate the decay time; if particle is stable or set to stable decayTime=0
  36     Double_t decayTime = GetDecayTime(*it, weakDecayLimit);
  37     it->SetLastInterTime(it->T() + decayTime);
  38     TVector3 shift(it->Mom().BoostVector());
  39     shift *= decayTime;
  40     it->SetDecayed();
  41
  42     // if decayTime>0 then apply the decay procedure (only 2 or 3 body decays)
  43     if(decayTime > 0.) {
  44       it->Pos(it->Pos() += TLorentzVector(shift, 0.));
  45       it->T(it->T() + decayTime);
  46       Decay(secondaries, *it, allocator, fDatabase);
  47     }
  48     // if particle is stable just continue
  49   }
  50
  51   it = secondaries.begin();
  52   Int_t npart = it->GetLastIndex();
  53   cout << "particles generated = " << npart << endl;
  54   for(Int_t count=0; count<npart; count++, it++) {
  55 //    cout << "=====================================================" << endl;
  56 //    cout << "InitialState::Evolve() particle count = " << count << endl;
  57 //    cout << "InitialState::Evolve() particle pdg = " << it->Encoding() << endl;
  58 //    cout << "InitialState::Evolve() index = " << it->GetIndex() << endl;
  59 //    cout << "InitialState::Evolve() mother index = " << it->GetMother() << endl;
  60 //    cout << "InitialState::Evolve() mother pdg = " << it->GetLastMotherPdg() << endl;
  61 //    cout << "InitialState::Evolve() n daughters = " << it->GetNDaughters() << endl;
  62 //    cout << "InitialState::Evolve() d (first, last) = " << it->GetFirstDaughterIndex() << ", "
  63 //       << it->GetLastDaughterIndex() << endl;
  64   }
  65 }