Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
57 TGenPhaseSpace::TGenPhaseSpace(const TGenPhaseSpace &gen) : TObject(gen)
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 
75 TGenPhaseSpace& TGenPhaseSpace::operator=(const TGenPhaseSpace &gen)
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 
98 Double_t TGenPhaseSpace::Generate()
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 
173 TLorentzVector *TGenPhaseSpace::GetDecay(Int_t n)
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 
192 Bool_t TGenPhaseSpace::SetDecay(TLorentzVector &P, Int_t nt,
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 }