#include "TGenPhaseSpace.h"
#include "TRandom.h"
#include "TMath.h"
#include <stdlib.h>
const Int_t kMAXP = 18;
ClassImp(TGenPhaseSpace)
Double_t TGenPhaseSpace::PDK(Double_t a, Double_t b, Double_t c)
{
Double_t x = (a-b-c)*(a+b+c)*(a-b+c)*(a+b-c);
x = TMath::Sqrt(x)/(2*a);
return x;
}
Int_t DoubleMax(const void *a, const void *b)
{
Double_t aa = * ((Double_t *) a);
Double_t bb = * ((Double_t *) b);
if (aa > bb) return 1;
if (aa < bb) return -1;
return 0;
}
TGenPhaseSpace::TGenPhaseSpace(const TGenPhaseSpace &gen) : TObject(gen)
{
fNt = gen.fNt;
fWtMax = gen.fWtMax;
fTeCmTm = gen.fTeCmTm;
fBeta[0] = gen.fBeta[0];
fBeta[1] = gen.fBeta[1];
fBeta[2] = gen.fBeta[2];
for (Int_t i=0;i<fNt;i++) {
fMass[i] = gen.fMass[i];
fDecPro[i] = gen.fDecPro[i];
}
}
Double_t TGenPhaseSpace::Generate()
{
Double_t rno[kMAXP];
rno[0] = 0;
Int_t n;
if (fNt>2) {
for (n=1; n<fNt-1; n++) rno[n]=gRandom->Rndm();
qsort(rno+1 ,fNt-2 ,sizeof(Double_t) ,DoubleMax);
}
rno[fNt-1] = 1;
Double_t invMas[kMAXP], sum=0;
for (n=0; n<fNt; n++) {
sum += fMass[n];
invMas[n] = rno[n]*fTeCmTm + sum;
}
Double_t wt=fWtMax;
Double_t pd[kMAXP];
for (n=0; n<fNt-1; n++) {
pd[n] = PDK(invMas[n+1],invMas[n],fMass[n+1]);
wt *= pd[n];
}
fDecPro[0].SetPxPyPzE(0, pd[0], 0 , TMath::Sqrt(pd[0]*pd[0]+fMass[0]*fMass[0]) );
Int_t i=1;
Int_t j;
while (1) {
fDecPro[i].SetPxPyPzE(0, -pd[i-1], 0 , TMath::Sqrt(pd[i-1]*pd[i-1]+fMass[i]*fMass[i]) );
Double_t cZ = 2*gRandom->Rndm() - 1;
Double_t sZ = TMath::Sqrt(1-cZ*cZ);
Double_t angY = 2*TMath::Pi() * gRandom->Rndm();
Double_t cY = TMath::Cos(angY);
Double_t sY = TMath::Sin(angY);
for (j=0; j<=i; j++) {
TLorentzVector *v = fDecPro+j;
Double_t x = v->Px();
Double_t y = v->Py();
v->SetPx( cZ*x - sZ*y );
v->SetPy( sZ*x + cZ*y );
x = v->Px();
Double_t z = v->Pz();
v->SetPx( cY*x - sY*z );
v->SetPz( sY*x + cY*z );
}
if (i == (fNt-1)) break;
Double_t beta = pd[i] / sqrt(pd[i]*pd[i] + invMas[i]*invMas[i]);
for (j=0; j<=i; j++) fDecPro[j].Boost(0,beta,0);
i++;
}
for (n=0;n<fNt;n++) fDecPro[n].Boost(fBeta[0],fBeta[1],fBeta[2]);
return wt;
}
TLorentzVector *TGenPhaseSpace::GetDecay(Int_t n)
{
if (n>fNt) return 0;
return fDecPro+n;
}
Bool_t TGenPhaseSpace::SetDecay(TLorentzVector &P, Int_t nt,
Double_t *mass, Option_t *opt)
{
Int_t n;
fNt = nt;
if (fNt<2 || fNt>18) return kFALSE;
fTeCmTm = P.Mag();
for (n=0;n<fNt;n++) {
fMass[n] = mass[n];
fTeCmTm -= mass[n];
}
if (fTeCmTm<=0) return kFALSE;
if (strcasecmp(opt,"fermi")==0) {
Double_t ffq[] = {0
,3.141592, 19.73921, 62.01255, 129.8788, 204.0131
,256.3704, 268.4705, 240.9780, 189.2637
,132.1308, 83.0202, 47.4210, 24.8295
,12.0006, 5.3858, 2.2560, 0.8859 };
fWtMax = TMath::Power(fTeCmTm,fNt-2) * ffq[fNt-1] / P.Mag();
} else {
Double_t emmax = fTeCmTm + fMass[0];
Double_t emmin = 0;
Double_t wtmax = 1;
for (n=1; n<fNt; n++) {
emmin += fMass[n-1];
emmax += fMass[n];
wtmax *= PDK(emmax, emmin, fMass[n]);
}
fWtMax = 1/wtmax;
}
if (P.Beta()) {
Double_t w = P.Beta()/P.Rho();
fBeta[0] = P(0)*w;
fBeta[1] = P(1)*w;
fBeta[2] = P(2)*w;
}
else fBeta[0]=fBeta[1]=fBeta[2]=0;
return kTRUE;
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.