Fixing error on gcc 4.5.1
[u/mrichter/AliRoot.git] / TTherminator / Therminator / ParticleDecayer.cxx
CommitLineData
2e967919 1/******************************************************************************
2 * T H E R M I N A T O R *
3 * THERMal heavy-IoN generATOR *
4 * version 1.0 *
5 * *
6 * Authors of the model: Wojciech Broniowski, Wojciech.Broniowski@ifj.edu.pl, *
7 * Wojciech Florkowski, Wojciech.Florkowski@ifj.edu.pl *
8 * Authors of the code: Adam Kisiel, kisiel@if.pw.edu.pl *
9 * Tomasz Taluc, ttaluc@if.pw.edu.pl *
10 * Code designers: Adam Kisiel, Tomasz Taluc, Wojciech Broniowski, *
11 * Wojciech Florkowski *
12 * *
13 * For the detailed description of the program and furhter references *
14 * to the description of the model plesase refer to: nucl-th/0504047, *
15 * accessibile at: http://www.arxiv.org/nucl-th/0504047 *
16 * *
17 * Homepage: http://hirg.if.pw.edu.pl/en/therminator/ *
18 * *
19 * This code can be freely used and redistributed. However if you decide to *
20 * make modifications to the code, please contact the authors, especially *
21 * if you plan to publish the results obtained with such modified code. *
22 * Any publication of results obtained using this code must include the *
23 * reference to nucl-th/0504047 and the published version of it, when *
24 * available. *
25 * *
26 *****************************************************************************/
27#include "THGlobal.h"
28#include "ParticleDecayer.h"
29#include "DecayChannel.h"
30#include "DecayTable.h"
31#include <TMath.h>
32#include <TDatime.h>
33
34inline Double_t ParticleDecayer::BreitWigner(Double_t Mass, Double_t Gamma) const{
35 Double_t x,y;
36
37 y=mRandom->Rndm();
38 x=Mass+Gamma/2*TMath::Tan(TMath::Pi()*(y-0.5));
39
40 return x;
41}
42
43ParticleDecayer::ParticleDecayer(ParticleDB *aDB)
44{
45 TDatime dat;
46
47 mDB = aDB;
48 mRandom = new TRandom2();
49 mRandom->SetSeed(dat.Get() / 2 * 3);
50}
51
52ParticleDecayer::~ParticleDecayer()
53{
54 delete mRandom;
55}
56
57void
58ParticleDecayer::DecayParticle(Particle *aFather, Particle** aDaughter1, Particle** aDaughter2, Particle** aDaughter3)
59{
60 ParticleType *tType = aFather->GetParticleType();
61#ifdef _RESCALE_CHANNELS_
62 int tChannelIndex = (aFather->GetParticleType())->GetTable()->ChooseDecayChannel(mRandom->Rndm());
63#else
64 Double_t tProb = mRandom->Rndm();
65 int tChannelIndex = (aFather->GetParticleType())->GetTable()->ChooseDecayChannelOrNot(tProb);
66 if (tChannelIndex == -1) {
67 (*aDaughter1) = NULL;
68 (*aDaughter2) = NULL;
69 (*aDaughter3) = NULL;
70
71 DecayTable *tab = (aFather->GetParticleType())->GetTable();
72 PRINT_DEBUG_3("Not decaying " << (aFather->GetParticleType())->GetName() << " for prob " << tProb);
73 for (int tIter=0; tIter<=tab->GetChannelCount(); tIter++) {
74 PRINT_DEBUG_3(mDB->GetParticleType(tab->GetDecayChannel(tIter)->GetParticle1())->GetName() << " ");
75 PRINT_DEBUG_3(mDB->GetParticleType(tab->GetDecayChannel(tIter)->GetParticle2())->GetName() << " ");
76 if (tab->GetDecayChannel(tIter)->GetParticle3()>-1)
77 PRINT_DEBUG_3(mDB->GetParticleType(tab->GetDecayChannel(tIter)->GetParticle3())->GetName() << " ");
78 PRINT_DEBUG_3(tab->GetDecayChannel(tIter)->GetBranchingRatio() << endl);
79 }
80
81 return;
82 }
83
84#endif
85 const DecayChannel *tChan = (aFather->GetParticleType())->GetTable()->GetDecayChannel(tChannelIndex);
86
87 if (tChan->Is3Particle())
88 {
89 // It is a 3-body decay channel
90 ParticleType *tType1 = mDB->GetParticleType(tChan->GetParticle1());
91 ParticleType *tType2 = mDB->GetParticleType(tChan->GetParticle2());
92 ParticleType *tType3 = mDB->GetParticleType(tChan->GetParticle3());
93
94 double tE = aFather->GetEnergy();
95 double tM = aFather->GetMass(); /*MCH uncommented*/
96 // double tFatherMass=aFather->GetMass();
97 // double tFatherGamma=tType->GetGamma();
98 double tM1 = tType1->GetMass();
99 double tM2 = tType2->GetMass();
100 double tM3 = tType3->GetMass();
101/* MCH commented begin
102 double tM;
103 do {
104 tM=BreitWigner(tFatherMass,tFatherGamma);
105 }
106 while (tM1+tM2+tM3>tM);
107MCH commented end*/
108 double tES1, tES2, tP1, tP2, tCos12, tZ;
109
110 do
111 {
112 // Generate E1 and E2 with the Momnte-Carlo method
113 do
114 {
115 tES1 = mRandom->Rndm() * (tM - tM2 - tM3 - tM1) + tM1;
116 tES2 = mRandom->Rndm() * (tM - tM1 - tM3 - tM2) + tM2;
117 }
118 while (tES1+tES2 > tM); // The sum of both energies must be smaller than the resonance mass
119
120 tP1 = TMath::Sqrt(tES1*tES1 - tM1*tM1);
121 tP2 = TMath::Sqrt(tES2*tES2 - tM2*tM2);
122
123 tZ = tM - tES1 - tES2;
124 tZ *= tZ;
125 tCos12 = (tZ - tP1*tP1 - tP2*tP2 - tM3*tM3)/(2*tP1*tP2);
126 }
127 while ((tCos12 < -1.0) || (tCos12 > 1.0)); // Cos Theta must exist (be within -1.0 to 1.0 )
128
129 double tTime;
130 if (tType->GetGamma() == 0.0)
131 tTime = 1.0e10;
132 else {
133 double tTau0 = tE/(tType->GetMass()*tType->GetGamma());
134 // When it decays
135 tTime = -tTau0*TMath::Log(mRandom->Rndm());
136 }
137
138 // Decay coordinates
139 double rxr = aFather->rx + (aFather->px/tE)*tTime;
140 double ryr = aFather->ry + (aFather->py/tE)*tTime;
141 double rzr = aFather->rz + (aFather->pz/tE)*tTime;
142 double rtr = aFather->rt + tTime;
143
144 double tPxr2 = tP2 * TMath::Sqrt(1-tCos12*tCos12);
145 double tPzr2 = tP2*tCos12;
146 double tPxr3 = - tPxr2;
147 double tPzr3 = - (tP1 + tPzr2);
148 double tP3 = TMath::Hypot(tPxr3, tPzr3);
149 double tES3 = TMath::Hypot(tM3, tP3);
150
151 // Generating Euler angles
152 double tPhi = mRandom->Rndm() * 2 * TMath::Pi();
153 double tKsi = mRandom->Rndm() * 2 * TMath::Pi();
154 double tCosTh = mRandom->Rndm() * 2.0 - 1.0;
155
156 double sp = TMath::Sin(tPhi);
157 double cp = TMath::Cos(tPhi);
158 double sk = TMath::Sin(tKsi);
159 double ck = TMath::Cos(tKsi);
160 double st = TMath::Sqrt(1.0-tCosTh*tCosTh);
161 double ct = tCosTh;
162
163 // Rotating the whole system
164 double tPxp1 = - st*ck * tP1;
165 double tPyp1 = st*sk * tP1;
166 double tPzp1 = ct * tP1;
167
168 double tPxp2 = (cp*ct*ck - sp*sk) * tPxr2 + (-st*ck) * tPzr2;
169 double tPyp2 = (-cp*ct*sk - sp*ck) * tPxr2 + (st*sk) * tPzr2;
170 double tPzp2 = cp*st * tPxr2 + ct * tPzr2;
171
172 double tPxp3 = (cp*ct*ck - sp*sk) * tPxr3 + (-st*ck) * tPzr3;
173 double tPyp3 = (-cp*ct*sk - sp*ck) * tPxr3 + (st*sk) * tPzr3;
174 double tPzp3 = cp*st * tPxr3 + ct * tPzr3;
175
176 double tVx = aFather->px/aFather->GetEnergy();
177 double tVy = aFather->py/aFather->GetEnergy();
178 double tVz = aFather->pz/aFather->GetEnergy();
179
180 tES1 = TMath::Sqrt(tM1*tM1+tPxp1*tPxp1+tPyp1*tPyp1+tPzp1*tPzp1);
181 tES2 = TMath::Sqrt(tM2*tM2+tPxp2*tPxp2+tPyp2*tPyp2+tPzp2*tPzp2);
182 tES3 = TMath::Sqrt(tM3*tM3+tPxp3*tPxp3+tPyp3*tPyp3+tPzp3*tPzp3);
183
184 double tV2 = tVx*tVx + tVy*tVy + tVz*tVz;
185 double tGamma = TMath::Power(1-tV2,-0.5);
186
187 // Boosting by the parent velocity
188 double tVP = tVx*tPxp1 + tVy*tPyp1 + tVz*tPzp1;
189 double tgvp = (tGamma - 1.0) * (1.0/tV2) * tVP;
190
191 double tPx1 = tPxp1 + (tgvp + tGamma * tES1) * tVx;
192 double tPy1 = tPyp1 + (tgvp + tGamma * tES1) * tVy;
193 double tPz1 = tPzp1 + (tgvp + tGamma * tES1) * tVz;
194
195 tVP = tVx*tPxp2 + tVy*tPyp2 + tVz*tPzp2;
196 tgvp = (tGamma - 1.0) * (1.0/tV2) * tVP;
197
198 double tPx2 = tPxp2 + (tgvp + tGamma * tES2) * tVx;
199 double tPy2 = tPyp2 + (tgvp + tGamma * tES2) * tVy;
200 double tPz2 = tPzp2 + (tgvp + tGamma * tES2) * tVz;
201
202 tVP = tVx*tPxp3 + tVy*tPyp3 + tVz*tPzp3;
203 tgvp = (tGamma - 1.0) * (1.0/tV2) * tVP;
204
205 double tPx3 = tPxp3 + (tgvp + tGamma * tES3) * tVx;
206 double tPy3 = tPyp3 + (tgvp + tGamma * tES3) * tVy;
207 double tPz3 = tPzp3 + (tgvp + tGamma * tES3) * tVz;
208
209 tES1 = TMath::Sqrt(tM1*tM1+tPx1*tPx1+tPy1*tPy1+tPz1*tPz1);
210 tES2 = TMath::Sqrt(tM2*tM2+tPx2*tPx2+tPy2*tPy2+tPz2*tPz2);
211 tES3 = TMath::Sqrt(tM3*tM3+tPx3*tPx3+tPy3*tPy3+tPz3*tPz3);
212
213 (*aDaughter1) = new Particle(tType1,
214 tPx1, tPy1, tPz1,
215 rxr, ryr, rzr,
216 rtr);
217 (*aDaughter2) = new Particle(tType2,
218 tPx2, tPy2, tPz2,
219 rxr, ryr, rzr,
220 rtr);
221 (*aDaughter3) = new Particle(tType3,
222 tPx3, tPy3, tPz3,
223 rxr, ryr, rzr,
224 rtr);
225
226 aFather->SetDecayed();
227 }
228 else
229 {
230 // It is a regular two-body decay channel
231 ParticleType *tType1 = mDB->GetParticleType(tChan->GetParticle1());
232 ParticleType *tType2 = mDB->GetParticleType(tChan->GetParticle2());
233
234 double tE = aFather->GetEnergy();
235 double tM = aFather->GetMass(); /*MCH uncommented*/
236 // double tFatherMass=aFather->GetMass();
237 // double tFatherGamma=tType->GetGamma();
238 double tM1 = tType1->GetMass();
239 double tM2 = tType2->GetMass();
240/* MCH commented begin
241 double tM;
242 do {
243 tM=BreitWigner(tFatherMass,tFatherGamma);
244 }
245 while (tM1+tM2>tM);
246MCH commented end*/
247
248 double tTime;
249 if (tType->GetGamma() == 0.0)
250 tTime = 1.0e10;
251 else {
252 double tTau0 = tE/(tType->GetMass()*tType->GetGamma());
253 // When it decays
254 tTime = -tTau0*TMath::Log(mRandom->Rndm());
255 }
256
257 // Decay coordinates
258 double rxr = aFather->rx + (aFather->px/tE)*tTime;
259 double ryr = aFather->ry + (aFather->py/tE)*tTime;
260 double rzr = aFather->rz + (aFather->pz/tE)*tTime;
261 double rtr = aFather->rt + tTime;
262
263 // Decay energy
264 double tMC1 = (tM*tM - (tM1+tM2)*(tM1+tM2));
265 double tMC2 = (tM*tM - (tM1-tM2)*(tM1-tM2));
266 double tMom = TMath::Sqrt(tMC1*tMC2)/(2*tM);
267 double tPhi = mRandom->Rndm() * 2 * TMath::Pi();
268 double tCosTh = mRandom->Rndm()*2.0-1.0;
269
270 double tPtr = tMom*TMath::Sqrt(1-tCosTh*tCosTh);
271 double tPxr1 = tPtr*TMath::Cos(tPhi);
272 double tPyr1 = tPtr*TMath::Sin(tPhi);
273 double tPzr1 = tMom*tCosTh;
274
275 double tVx = aFather->px/aFather->GetEnergy();
276 double tVy = aFather->py/aFather->GetEnergy();
277 double tVz = aFather->pz/aFather->GetEnergy();
278
279 double tES1 = TMath::Sqrt(tM1*tM1+tPxr1*tPxr1+tPyr1*tPyr1+tPzr1*tPzr1);
280 double tES2 = TMath::Sqrt(tM2*tM2+tPxr1*tPxr1+tPyr1*tPyr1+tPzr1*tPzr1);
281
282 double tV2 = tVx*tVx + tVy*tVy + tVz*tVz;
283 double tGamma = TMath::Power(1-tV2,-0.5);
284 double tVP = tVx*tPxr1 + tVy*tPyr1 + tVz*tPzr1;
285 double tgvp = (tGamma - 1.0) * (1.0/tV2) * tVP;
286
287 double tPx1 = tPxr1 + (tgvp + tGamma * tES1) * tVx;
288 double tPy1 = tPyr1 + (tgvp + tGamma * tES1) * tVy;
289 double tPz1 = tPzr1 + (tgvp + tGamma * tES1) * tVz;
290
291 double tPx2 = - tPxr1 + (-tgvp + tGamma * tES2) * tVx;
292 double tPy2 = - tPyr1 + (-tgvp + tGamma * tES2) * tVy;
293 double tPz2 = - tPzr1 + (-tgvp + tGamma * tES2) * tVz;
294
295 (*aDaughter1) = new Particle(tType1,
296 tPx1, tPy1, tPz1,
297 rxr, ryr, rzr,
298 rtr);
299 (*aDaughter2) = new Particle(tType2,
300 tPx2, tPy2, tPz2,
301 rxr, ryr, rzr,
302 rtr);
303 (*aDaughter3) = 0;
304
305 aFather->SetDecayed();
306
307 }
308}
309
310void ParticleDecayer::SeedSet(int aSeed)
311{
312 // mRandom->SetSeed2(aSeed, aSeed * 11 % 9);
313 mRandom->SetSeed(aSeed);
314}
315