[u/mrichter/AliRoot.git] / TTherminator / Therminator / ParticleDecayer.cxx

/******************************************************************************
 *                      T H E R M I N A T O R                                 *
 *                   THERMal heavy-IoN generATOR                              *
 *                           version 1.0                                      *
 *                                                                            *
 * Authors of the model: Wojciech Broniowski, Wojciech.Broniowski@ifj.edu.pl, *
 *                       Wojciech Florkowski, Wojciech.Florkowski@ifj.edu.pl  *
 * Authors of the code:  Adam Kisiel, kisiel@if.pw.edu.pl                     *
 *                       Tomasz Taluc, ttaluc@if.pw.edu.pl                    *
 * Code designers: Adam Kisiel, Tomasz Taluc, Wojciech Broniowski,            *
 *                 Wojciech Florkowski                                        *
 *                                                                            *
 * For the detailed description of the program and furhter references         * 
 * to the description of the model plesase refer to: nucl-th/0504047,         *
 * accessibile at: http://www.arxiv.org/nucl-th/0504047                       *
 *                                                                            *
 * Homepage: http://hirg.if.pw.edu.pl/en/therminator/                         *
 *                                                                            *
 * This code can be freely used and redistributed. However if you decide to   *
 * make modifications to the code, please contact the authors, especially     *
 * if you plan to publish the results obtained with such modified code.       *
 * Any publication of results obtained using this code must include the       *
 * reference to nucl-th/0504047 and the published version of it, when         *
 * available.                                                                 *
 *                                                                            *
 *****************************************************************************/
#include "THGlobal.h"
#include "ParticleDecayer.h"
#include "DecayChannel.h"
#include "DecayTable.h"
#include <TMath.h>
#include <TDatime.h>

inline Double_t ParticleDecayer::BreitWigner(Double_t Mass, Double_t Gamma) const{
  Double_t x,y;

  y=mRandom->Rndm();
  x=Mass+Gamma/2*TMath::Tan(TMath::Pi()*(y-0.5));

  return x;
}

ParticleDecayer::ParticleDecayer(ParticleDB *aDB)
{
  TDatime dat;

  mDB = aDB;
  mRandom = new TRandom2();
  mRandom->SetSeed(dat.Get() / 2 * 3);
}

ParticleDecayer::~ParticleDecayer()
{
  delete mRandom;
}

void
ParticleDecayer::DecayParticle(Particle *aFather, Particle** aDaughter1, Particle** aDaughter2, Particle** aDaughter3)
{
  ParticleType *tType = aFather->GetParticleType();
#ifdef _RESCALE_CHANNELS_
  int tChannelIndex = (aFather->GetParticleType())->GetTable()->ChooseDecayChannel(mRandom->Rndm());
#else
  Double_t tProb = mRandom->Rndm();
  int tChannelIndex = (aFather->GetParticleType())->GetTable()->ChooseDecayChannelOrNot(tProb);
  if (tChannelIndex == -1) {
    (*aDaughter1) =  NULL;
    (*aDaughter2) =  NULL;
    (*aDaughter3) =  NULL; 
    
    DecayTable *tab = (aFather->GetParticleType())->GetTable();
    PRINT_DEBUG_3("Not decaying " << (aFather->GetParticleType())->GetName() << " for prob " << tProb);
    for (int tIter=0; tIter<=tab->GetChannelCount(); tIter++) {
      PRINT_DEBUG_3(mDB->GetParticleType(tab->GetDecayChannel(tIter)->GetParticle1())->GetName() << " ");
      PRINT_DEBUG_3(mDB->GetParticleType(tab->GetDecayChannel(tIter)->GetParticle2())->GetName() << " ");
      if (tab->GetDecayChannel(tIter)->GetParticle3()>-1)
	PRINT_DEBUG_3(mDB->GetParticleType(tab->GetDecayChannel(tIter)->GetParticle3())->GetName() << " ");
      PRINT_DEBUG_3(tab->GetDecayChannel(tIter)->GetBranchingRatio() << endl);
    }
    
    return;
  }

#endif
  const DecayChannel *tChan = (aFather->GetParticleType())->GetTable()->GetDecayChannel(tChannelIndex);

  if (tChan->Is3Particle())
    {
      // It is a 3-body decay channel
      ParticleType *tType1 = mDB->GetParticleType(tChan->GetParticle1());
      ParticleType *tType2 = mDB->GetParticleType(tChan->GetParticle2());
      ParticleType *tType3 = mDB->GetParticleType(tChan->GetParticle3());
      
      double tE = aFather->GetEnergy();
      double tM = aFather->GetMass();		/*MCH uncommented*/
      //      double tFatherMass=aFather->GetMass();
      //      double tFatherGamma=tType->GetGamma();
      double tM1 = tType1->GetMass();
      double tM2 = tType2->GetMass();
      double tM3 = tType3->GetMass();
/* MCH commented begin
      double tM;
      do {
	tM=BreitWigner(tFatherMass,tFatherGamma);
      }
      while (tM1+tM2+tM3>tM);
MCH commented end*/
      double tES1, tES2, tP1, tP2, tCos12, tZ;
      
      do 
	{
	  // Generate E1 and E2 with the Momnte-Carlo method
	  do 
	    {
	      tES1 = mRandom->Rndm() * (tM - tM2 - tM3 - tM1) + tM1;
	      tES2 = mRandom->Rndm() * (tM - tM1 - tM3 - tM2) + tM2;
	    }
	  while (tES1+tES2 > tM); // The sum of both energies must be smaller than the resonance mass
	  
	  tP1  = TMath::Sqrt(tES1*tES1 - tM1*tM1);
	  tP2  = TMath::Sqrt(tES2*tES2 - tM2*tM2);
	  
	  tZ = tM - tES1 - tES2;
	  tZ *= tZ;
	  tCos12 = (tZ - tP1*tP1 - tP2*tP2 - tM3*tM3)/(2*tP1*tP2);
	}
      while ((tCos12 < -1.0) || (tCos12 > 1.0)); // Cos Theta must exist (be within -1.0 to 1.0 )
      
      double tTime;
      if (tType->GetGamma() == 0.0)
	tTime = 1.0e10;
      else {
	double tTau0 = tE/(tType->GetMass()*tType->GetGamma());
	// When it decays
	tTime = -tTau0*TMath::Log(mRandom->Rndm());
      }
      
      // Decay coordinates
      double rxr = aFather->rx + (aFather->px/tE)*tTime;
      double ryr = aFather->ry + (aFather->py/tE)*tTime;
      double rzr = aFather->rz + (aFather->pz/tE)*tTime;
      double rtr = aFather->rt + tTime;
      
      double tPxr2 = tP2 * TMath::Sqrt(1-tCos12*tCos12);
      double tPzr2 = tP2*tCos12;
      double tPxr3 = - tPxr2;
      double tPzr3 = - (tP1 + tPzr2);
      double tP3 = TMath::Hypot(tPxr3, tPzr3);
      double tES3 = TMath::Hypot(tM3, tP3);

      // Generating Euler angles
      double tPhi = mRandom->Rndm() * 2 * TMath::Pi();
      double tKsi = mRandom->Rndm() * 2 * TMath::Pi();
      double tCosTh = mRandom->Rndm() * 2.0 - 1.0;

      double sp = TMath::Sin(tPhi);
      double cp = TMath::Cos(tPhi);
      double sk = TMath::Sin(tKsi);
      double ck = TMath::Cos(tKsi);
      double st = TMath::Sqrt(1.0-tCosTh*tCosTh);
      double ct = tCosTh;

      // Rotating the whole system
      double tPxp1 = - st*ck * tP1;
      double tPyp1 = st*sk * tP1;
      double tPzp1 = ct * tP1;
      
      double tPxp2 = (cp*ct*ck - sp*sk)  * tPxr2 + (-st*ck) * tPzr2;
      double tPyp2 = (-cp*ct*sk - sp*ck) * tPxr2 + (st*sk)  * tPzr2;
      double tPzp2 = cp*st               * tPxr2 + ct       * tPzr2;

      double tPxp3 = (cp*ct*ck - sp*sk)  * tPxr3 + (-st*ck) * tPzr3;
      double tPyp3 = (-cp*ct*sk - sp*ck) * tPxr3 + (st*sk)  * tPzr3;
      double tPzp3 = cp*st               * tPxr3 + ct       * tPzr3;
      
      double tVx = aFather->px/aFather->GetEnergy();
      double tVy = aFather->py/aFather->GetEnergy();
      double tVz = aFather->pz/aFather->GetEnergy();
      
       tES1 = TMath::Sqrt(tM1*tM1+tPxp1*tPxp1+tPyp1*tPyp1+tPzp1*tPzp1);
       tES2 = TMath::Sqrt(tM2*tM2+tPxp2*tPxp2+tPyp2*tPyp2+tPzp2*tPzp2);
       tES3 = TMath::Sqrt(tM3*tM3+tPxp3*tPxp3+tPyp3*tPyp3+tPzp3*tPzp3);

      double tV2 = tVx*tVx + tVy*tVy + tVz*tVz;
      double tGamma = TMath::Power(1-tV2,-0.5);
      
      // Boosting by the parent velocity
      double tVP = tVx*tPxp1 + tVy*tPyp1 + tVz*tPzp1;
      double tgvp = (tGamma - 1.0) * (1.0/tV2) * tVP;

      double tPx1 =   tPxp1 + (tgvp + tGamma * tES1) * tVx;
      double tPy1 =   tPyp1 + (tgvp + tGamma * tES1) * tVy;
      double tPz1 =   tPzp1 + (tgvp + tGamma * tES1) * tVz;
  
      tVP = tVx*tPxp2 + tVy*tPyp2 + tVz*tPzp2;
      tgvp = (tGamma - 1.0) * (1.0/tV2) * tVP;

      double tPx2 =   tPxp2 + (tgvp + tGamma * tES2) * tVx;
      double tPy2 =   tPyp2 + (tgvp + tGamma * tES2) * tVy;
      double tPz2 =   tPzp2 + (tgvp + tGamma * tES2) * tVz;
  
      tVP = tVx*tPxp3 + tVy*tPyp3 + tVz*tPzp3;
      tgvp = (tGamma - 1.0) * (1.0/tV2) * tVP;

      double tPx3 =   tPxp3 + (tgvp + tGamma * tES3) * tVx;
      double tPy3 =   tPyp3 + (tgvp + tGamma * tES3) * tVy;
      double tPz3 =   tPzp3 + (tgvp + tGamma * tES3) * tVz;
      
      tES1 = TMath::Sqrt(tM1*tM1+tPx1*tPx1+tPy1*tPy1+tPz1*tPz1);
      tES2 = TMath::Sqrt(tM2*tM2+tPx2*tPx2+tPy2*tPy2+tPz2*tPz2);
      tES3 = TMath::Sqrt(tM3*tM3+tPx3*tPx3+tPy3*tPy3+tPz3*tPz3);
      
      (*aDaughter1) = new Particle(tType1, 
				   tPx1, tPy1, tPz1,
				   rxr, ryr, rzr, 
				   rtr);
      (*aDaughter2) = new Particle(tType2, 
				   tPx2, tPy2, tPz2,
				   rxr, ryr, rzr, 
				   rtr);
      (*aDaughter3) = new Particle(tType3, 
				   tPx3, tPy3, tPz3,
				   rxr, ryr, rzr, 
				   rtr);

      aFather->SetDecayed();
    }
  else
    {
      // It is a regular two-body decay channel
      ParticleType *tType1 = mDB->GetParticleType(tChan->GetParticle1());
      ParticleType *tType2 = mDB->GetParticleType(tChan->GetParticle2());
      
      double tE = aFather->GetEnergy();
      double tM = aFather->GetMass();		/*MCH uncommented*/
      //      double tFatherMass=aFather->GetMass();
      //      double tFatherGamma=tType->GetGamma();
      double tM1 = tType1->GetMass();
      double tM2 = tType2->GetMass();
/* MCH commented begin
      double tM;
      do {
	tM=BreitWigner(tFatherMass,tFatherGamma);
      }
      while (tM1+tM2>tM);
MCH commented end*/

      double tTime;
      if (tType->GetGamma() == 0.0)
	tTime = 1.0e10;
      else {
	double tTau0 = tE/(tType->GetMass()*tType->GetGamma());
	// When it decays
	tTime = -tTau0*TMath::Log(mRandom->Rndm());
      }
      
      // Decay coordinates
      double rxr = aFather->rx + (aFather->px/tE)*tTime;
      double ryr = aFather->ry + (aFather->py/tE)*tTime;
      double rzr = aFather->rz + (aFather->pz/tE)*tTime;
      double rtr = aFather->rt + tTime;
      
      // Decay energy
      double tMC1 = (tM*tM - (tM1+tM2)*(tM1+tM2));
      double tMC2 = (tM*tM - (tM1-tM2)*(tM1-tM2));
      double tMom = TMath::Sqrt(tMC1*tMC2)/(2*tM);
      double tPhi = mRandom->Rndm() * 2 * TMath::Pi();
      double tCosTh = mRandom->Rndm()*2.0-1.0;
      
      double tPtr = tMom*TMath::Sqrt(1-tCosTh*tCosTh);
      double tPxr1 = tPtr*TMath::Cos(tPhi);
      double tPyr1 = tPtr*TMath::Sin(tPhi);
      double tPzr1 = tMom*tCosTh;
      
      double tVx = aFather->px/aFather->GetEnergy();
      double tVy = aFather->py/aFather->GetEnergy();
      double tVz = aFather->pz/aFather->GetEnergy();
      
      double tES1 = TMath::Sqrt(tM1*tM1+tPxr1*tPxr1+tPyr1*tPyr1+tPzr1*tPzr1);
      double tES2 = TMath::Sqrt(tM2*tM2+tPxr1*tPxr1+tPyr1*tPyr1+tPzr1*tPzr1);
      
      double tV2 = tVx*tVx + tVy*tVy + tVz*tVz;
      double tGamma = TMath::Power(1-tV2,-0.5);
      double tVP = tVx*tPxr1 + tVy*tPyr1 + tVz*tPzr1;
      double tgvp = (tGamma - 1.0) * (1.0/tV2) * tVP;
      
      double tPx1 =   tPxr1 + (tgvp + tGamma * tES1) * tVx;
      double tPy1 =   tPyr1 + (tgvp + tGamma * tES1) * tVy;
      double tPz1 =   tPzr1 + (tgvp + tGamma * tES1) * tVz;
  
      double tPx2 = - tPxr1 + (-tgvp + tGamma * tES2) * tVx;
      double tPy2 = - tPyr1 + (-tgvp + tGamma * tES2) * tVy;
      double tPz2 = - tPzr1 + (-tgvp + tGamma * tES2) * tVz;

      (*aDaughter1) = new Particle(tType1, 
				   tPx1, tPy1, tPz1,
				   rxr, ryr, rzr, 
				   rtr);
      (*aDaughter2) = new Particle(tType2, 
				   tPx2, tPy2, tPz2,
				   rxr, ryr, rzr, 
				   rtr);
      (*aDaughter3) = 0;

      aFather->SetDecayed();
  
    }
}

void ParticleDecayer::SeedSet(int aSeed)
{
  //  mRandom->SetSeed2(aSeed, aSeed * 11 % 9);
  mRandom->SetSeed(aSeed);
}
Commit	Line	Data
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
	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/2TMath::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(tES1tES1 - tM1tM1);
121	tP2 = TMath::Sqrt(tES2tES2 - tM2tM2);
122
123	tZ = tM - tES1 - tES2;
124	tZ *= tZ;
125	tCos12 = (tZ - tP1tP1 - tP2tP2 - tM3tM3)/(2tP1*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 = - stck tP1;
165	double tPyp1 = stsk tP1;
166	double tPzp1 = ct * tP1;
167
168	double tPxp2 = (cpctck - spsk) tPxr2 + (-stck) tPzr2;
169	double tPyp2 = (-cpctsk - spck) tPxr2 + (stsk) tPzr2;
170	double tPzp2 = cpst tPxr2 + ct * tPzr2;
171
172	double tPxp3 = (cpctck - spsk) tPxr3 + (-stck) tPzr3;
173	double tPyp3 = (-cpctsk - spck) tPxr3 + (stsk) tPzr3;
174	double tPzp3 = cpst 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(tM1tM1+tPxp1tPxp1+tPyp1tPyp1+tPzp1tPzp1);
181	tES2 = TMath::Sqrt(tM2tM2+tPxp2tPxp2+tPyp2tPyp2+tPzp2tPzp2);
182	tES3 = TMath::Sqrt(tM3tM3+tPxp3tPxp3+tPyp3tPyp3+tPzp3tPzp3);
183
184	double tV2 = tVxtVx + tVytVy + tVz*tVz;
185	double tGamma = TMath::Power(1-tV2,-0.5);
186
187	// Boosting by the parent velocity
188	double tVP = tVxtPxp1 + tVytPyp1 + 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 = tVxtPxp2 + tVytPyp2 + 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 = tVxtPxp3 + tVytPyp3 + 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(tM1tM1+tPx1tPx1+tPy1tPy1+tPz1tPz1);
210	tES2 = TMath::Sqrt(tM2tM2+tPx2tPx2+tPy2tPy2+tPz2tPz2);
211	tES3 = TMath::Sqrt(tM3tM3+tPx3tPx3+tPy3tPy3+tPz3tPz3);
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 = (tMtM - (tM1+tM2)(tM1+tM2));
265	double tMC2 = (tMtM - (tM1-tM2)(tM1-tM2));
266	double tMom = TMath::Sqrt(tMC1tMC2)/(2tM);
267	double tPhi = mRandom->Rndm() * 2 * TMath::Pi();
268	double tCosTh = mRandom->Rndm()*2.0-1.0;
269
270	double tPtr = tMomTMath::Sqrt(1-tCosThtCosTh);
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(tM1tM1+tPxr1tPxr1+tPyr1tPyr1+tPzr1tPzr1);
280	double tES2 = TMath::Sqrt(tM2tM2+tPxr1tPxr1+tPyr1tPyr1+tPzr1tPzr1);
281
282	double tV2 = tVxtVx + tVytVy + tVz*tVz;
283	double tGamma = TMath::Power(1-tV2,-0.5);
284	double tVP = tVxtPxr1 + tVytPyr1 + 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