Coding conventions (Ionut)
[u/mrichter/AliRoot.git] / 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 }