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
30
31////////////////////////////////////////////////////////////////////////////////
32/// The PDK function.
33
35{
36 Double_t x = (a-b-c)*(a+b+c)*(a-b+c)*(a+b-c);
37 x = TMath::Sqrt(x)/(2*a);
38 return x;
39}
40
41////////////////////////////////////////////////////////////////////////////////
42/// Special max function
43
44Int_t DoubleMax(const void *a, const void *b)
45{
46 Double_t aa = * ((Double_t *) a);
47 Double_t bb = * ((Double_t *) b);
48 if (aa > bb) return 1;
49 if (aa < bb) return -1;
50 return 0;
51
52}
53
54////////////////////////////////////////////////////////////////////////////////
55/// Copy constructor
56
58{
59 fNt = gen.fNt;
60 fWtMax = gen.fWtMax;
61 fTeCmTm = gen.fTeCmTm;
62 fBeta[0] = gen.fBeta[0];
63 fBeta[1] = gen.fBeta[1];
64 fBeta[2] = gen.fBeta[2];
65 for (Int_t i=0;i<fNt;i++) {
66 fMass[i] = gen.fMass[i];
67 fDecPro[i] = gen.fDecPro[i];
68 }
69}
70
71
72////////////////////////////////////////////////////////////////////////////////
73/// Assignment operator
74
76{
78 fNt = gen.fNt;
79 fWtMax = gen.fWtMax;
80 fTeCmTm = gen.fTeCmTm;
81 fBeta[0] = gen.fBeta[0];
82 fBeta[1] = gen.fBeta[1];
83 fBeta[2] = gen.fBeta[2];
84 for (Int_t i=0;i<fNt;i++) {
85 fMass[i] = gen.fMass[i];
86 fDecPro[i] = gen.fDecPro[i];
87 }
88 return *this;
89}
90
91////////////////////////////////////////////////////////////////////////////////
92/// Generate a random final state.
93/// The function returns the un-normalized weight of the current event.
94/// The TLorentzVector of each decay product can be obtained using GetDecay(n).
95///
96/// Note that Momentum, Energy units are Gev/C, GeV
97
99{
101 rno[0] = 0;
102 Int_t n;
103 if (fNt>2) {
104 for (n=1; n<fNt-1; n++) rno[n]=gRandom->Rndm(); // fNt-2 random numbers
105 qsort(rno+1 ,fNt-2 ,sizeof(Double_t) ,DoubleMax); // sort them
106 }
107 rno[fNt-1] = 1;
108
110 for (n=0; n<fNt; n++) {
111 sum += fMass[n];
112 invMas[n] = rno[n]*fTeCmTm + sum;
113 }
114
115 //
116 //-----> compute the un-normalized weight of the current event
117 //
120 for (n=0; n<fNt-1; n++) {
121 pd[n] = PDK(invMas[n+1],invMas[n],fMass[n+1]);
122 wt *= pd[n];
123 }
124
125 //
126 //-----> complete specification of event (Raubold-Lynch method)
127 //
128 fDecPro[0].SetPxPyPzE(0, pd[0], 0 , TMath::Sqrt(pd[0]*pd[0]+fMass[0]*fMass[0]) );
129
130 Int_t i=1;
131 Int_t j;
132 while (true) {
133 fDecPro[i].SetPxPyPzE(0, -pd[i-1], 0 , TMath::Sqrt(pd[i-1]*pd[i-1]+fMass[i]*fMass[i]) );
134
135 Double_t cZ = 2*gRandom->Rndm() - 1;
140 for (j=0; j<=i; j++) {
142 Double_t x = v->Px();
143 Double_t y = v->Py();
144 v->SetPx( cZ*x - sZ*y );
145 v->SetPy( sZ*x + cZ*y ); // rotation around Z
146 x = v->Px();
147 Double_t z = v->Pz();
148 v->SetPx( cY*x - sY*z );
149 v->SetPz( sY*x + cY*z ); // rotation around Y
150 }
151
152 if (i == (fNt-1)) break;
153
154 Double_t beta = pd[i] / sqrt(pd[i]*pd[i] + invMas[i]*invMas[i]);
155 for (j=0; j<=i; j++) fDecPro[j].Boost(0,beta,0);
156 i++;
157 }
158
159 //
160 //---> final boost of all particles
161 //
162 for (n=0;n<fNt;n++) fDecPro[n].Boost(fBeta[0],fBeta[1],fBeta[2]);
163
164 //
165 //---> return the un-normalized weight of event
166 //
167 return wt;
168}
169
170////////////////////////////////////////////////////////////////////////////////
171/// Return Lorentz vector corresponding to decay n
172
174{
175 if (n>fNt) return nullptr;
176 return fDecPro+n;
177}
178
179
180////////////////////////////////////////////////////////////////////////////////
181/// Input:
182/// - TLorentzVector &P: decay particle (Momentum, Energy units are Gev/C, GeV)
183/// - Int_t nt: number of decay products
184/// - Double_t *mass: array of decay product masses
185/// - Option_t *opt: default -> constant cross section
186/// "Fermi" -> Fermi energy dependence
187/// Return value:
188/// - kTRUE: the decay is permitted by kinematics
189/// - kFALSE: the decay is forbidden by kinematics
190///
191
193 const Double_t *mass, Option_t *opt)
194{
195 Int_t n;
196 fNt = nt;
197 if (fNt<2 || fNt>18) return kFALSE; // no more then 18 particle
198
199 //
200 //
201 //
202 fTeCmTm = P.Mag(); // total energy in C.M. minus the sum of the masses
203 for (n=0;n<fNt;n++) {
204 fMass[n] = mass[n];
205 fTeCmTm -= mass[n];
206 }
207
208 if (fTeCmTm<=0) return kFALSE; // not enough energy for this decay
209
210 //
211 //------> the max weight depends on opt:
212 // opt == "Fermi" --> fermi energy dependence for cross section
213 // else --> constant cross section as function of TECM (default)
214 //
215 if (strcasecmp(opt,"fermi")==0) {
216 // ffq[] = pi * (2*pi)**(FNt-2) / (FNt-2)!
217 Double_t ffq[] = {0
218 ,3.141592, 19.73921, 62.01255, 129.8788, 204.0131
219 ,256.3704, 268.4705, 240.9780, 189.2637
220 ,132.1308, 83.0202, 47.4210, 24.8295
221 ,12.0006, 5.3858, 2.2560, 0.8859 };
222 fWtMax = TMath::Power(fTeCmTm,fNt-2) * ffq[fNt-1] / P.Mag();
223
224 } else {
226 Double_t emmin = 0;
227 Double_t wtmax = 1;
228 for (n=1; n<fNt; n++) {
229 emmin += fMass[n-1];
230 emmax += fMass[n];
231 wtmax *= PDK(emmax, emmin, fMass[n]);
232 }
233 fWtMax = 1/wtmax;
234 }
235
236 //
237 //----> save the betas of the decaying particle
238 //
239 if (P.Beta()) {
240 Double_t w = P.Beta()/P.Rho();
241 fBeta[0] = P(0)*w;
242 fBeta[1] = P(1)*w;
243 fBeta[2] = P(2)*w;
244 }
245 else fBeta[0]=fBeta[1]=fBeta[2]=0;
246
247 return kTRUE;
248}
#define b(i)
Definition RSha256.hxx:100
#define c(i)
Definition RSha256.hxx:101
#define a(i)
Definition RSha256.hxx:99
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:59
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
const char Option_t
Option string (const char)
Definition RtypesCore.h:80
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
const Int_t kMAXP
Int_t DoubleMax(const void *a, const void *b)
Special max function.
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
Phase Space Generator, based on the GENBOD routine of CERNLIB.
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) noexcept
TObject assignment operator.
Definition TObject.h:299
Double_t Rndm() override
Machine independent random number generator.
Definition TRandom.cxx:558
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:673
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:732
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:605
constexpr Double_t Pi()
Definition TMath.h:40
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:599
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2339