ROOT   Reference Guide
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 \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 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
const Bool_t kFALSE
Definition RtypesCore.h:92
double Double_t
Definition RtypesCore.h:59
const Bool_t kTRUE
Definition RtypesCore.h:91
const char Option_t
Definition RtypesCore.h:66
#define ClassImp(name)
Definition Rtypes.h:364
const Int_t kMAXP
Int_t DoubleMax(const void *a, const void *b)
Special max function.
double sqrt(double)
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 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:37
TObject & operator=(const TObject &rhs)
TObject assignment operator.
Definition TObject.h:283
virtual Double_t Rndm()
Machine independent random number generator.
Definition TRandom.cxx:552
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)
Definition TMath.h:691
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition TMath.h:735
Double_t Cos(Double_t)
Definition TMath.h:643
constexpr Double_t Pi()
Definition TMath.h:37
Double_t Sin(Double_t)
Definition TMath.h:639
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2345