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