CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtGenKine.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of the EvtGen package developed jointly
5// for the BaBar and CLEO collaborations. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/COPYRIGHT
9// Copyright (C) 1998 Caltech, UCSB
10//
11// Module: EvtGenKine.cc
12//
13// Description: Tools for generating distributions of four vectors in
14// phasespace
15//
16// Modification history:
17//
18// DJL/RYD September 25, 1996 Module created
19//
20//------------------------------------------------------------------------
21//
22#include "EvtGenBase/EvtPatches.hh"
23#include <iostream>
24#include "EvtGenBase/EvtGenKine.hh"
25#include "EvtGenBase/EvtRandom.hh"
26#include "EvtGenBase/EvtVector4R.hh"
27#include "EvtGenBase/EvtReport.hh"
28#include "EvtGenBase/EvtConst.hh"
29#include <math.h>
30using std::endl;
31
32
33double EvtPawt(double a,double b,double c)
34{
35 double temp=(a*a-(b+c)*(b+c))*(a*a-(b-c)*(b-c));
36
37 if (temp<=0) {
38
39 //report(ERROR,"EvtGen")<<"Sqrt of negative number in EvtPhaseSpace\n"<<
40 // "This seems to happen on AIX but I do not know why yet!"<<endl;
41
42 return 0.0;
43
44 }
45
46 return sqrt(temp)/(2.0*a);
47}
48
49
50double EvtGenKine::PhaseSpace( int ndaug, double mass[30], EvtVector4R p4[30],
51 double mp )
52
53// N body phase space routine. Send parent with
54// daughters already defined ( Number and masses )
55// Returns four vectors in parent frame.
56
57{
58
59 double energy, p3, alpha, beta;
60
61 if ( ndaug == 1 ) {
62 p4[0].set(mass[0],0.0,0.0,0.0);
63 return 1.0;
64 }
65
66
67 if ( ndaug == 2 ) {
68
69// Two body phase space
70
71 energy = ( mp*mp + mass[0]*mass[0] -
72 mass[1]*mass[1] ) / ( 2.0 * mp );
73
74 p3 = sqrt( energy*energy - mass[0]*mass[0] );
75
76 p4[0].set( energy, 0.0, 0.0, p3 );
77
78 energy = mp - energy;
79 p3 = -1.0*p3;
80 p4[1].set( energy, 0.0, 0.0, p3 );
81
82// Now rotate four vectors.
83
84 alpha = EvtRandom::Flat( EvtConst::twoPi );
85 beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
86
87 p4[0].applyRotateEuler( alpha, beta, -alpha );
88 p4[1].applyRotateEuler( alpha, beta, -alpha );
89
90 return 1.0;
91 }
92
93 if ( ndaug != 2 ) {
94
95 double wtmax=0.0;
96 double pm[5][30],to[4],pmin,pmax,psum,rnd[30];
97 double ran,wt,pa,costh,sinth,phi,p[4][30],be[4],bep,temp;
98 int i,il,ilr,i1,il1u,il1,il2r,ilu;
99 int il2=0;
100
101 for(i=0;i<ndaug;i++){
102 pm[4][i]=0.0;
103 rnd[i]=0.0;
104 }
105
106 pm[0][0]=mp;
107 pm[1][0]=0.0;
108 pm[2][0]=0.0;
109 pm[3][0]=0.0;
110 pm[4][0]=mp;
111
112 to[0]=mp;
113 to[1]=0.0;
114 to[2]=0.0;
115 to[3]=0.0;
116
117 psum=0.0;
118 for(i=1;i<ndaug+1;i++){
119 psum=psum+mass[i-1];
120 }
121
122 pm[4][ndaug-1]=mass[ndaug-1];
123
124 switch (ndaug) {
125
126 case 1:
127 wtmax=1.0/16.0;
128 break;
129 case 2:
130 wtmax=1.0/150.0;
131 break;
132 case 3:
133 wtmax=1.0/2.0;
134 break;
135 case 4:
136 wtmax=1.0/5.0;
137 break;
138 case 5:
139 wtmax=1.0/15.0;
140 break;
141 case 6:
142 wtmax=1.0/15.0;
143 break;
144 case 7:
145 wtmax=1.0/15.0;
146 break;
147 case 8:
148 wtmax=1.0/15.0;
149 break;
150 case 9:
151 wtmax=1.0/15.0;
152 break;
153 case 10:
154 wtmax=1.0/15.0;
155 break;
156 case 11:
157 wtmax=1.0/15.0;
158 break;
159 case 12:
160 wtmax=1.0/15.0;
161 break;
162 case 13:
163 wtmax=1.0/15.0;
164 break;
165 case 14:
166 wtmax=1.0/15.0;
167 break;
168 case 15:
169 wtmax=1.0/15.0;
170 break;
171 default:
172 report(ERROR,"EvtGen") << "too many daughters for phase space..." << ndaug << " "<< mp <<endl;;
173 break;
174 }
175
176 pmax=mp-psum+mass[ndaug-1];
177
178 pmin=0.0;
179
180 for(ilr=2;ilr<ndaug+1;ilr++){
181 il=ndaug+1-ilr;
182 pmax=pmax+mass[il-1];
183 pmin=pmin+mass[il+1-1];
184 wtmax=wtmax*EvtPawt(pmax,pmin,mass[il-1]);
185 }
186
187 //report(INFO,"EvtGen") << "wtmax:"<<wtmax<<endl;
188
189 do{
190
191 rnd[0]=1.0;
192 il1u=ndaug-1;
193
194 for (il1=2;il1<il1u+1;il1++){
195 ran=EvtRandom::Flat();
196 for (il2r=2;il2r<il1+1;il2r++){
197 il2=il1+1-il2r;
198 if (ran<=rnd[il2-1]) goto two39;
199 rnd[il2+1-1]=rnd[il2-1];
200 }
201 two39:
202 rnd[il2+1-1]=ran;
203 }
204
205 rnd[ndaug-1]=0.0;
206 wt=1.0;
207 for(ilr=2;ilr<ndaug+1;ilr++){
208 il=ndaug+1-ilr;
209 pm[4][il-1]=pm[4][il+1-1]+mass[il-1]+
210 (rnd[il-1]-rnd[il+1-1])*(mp-psum);
211 wt=wt*EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
212 }
213 if (wt>wtmax) {
214 report(ERROR,"EvtGen") << "wtmax to small in EvtPhaseSpace with "
215 << ndaug <<" daughters"<<endl;;
216 }
217 } while (wt<EvtRandom::Flat(wtmax));
218
219 ilu=ndaug-1;
220
221 for (il=1;il<ilu+1;il++){
222 pa=EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
223 costh=EvtRandom::Flat(-1.0,1.0);
224 sinth=sqrt(1.0-costh*costh);
225 phi=EvtRandom::Flat(EvtConst::twoPi);
226 p[1][il-1]=pa*sinth*cos(phi);
227 p[2][il-1]=pa*sinth*sin(phi);
228 p[3][il-1]=pa*costh;
229 pm[1][il+1-1]=-p[1][il-1];
230 pm[2][il+1-1]=-p[2][il-1];
231 pm[3][il+1-1]=-p[3][il-1];
232 p[0][il-1]=sqrt(pa*pa+mass[il-1]*mass[il-1]);
233 pm[0][il+1-1]=sqrt(pa*pa+pm[4][il+1-1]*pm[4][il+1-1]);
234 }
235
236 p[0][ndaug-1]=pm[0][ndaug-1];
237 p[1][ndaug-1]=pm[1][ndaug-1];
238 p[2][ndaug-1]=pm[2][ndaug-1];
239 p[3][ndaug-1]=pm[3][ndaug-1];
240
241 for (ilr=2;ilr<ndaug+1;ilr++){
242 il=ndaug+1-ilr;
243 be[0]=pm[0][il-1]/pm[4][il-1];
244 be[1]=pm[1][il-1]/pm[4][il-1];
245 be[2]=pm[2][il-1]/pm[4][il-1];
246 be[3]=pm[3][il-1]/pm[4][il-1];
247
248 for (i1=il;i1<ndaug+1;i1++){
249 bep=be[1]*p[1][i1-1]+be[2]*p[2][i1-1]+
250 be[3]*p[3][i1-1]+be[0]*p[0][i1-1];
251 temp=(p[0][i1-1]+bep)/(be[0]+1.0);
252 p[1][i1-1]=p[1][i1-1]+temp*be[1];
253 p[2][i1-1]=p[2][i1-1]+temp*be[2];
254 p[3][i1-1]=p[3][i1-1]+temp*be[3];
255 p[0][i1-1]=bep;
256
257 }
258 }
259
260 for (ilr=0;ilr<ndaug;ilr++){
261 p4[ilr].set(p[0][ilr],p[1][ilr],p[2][ilr],p[3][ilr]);
262 }
263
264 return 1.0;
265 }
266
267 return 1.0;
268
269}
270
271
272double EvtGenKine::PhaseSpacePole(double M, double m1, double m2, double m3,
273 double a,EvtVector4R p4[10])
274
275// generate kinematics for 3 body decays, pole for the m1,m2 mass.
276
277{
278
279 //f1 = 1 (phasespace)
280 //f2 = a*(1/m12sq)^2
281
282 double m12sqmax=(M-m3)*(M-m3);
283 double m12sqmin=(m1+m2)*(m1+m2);
284
285 double m13sqmax=(M-m2)*(M-m2);
286 double m13sqmin=(m1+m3)*(m1+m3);
287
288 double v1=(m12sqmax-m12sqmin)*(m13sqmax-m13sqmin);
289 double v2= a*(1.0/m12sqmin-1.0/m12sqmax)*(m13sqmax-m13sqmin);
290
291 double m12sq,m13sq;
292
293 double r=v1/(v1+v2);
294
295 //report(INFO,"EvtGen") << "v1,v2:"<<v1<<","<<v2<<endl;
296
297 double m13min,m13max;
298
299 do{
300
301 m13sq=EvtRandom::Flat(m13sqmin,m13sqmax);
302
303 if (r>EvtRandom::Flat()){
304 m12sq=EvtRandom::Flat(m12sqmin,m12sqmax);
305 }
306 else{
307 m12sq=1.0/(1.0/m12sqmin-EvtRandom::Flat()*(1.0/m12sqmin-1.0/m12sqmax));
308 }
309
310 //kinematically allowed?
311 double E3star=(M*M-m12sq-m3*m3)/sqrt(4*m12sq);
312 double E1star=(m12sq+m1*m1-m2*m2)/sqrt(4*m12sq);
313 double p3star=sqrt(E3star*E3star-m3*m3);
314 double p1star=sqrt(E1star*E1star-m1*m1);
315 m13max=(E3star+E1star)*(E3star+E1star)-
316 (p3star-p1star)*(p3star-p1star);
317 m13min=(E3star+E1star)*(E3star+E1star)-
318 (p3star+p1star)*(p3star+p1star);
319
320 }while(m13sq<m13min||m13sq>m13max);
321
322
323 double E2=(M*M+m2*m2-m13sq)/(2.0*M);
324 double E3=(M*M+m3*m3-m12sq)/(2.0*M);
325 double E1=M-E2-E3;
326 double p1mom=sqrt(E1*E1-m1*m1);
327 double p3mom=sqrt(E3*E3-m3*m3);
328 double cost13=(2.0*E1*E3+m1*m1+m3*m3-m13sq)/(2.0*p1mom*p3mom);
329
330 //report(INFO,"EvtGen") << m13sq << endl;
331 //report(INFO,"EvtGen") << m12sq << endl;
332 //report(INFO,"EvtGen") << E1 << endl;
333 //report(INFO,"EvtGen") << E2 << endl;
334 //report(INFO,"EvtGen") << E3 << endl;
335 //report(INFO,"EvtGen") << p1mom << endl;
336 //report(INFO,"EvtGen") << p3mom << endl;
337 //report(INFO,"EvtGen") << cost13 << endl;
338
339
340 p4[2].set(E3,0.0,0.0,p3mom);
341 p4[0].set(E1,p1mom*sqrt(1.0-cost13*cost13),0.0,p1mom*cost13);
342 p4[1].set(E2,-p1mom*sqrt(1.0-cost13*cost13),0.0,-p1mom*cost13-p3mom);
343
344 //report(INFO,"EvtGen") << "p4:"<<p4[0]<<p4[1]<<p4[2]<<endl;
345
346 double alpha = EvtRandom::Flat( EvtConst::twoPi );
347 double beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
348 double gamma = EvtRandom::Flat( EvtConst::twoPi );
349
350 p4[0].applyRotateEuler( alpha, beta, gamma );
351 p4[1].applyRotateEuler( alpha, beta, gamma );
352 p4[2].applyRotateEuler( alpha, beta, gamma );
353
354 return 1.0+a/(m12sq*m12sq);
355
356}
357
sin
double sin(const BesAngle a)
Definition BesAngle.h:210
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
mass
double mass
Definition CosmicGenerator.cxx:138
EvtPawt
double EvtPawt(double a, double b, double c)
Definition EvtGenKine.cc:33
report
ostream & report(Severity severity, const char *facility)
Definition EvtReport.cc:36
ERROR
@ ERROR
Definition EvtReport.hh:49
alpha
const double alpha
Definition FastVertexFit.cxx:4
energy
************Class m_ypar INTEGER m_KeyWgt INTEGER m_nphot INTEGER m_KeyGPS INTEGER m_IsBeamPolarized INTEGER m_EvtGenInterface DOUBLE PRECISION m_Emin DOUBLE PRECISION m_sphot DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_q2 DOUBLE PRECISION m_PolBeam2 DOUBLE PRECISION m_xErrPb *COMMON c_KK2f $ !CMS energy average $ !Spin Polarization vector first beam $ !Spin Polarization vector second beam $ !Beam energy spread[GeV] $ !minimum hadronization energy[GeV] $ !input READ never touch them !$ !debug facility $ !maximum weight $ !inverse alfaQED $ !minimum real photon energy
Definition KK2f.h:50
EvtConst::twoPi
static const double twoPi
Definition EvtConst.hh:29
EvtGenKine::PhaseSpacePole
static double PhaseSpacePole(double M, double m1, double m2, double m3, double a, EvtVector4R p4[10])
Definition EvtGenKine.cc:272
EvtGenKine::PhaseSpace
static double PhaseSpace(int ndaug, double mass[30], EvtVector4R p4[30], double mp)
Definition EvtGenKine.cc:50
EvtRandom::Flat
static double Flat()
Definition EvtRandom.cc:73
EvtVector4R::applyRotateEuler
void applyRotateEuler(double alpha, double beta, double gamma)
Definition EvtVector4R.cc:81
EvtVector4R::set
void set(int i, double d)
Definition EvtVector4R.hh:183
mp
const double mp
Definition incllambda.cxx:45