Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TGenPhaseSpace.cxx
Go to the documentation of this file.
1// @(#)root/physics:$Id$
2// Author: Rene Brun , Valerio Filippini 06/09/2000
3
4/** \class TGenPhaseSpace
5 \legacy{TGenPhaseSpace, No alternative classes are currently available.}
6 \ingroup Physics
7
8 Utility class to generate n-body event,
9 with constant cross-section (default)
10 or with Fermi energy dependence (opt="Fermi").
11 The event is generated in the center-of-mass frame,
12 but the decay products are finally boosted
13 using the betas of the original particle.
14
15 The code is based on the GENBOD function (W515 from CERNLIB)
16 using the Raubold and Lynch method
17 F. James, Monte Carlo Phase Space, CERN 68-15 (1968)
18
19see example of use in PhaseSpace.C
20
21Note that Momentum, Energy units are Gev/C, GeV
22*/
23
24#include "TGenPhaseSpace.h"
25#include "TRandom.h"
26#include "TMath.h"
27
28const Int_t kMAXP = 18;
29
31
32////////////////////////////////////////////////////////////////////////////////
33/// The PDK function.
34
36{
37 Double_t x = (a-b-c)*(a+b+c)*(a-b+c)*(a+b-c);
38 x = TMath::Sqrt(x)/(2*a);
39 return x;
40}
41
42////////////////////////////////////////////////////////////////////////////////
43/// Special max function
44
45Int_t DoubleMax(const void *a, const void *b)
46{
47 Double_t aa = * ((Double_t *) a);
48 Double_t bb = * ((Double_t *) b);
49 if (aa > bb) return 1;
50 if (aa < bb) return -1;
51 return 0;
52
53}
54
55////////////////////////////////////////////////////////////////////////////////
56/// Copy constructor
57
59{
60 fNt = gen.fNt;
61 fWtMax = gen.fWtMax;
62 fTeCmTm = gen.fTeCmTm;
63 fBeta[0] = gen.fBeta[0];
64 fBeta[1] = gen.fBeta[1];
65 fBeta[2] = gen.fBeta[2];
66 for (Int_t i=0;i<fNt;i++) {
67 fMass[i] = gen.fMass[i];
68 fDecPro[i] = gen.fDecPro[i];
69 }
70}
71
72
73////////////////////////////////////////////////////////////////////////////////
74/// Assignment operator
75
77{
79 fNt = gen.fNt;
80 fWtMax = gen.fWtMax;
81 fTeCmTm = gen.fTeCmTm;
82 fBeta[0] = gen.fBeta[0];
83 fBeta[1] = gen.fBeta[1];
84 fBeta[2] = gen.fBeta[2];
85 for (Int_t i=0;i<fNt;i++) {
86 fMass[i] = gen.fMass[i];
87 fDecPro[i] = gen.fDecPro[i];
88 }
89 return *this;
90}
91
92////////////////////////////////////////////////////////////////////////////////
93/// Generate a random final state.
94/// The function returns the weight of the current event.
95/// The TLorentzVector of each decay product can be obtained using GetDecay(n).
96///
97/// Note that Momentum, Energy units are Gev/C, GeV
98
100{
101 Double_t rno[kMAXP];
102 rno[0] = 0;
103 Int_t n;
104 if (fNt>2) {
105 for (n=1; n<fNt-1; n++) rno[n]=gRandom->Rndm(); // fNt-2 random numbers
106 qsort(rno+1 ,fNt-2 ,sizeof(Double_t) ,DoubleMax); // sort them
107 }
108 rno[fNt-1] = 1;
109
110 Double_t invMas[kMAXP], sum=0;
111 for (n=0; n<fNt; n++) {
112 sum += fMass[n];
113 invMas[n] = rno[n]*fTeCmTm + sum;
114 }
115
116 //
117 //-----> compute the weight of the current event
118 //
119 Double_t wt=fWtMax;
120 Double_t pd[kMAXP];
121 for (n=0; n<fNt-1; n++) {
122 pd[n] = PDK(invMas[n+1],invMas[n],fMass[n+1]);
123 wt *= pd[n];
124 }
125
126 //
127 //-----> complete specification of event (Raubold-Lynch method)
128 //
129 fDecPro[0].SetPxPyPzE(0, pd[0], 0 , TMath::Sqrt(pd[0]*pd[0]+fMass[0]*fMass[0]) );
130
131 Int_t i=1;
132 Int_t j;
133 while (true) {
134 fDecPro[i].SetPxPyPzE(0, -pd[i-1], 0 , TMath::Sqrt(pd[i-1]*pd[i-1]+fMass[i]*fMass[i]) );
135
136 Double_t cZ = 2*gRandom->Rndm() - 1;
137 Double_t sZ = TMath::Sqrt(1-cZ*cZ);
138 Double_t angY = 2*TMath::Pi() * gRandom->Rndm();
139 Double_t cY = TMath::Cos(angY);
140 Double_t sY = TMath::Sin(angY);
141 for (j=0; j<=i; j++) {
143 Double_t x = v->Px();
144 Double_t y = v->Py();
145 v->SetPx( cZ*x - sZ*y );
146 v->SetPy( sZ*x + cZ*y ); // rotation around Z
147 x = v->Px();
148 Double_t z = v->Pz();
149 v->SetPx( cY*x - sY*z );
150 v->SetPz( sY*x + cY*z ); // rotation around Y
151 }
152
153 if (i == (fNt-1)) break;
154
155 Double_t beta = pd[i] / sqrt(pd[i]*pd[i] + invMas[i]*invMas[i]);
156 for (j=0; j<=i; j++) fDecPro[j].Boost(0,beta,0);
157 i++;
158 }
159
160 //
161 //---> final boost of all particles
162 //
163 for (n=0;n<fNt;n++) fDecPro[n].Boost(fBeta[0],fBeta[1],fBeta[2]);
164
165 //
166 //---> return the weight of event
167 //
168 return wt;
169}
170
171////////////////////////////////////////////////////////////////////////////////
172/// Return Lorentz vector corresponding to decay n
173
175{
176 if (n>fNt) return nullptr;
177 return fDecPro+n;
178}
179
180
181////////////////////////////////////////////////////////////////////////////////
182/// Input:
183/// - TLorentzVector &P: decay particle (Momentum, Energy units are Gev/C, GeV)
184/// - Int_t nt: number of decay products
185/// - Double_t *mass: array of decay product masses
186/// - Option_t *opt: default -> constant cross section
187/// "Fermi" -> Fermi energy dependence
188/// Return value:
189/// - kTRUE: the decay is permitted by kinematics
190/// - kFALSE: the decay is forbidden by kinematics
191///
192
194 const Double_t *mass, Option_t *opt)
195{
196 Int_t n;
197 fNt = nt;
198 if (fNt<2 || fNt>18) return kFALSE; // no more then 18 particle
199
200 //
201 //
202 //
203 fTeCmTm = P.Mag(); // total energy in C.M. minus the sum of the masses
204 for (n=0;n<fNt;n++) {
205 fMass[n] = mass[n];
206 fTeCmTm -= mass[n];
207 }
208
209 if (fTeCmTm<=0) return kFALSE; // not enough energy for this decay
210
211 //
212 //------> the max weight depends on opt:
213 // opt == "Fermi" --> fermi energy dependence for cross section
214 // else --> constant cross section as function of TECM (default)
215 //
216 if (strcasecmp(opt,"fermi")==0) {
217 // ffq[] = pi * (2*pi)**(FNt-2) / (FNt-2)!
218 Double_t ffq[] = {0
219 ,3.141592, 19.73921, 62.01255, 129.8788, 204.0131
220 ,256.3704, 268.4705, 240.9780, 189.2637
221 ,132.1308, 83.0202, 47.4210, 24.8295
222 ,12.0006, 5.3858, 2.2560, 0.8859 };
223 fWtMax = TMath::Power(fTeCmTm,fNt-2) * ffq[fNt-1] / P.Mag();
224
225 } else {
226 Double_t emmax = fTeCmTm + fMass[0];
227 Double_t emmin = 0;
228 Double_t wtmax = 1;
229 for (n=1; n<fNt; n++) {
230 emmin += fMass[n-1];
231 emmax += fMass[n];
232 wtmax *= PDK(emmax, emmin, fMass[n]);
233 }
234 fWtMax = 1/wtmax;
235 }
236
237 //
238 //----> save the betas of the decaying particle
239 //
240 if (P.Beta()) {
241 Double_t w = P.Beta()/P.Rho();
242 fBeta[0] = P(0)*w;
243 fBeta[1] = P(1)*w;
244 fBeta[2] = P(2)*w;
245 }
246 else fBeta[0]=fBeta[1]=fBeta[2]=0;
247
248 return kTRUE;
249}
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
int Int_t
Definition RtypesCore.h:45
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
double Double_t
Definition RtypesCore.h:59
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:377
const Int_t kMAXP
Int_t DoubleMax(const void *a, const void *b)
Special max function.
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
<div class="legacybox"><h2>Legacy Code</h2> TGenPhaseSpace is a legacy interface: there will be no bu...
Double_t PDK(Double_t a, Double_t b, Double_t c)
The PDK function.
TLorentzVector fDecPro[18]
Double_t fBeta[3]
TGenPhaseSpace & operator=(const TGenPhaseSpace &gen)
Assignment operator.
Double_t Generate()
Generate a random final state.
TLorentzVector * GetDecay(Int_t n)
Return Lorentz vector corresponding to decay n.
Bool_t SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass, Option_t *opt="")
Input:
Double_t fMass[18]
Double_t Rho() const
Double_t Beta() const
Double_t Mag() const
void SetPxPyPzE(Double_t px, Double_t py, Double_t pz, Double_t e)
Mother of all ROOT objects.
Definition TObject.h:41
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition TObject.h:298
Double_t Rndm() override
Machine independent random number generator.
Definition TRandom.cxx:559
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:721
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:594
constexpr Double_t Pi()
Definition TMath.h:37
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:588
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345