Logo ROOT  
Reference Guide
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 \ingroup Physics
6
7 Utility class to generate n-body event,
8 with constant cross-section (default)
9 or with Fermi energy dependence (opt="Fermi").
10 The event is generated in the center-of-mass frame,
11 but the decay products are finally boosted
12 using the betas of the original particle.
13
14 The code is based on the GENBOD function (W515 from CERNLIB)
15 using the Raubold and Lynch method
16 F. James, Monte Carlo Phase Space, CERN 68-15 (1968)
17
18see example of use in PhaseSpace.C
19
20Note that Momentum, Energy units are Gev/C, GeV
21*/
22
23#include "TGenPhaseSpace.h"
24#include "TRandom.h"
25#include "TMath.h"
26
27const Int_t kMAXP = 18;
28
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 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{
100 Double_t rno[kMAXP];
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
109 Double_t invMas[kMAXP], sum=0;
110 for (n=0; n<fNt; n++) {
111 sum += fMass[n];
112 invMas[n] = rno[n]*fTeCmTm + sum;
113 }
114
115 //
116 //-----> compute the weight of the current event
117 //
118 Double_t wt=fWtMax;
119 Double_t pd[kMAXP];
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 (1) {
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;
136 Double_t sZ = TMath::Sqrt(1-cZ*cZ);
137 Double_t angY = 2*TMath::Pi() * gRandom->Rndm();
138 Double_t cY = TMath::Cos(angY);
139 Double_t sY = TMath::Sin(angY);
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 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 0;
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 {
225 Double_t emmax = fTeCmTm + fMass[0];
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 c(i)
Definition: RSha256.hxx:101
int Int_t
Definition: RtypesCore.h:45
const Bool_t kFALSE
Definition: RtypesCore.h:101
double Double_t
Definition: RtypesCore.h:59
const Bool_t kTRUE
Definition: RtypesCore.h:100
const char Option_t
Definition: RtypesCore.h:66
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
const Int_t kMAXP
ClassImp(TGenPhaseSpace)
Int_t DoubleMax(const void *a, const void *b)
Special max function.
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
Utility class to generate n-body event, with constant cross-section (default) or with Fermi energy de...
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 fTeCmTm
Double_t fMass[18]
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:552
double beta(double x, double y)
Calculates the beta function.
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
static double P[]
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition: TMath.h:660
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition: TMath.h:719
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition: TMath.h:592
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:586
TArc a
Definition: textangle.C:12
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345