Logo ROOT   6.13/01
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 
18 see example of use in PhaseSpace.C
19 
20 Note that Momentum, Energy units are Gev/C, GeV
21 */
22 
23 #include "TGenPhaseSpace.h"
24 #include "TRandom.h"
25 #include "TMath.h"
26 
27 const Int_t kMAXP = 18;
28 
29 ClassImp(TGenPhaseSpace);
30 
31 ////////////////////////////////////////////////////////////////////////////////
32 /// The PDK function.
33 
34 Double_t TGenPhaseSpace::PDK(Double_t a, Double_t b, Double_t c)
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 
44 Int_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 {
77  TObject::operator=(gen);
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++) {
141  TLorentzVector *v = fDecPro+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 }
void SetPxPyPzE(Double_t px, Double_t py, Double_t pz, Double_t e)
Double_t Rho() const
void SetPx(Double_t a)
Double_t fMass[18]
const Int_t kMAXP
Double_t Generate()
Generate a random final state.
you should not use this method at all Int_t y
Definition: TRolke.cxx:630
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:734
double beta(double x, double y)
Calculates the beta function.
double sqrt(double)
Double_t Beta() const
Double_t Pz() const
TGenPhaseSpace & operator=(const TGenPhaseSpace &gen)
Assignment operator.
TLorentzVector fDecPro[18]
void SetPy(Double_t a)
constexpr Double_t Pi()
Definition: TMath.h:38
Int_t DoubleMax(const void *a, const void *b)
Special max function.
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:533
TLorentzVector is a general four-vector class, which can be used either for the description of positi...
SVector< double, 2 > v
Definition: Dict.h:5
static double P[]
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
void SetPz(Double_t a)
Utility class to generate n-body event, with constant cross-section (default) or with Fermi energy de...
* x
Deprecated and error prone model selection interface.
Definition: TRolke.cxx:630
Double_t Cos(Double_t)
Definition: TMath.h:640
Double_t Py() const
Bool_t SetDecay(TLorentzVector &P, Int_t nt, const Double_t *mass, Option_t *opt="")
Input:
Double_t Px() const
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
Double_t PDK(Double_t a, Double_t b, Double_t c)
The PDK function.
TLorentzVector * GetDecay(Int_t n)
Return Lorentz vector corresponding to decay n.
Double_t Sin(Double_t)
Definition: TMath.h:636
Double_t fTeCmTm
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
Double_t Mag() const
Double_t Sqrt(Double_t x)
Definition: TMath.h:690
Double_t fBeta[3]