Update master to aliroot
[u/mrichter/AliRoot.git] / TTherminator / Therminator / ParticleDecayer.cxx
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
34 inline 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
43 ParticleDecayer::ParticleDecayer(ParticleDB *aDB)
44 {
45   TDatime dat;
46
47   mDB = aDB;
48   mRandom = new TRandom2();
49   mRandom->SetSeed(dat.Get() / 2 * 3);
50 }
51
52 ParticleDecayer::~ParticleDecayer()
53 {
54   delete mRandom;
55 }
56
57 void
58 ParticleDecayer::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);
107 MCH 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);
246 MCH 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
310 void ParticleDecayer::SeedSet(int aSeed)
311 {
312   //  mRandom->SetSeed2(aSeed, aSeed * 11 % 9);
313   mRandom->SetSeed(aSeed);
314 }
315