```// @(#)root/physics:\$Id\$
// Author: Rene Brun , Valerio Filippini  06/09/2000

//_____________________________________________________________________________________
//
//  Utility class to generate n-body event,
//  with constant cross-section (default)
//  or with Fermi energy dependence (opt="Fermi").
//  The event is generated in the center-of-mass frame,
//  but the decay products are finally boosted
//  using the betas of the original particle.
//
//  The code is based on the GENBOD function (W515 from CERNLIB)
//  using the Raubold and Lynch method
//      F. James, Monte Carlo Phase Space, CERN 68-15 (1968)
//
// see example of use in \$ROOTSYS/tutorials/physics/PhaseSpace.C
//
// Note that Momentum, Energy units are Gev/C, GeV

#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)
{
//the PDK function
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)
{
//special max function
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)
{
//copy constructor
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];
}
}

//__________________________________________________________________________________________________
TGenPhaseSpace& TGenPhaseSpace::operator=(const TGenPhaseSpace &gen)
{
// Assignment operator
TObject::operator=(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];
}
return *this;
}

//__________________________________________________________________________________________________
Double_t TGenPhaseSpace::Generate()
{
//  Generate a random final state.
//  The function returns the weigth of the current event.
//  The TLorentzVector of each decay product can be obtained using GetDecay(n).
//
// Note that Momentum, Energy units are Gev/C, GeV

Double_t rno[kMAXP];
rno[0] = 0;
Int_t n;
if (fNt>2) {
for (n=1; n<fNt-1; n++)  rno[n]=gRandom->Rndm();   // fNt-2 random numbers
qsort(rno+1 ,fNt-2 ,sizeof(Double_t) ,DoubleMax);  // sort them
}
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;
}

//
//-----> compute the weight of the current event
//
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];
}

//
//-----> complete specification of event (Raubold-Lynch method)
//
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 );   // rotation around Z
x = v->Px();
Double_t z = v->Pz();
v->SetPx( cY*x - sY*z );
v->SetPz( sY*x + cY*z );   // rotation around Y
}

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++;
}

//
//---> final boost of all particles
//
for (n=0;n<fNt;n++) fDecPro[n].Boost(fBeta[0],fBeta[1],fBeta[2]);

//
//---> return the weigth of event
//
return wt;
}

//__________________________________________________________________________________
TLorentzVector *TGenPhaseSpace::GetDecay(Int_t n)
{
//return Lorentz vector corresponding to decay n
if (n>fNt) return 0;
return fDecPro+n;
}

//_____________________________________________________________________________________
Bool_t TGenPhaseSpace::SetDecay(TLorentzVector &P, Int_t nt,
const Double_t *mass, Option_t *opt)
{
// input:
// TLorentzVector &P:    decay particle (Momentum, Energy units are Gev/C, GeV)
// Int_t nt:             number of decay products
// Double_t *mass:       array of decay product masses
// Option_t *opt:        default -> constant cross section
//                       "Fermi" -> Fermi energy dependece
// return value:
// kTRUE:      the decay is permitted by kinematics
// kFALSE:     the decay is forbidden by kinematics
//

Int_t n;
fNt = nt;
if (fNt<2 || fNt>18) return kFALSE;  // no more then 18 particle

//
//
//
fTeCmTm = P.Mag();           // total energy in C.M. minus the sum of the masses
for (n=0;n<fNt;n++) {
fMass[n]  = mass[n];
fTeCmTm  -= mass[n];
}

if (fTeCmTm<=0) return kFALSE;    // not enough energy for this decay

//
//------> the max weigth depends on opt:
//   opt == "Fermi"  --> fermi energy dependence for cross section
//   else            --> constant cross section as function of TECM (default)
//
if (strcasecmp(opt,"fermi")==0) {
// ffq[] = pi * (2*pi)**(FNt-2) / (FNt-2)!
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;
}

//
//---->  save the betas of the decaying particle
//
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;
}
```
TGenPhaseSpace.cxx:1
TGenPhaseSpace.cxx:2
TGenPhaseSpace.cxx:3
TGenPhaseSpace.cxx:4
TGenPhaseSpace.cxx:5
TGenPhaseSpace.cxx:6
TGenPhaseSpace.cxx:7
TGenPhaseSpace.cxx:8
TGenPhaseSpace.cxx:9
TGenPhaseSpace.cxx:10
TGenPhaseSpace.cxx:11
TGenPhaseSpace.cxx:12
TGenPhaseSpace.cxx:13
TGenPhaseSpace.cxx:14
TGenPhaseSpace.cxx:15
TGenPhaseSpace.cxx:16
TGenPhaseSpace.cxx:17
TGenPhaseSpace.cxx:18
TGenPhaseSpace.cxx:19
TGenPhaseSpace.cxx:20
TGenPhaseSpace.cxx:21
TGenPhaseSpace.cxx:22
TGenPhaseSpace.cxx:23
TGenPhaseSpace.cxx:24
TGenPhaseSpace.cxx:25
TGenPhaseSpace.cxx:26
TGenPhaseSpace.cxx:27
TGenPhaseSpace.cxx:28
TGenPhaseSpace.cxx:29
TGenPhaseSpace.cxx:30
TGenPhaseSpace.cxx:31
TGenPhaseSpace.cxx:32
TGenPhaseSpace.cxx:33
TGenPhaseSpace.cxx:34
TGenPhaseSpace.cxx:35
TGenPhaseSpace.cxx:36
TGenPhaseSpace.cxx:37
TGenPhaseSpace.cxx:38
TGenPhaseSpace.cxx:39
TGenPhaseSpace.cxx:40
TGenPhaseSpace.cxx:41
TGenPhaseSpace.cxx:42
TGenPhaseSpace.cxx:43
TGenPhaseSpace.cxx:44
TGenPhaseSpace.cxx:45
TGenPhaseSpace.cxx:46
TGenPhaseSpace.cxx:47
TGenPhaseSpace.cxx:48
TGenPhaseSpace.cxx:49
TGenPhaseSpace.cxx:50
TGenPhaseSpace.cxx:51
TGenPhaseSpace.cxx:52
TGenPhaseSpace.cxx:53
TGenPhaseSpace.cxx:54
TGenPhaseSpace.cxx:55
TGenPhaseSpace.cxx:56
TGenPhaseSpace.cxx:57
TGenPhaseSpace.cxx:58
TGenPhaseSpace.cxx:59
TGenPhaseSpace.cxx:60
TGenPhaseSpace.cxx:61
TGenPhaseSpace.cxx:62
TGenPhaseSpace.cxx:63
TGenPhaseSpace.cxx:64
TGenPhaseSpace.cxx:65
TGenPhaseSpace.cxx:66
TGenPhaseSpace.cxx:67
TGenPhaseSpace.cxx:68
TGenPhaseSpace.cxx:69
TGenPhaseSpace.cxx:70
TGenPhaseSpace.cxx:71
TGenPhaseSpace.cxx:72
TGenPhaseSpace.cxx:73
TGenPhaseSpace.cxx:74
TGenPhaseSpace.cxx:75
TGenPhaseSpace.cxx:76
TGenPhaseSpace.cxx:77
TGenPhaseSpace.cxx:78
TGenPhaseSpace.cxx:79
TGenPhaseSpace.cxx:80
TGenPhaseSpace.cxx:81
TGenPhaseSpace.cxx:82
TGenPhaseSpace.cxx:83
TGenPhaseSpace.cxx:84
TGenPhaseSpace.cxx:85
TGenPhaseSpace.cxx:86
TGenPhaseSpace.cxx:87
TGenPhaseSpace.cxx:88
TGenPhaseSpace.cxx:89
TGenPhaseSpace.cxx:90
TGenPhaseSpace.cxx:91
TGenPhaseSpace.cxx:92
TGenPhaseSpace.cxx:93
TGenPhaseSpace.cxx:94
TGenPhaseSpace.cxx:95
TGenPhaseSpace.cxx:96
TGenPhaseSpace.cxx:97
TGenPhaseSpace.cxx:98
TGenPhaseSpace.cxx:99
TGenPhaseSpace.cxx:100
TGenPhaseSpace.cxx:101
TGenPhaseSpace.cxx:102
TGenPhaseSpace.cxx:103
TGenPhaseSpace.cxx:104
TGenPhaseSpace.cxx:105
TGenPhaseSpace.cxx:106
TGenPhaseSpace.cxx:107
TGenPhaseSpace.cxx:108
TGenPhaseSpace.cxx:109
TGenPhaseSpace.cxx:110
TGenPhaseSpace.cxx:111
TGenPhaseSpace.cxx:112
TGenPhaseSpace.cxx:113
TGenPhaseSpace.cxx:114
TGenPhaseSpace.cxx:115
TGenPhaseSpace.cxx:116
TGenPhaseSpace.cxx:117
TGenPhaseSpace.cxx:118
TGenPhaseSpace.cxx:119
TGenPhaseSpace.cxx:120
TGenPhaseSpace.cxx:121
TGenPhaseSpace.cxx:122
TGenPhaseSpace.cxx:123
TGenPhaseSpace.cxx:124
TGenPhaseSpace.cxx:125
TGenPhaseSpace.cxx:126
TGenPhaseSpace.cxx:127
TGenPhaseSpace.cxx:128
TGenPhaseSpace.cxx:129
TGenPhaseSpace.cxx:130
TGenPhaseSpace.cxx:131
TGenPhaseSpace.cxx:132
TGenPhaseSpace.cxx:133
TGenPhaseSpace.cxx:134
TGenPhaseSpace.cxx:135
TGenPhaseSpace.cxx:136
TGenPhaseSpace.cxx:137
TGenPhaseSpace.cxx:138
TGenPhaseSpace.cxx:139
TGenPhaseSpace.cxx:140
TGenPhaseSpace.cxx:141
TGenPhaseSpace.cxx:142
TGenPhaseSpace.cxx:143
TGenPhaseSpace.cxx:144
TGenPhaseSpace.cxx:145
TGenPhaseSpace.cxx:146
TGenPhaseSpace.cxx:147
TGenPhaseSpace.cxx:148
TGenPhaseSpace.cxx:149
TGenPhaseSpace.cxx:150
TGenPhaseSpace.cxx:151
TGenPhaseSpace.cxx:152
TGenPhaseSpace.cxx:153
TGenPhaseSpace.cxx:154
TGenPhaseSpace.cxx:155
TGenPhaseSpace.cxx:156
TGenPhaseSpace.cxx:157
TGenPhaseSpace.cxx:158
TGenPhaseSpace.cxx:159
TGenPhaseSpace.cxx:160
TGenPhaseSpace.cxx:161
TGenPhaseSpace.cxx:162
TGenPhaseSpace.cxx:163
TGenPhaseSpace.cxx:164
TGenPhaseSpace.cxx:165
TGenPhaseSpace.cxx:166
TGenPhaseSpace.cxx:167
TGenPhaseSpace.cxx:168
TGenPhaseSpace.cxx:169
TGenPhaseSpace.cxx:170
TGenPhaseSpace.cxx:171
TGenPhaseSpace.cxx:172
TGenPhaseSpace.cxx:173
TGenPhaseSpace.cxx:174
TGenPhaseSpace.cxx:175
TGenPhaseSpace.cxx:176
TGenPhaseSpace.cxx:177
TGenPhaseSpace.cxx:178
TGenPhaseSpace.cxx:179
TGenPhaseSpace.cxx:180
TGenPhaseSpace.cxx:181
TGenPhaseSpace.cxx:182
TGenPhaseSpace.cxx:183
TGenPhaseSpace.cxx:184
TGenPhaseSpace.cxx:185
TGenPhaseSpace.cxx:186
TGenPhaseSpace.cxx:187
TGenPhaseSpace.cxx:188
TGenPhaseSpace.cxx:189
TGenPhaseSpace.cxx:190
TGenPhaseSpace.cxx:191
TGenPhaseSpace.cxx:192
TGenPhaseSpace.cxx:193
TGenPhaseSpace.cxx:194
TGenPhaseSpace.cxx:195
TGenPhaseSpace.cxx:196
TGenPhaseSpace.cxx:197
TGenPhaseSpace.cxx:198
TGenPhaseSpace.cxx:199
TGenPhaseSpace.cxx:200
TGenPhaseSpace.cxx:201
TGenPhaseSpace.cxx:202
TGenPhaseSpace.cxx:203
TGenPhaseSpace.cxx:204
TGenPhaseSpace.cxx:205
TGenPhaseSpace.cxx:206
TGenPhaseSpace.cxx:207
TGenPhaseSpace.cxx:208
TGenPhaseSpace.cxx:209
TGenPhaseSpace.cxx:210
TGenPhaseSpace.cxx:211
TGenPhaseSpace.cxx:212
TGenPhaseSpace.cxx:213
TGenPhaseSpace.cxx:214
TGenPhaseSpace.cxx:215
TGenPhaseSpace.cxx:216
TGenPhaseSpace.cxx:217
TGenPhaseSpace.cxx:218
TGenPhaseSpace.cxx:219
TGenPhaseSpace.cxx:220
TGenPhaseSpace.cxx:221
TGenPhaseSpace.cxx:222
TGenPhaseSpace.cxx:223
TGenPhaseSpace.cxx:224
TGenPhaseSpace.cxx:225
TGenPhaseSpace.cxx:226
TGenPhaseSpace.cxx:227
TGenPhaseSpace.cxx:228
TGenPhaseSpace.cxx:229
TGenPhaseSpace.cxx:230
TGenPhaseSpace.cxx:231
TGenPhaseSpace.cxx:232
TGenPhaseSpace.cxx:233
TGenPhaseSpace.cxx:234
TGenPhaseSpace.cxx:235
TGenPhaseSpace.cxx:236
TGenPhaseSpace.cxx:237
TGenPhaseSpace.cxx:238
TGenPhaseSpace.cxx:239
TGenPhaseSpace.cxx:240
TGenPhaseSpace.cxx:241
TGenPhaseSpace.cxx:242