Line data Source code
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>
30 : using std::endl;
31 :
32 :
33 : double EvtPawt(double a,double b,double c)
34 : {
35 0 : double temp=(a*a-(b+c)*(b+c))*(a*a-(b-c)*(b-c));
36 :
37 0 : if (temp<=0) {
38 0 : return 0.0;
39 : }
40 :
41 0 : return sqrt(temp)/(2.0*a);
42 0 : }
43 :
44 :
45 : double EvtGenKine::PhaseSpace( int ndaug, double mass[30], EvtVector4R p4[30],
46 : double mp )
47 :
48 : // N body phase space routine. Send parent with
49 : // daughters already defined ( Number and masses )
50 : // Returns four vectors in parent frame.
51 :
52 : {
53 :
54 : double energy, p3, alpha, beta;
55 :
56 0 : if ( ndaug == 1 ) {
57 0 : p4[0].set(mass[0],0.0,0.0,0.0);
58 0 : return 1.0;
59 : }
60 :
61 :
62 0 : if ( ndaug == 2 ) {
63 :
64 : //Two body phase space
65 :
66 0 : energy = ( mp*mp + mass[0]*mass[0] -
67 0 : mass[1]*mass[1] ) / ( 2.0 * mp );
68 :
69 : p3 = 0.0;
70 0 : if (energy > mass[0]) {
71 0 : p3 = sqrt( energy*energy - mass[0]*mass[0] );
72 0 : }
73 :
74 0 : p4[0].set( energy, 0.0, 0.0, p3 );
75 :
76 0 : energy = mp - energy;
77 0 : p3 = -1.0*p3;
78 0 : p4[1].set( energy, 0.0, 0.0, p3 );
79 :
80 : //Now rotate four vectors.
81 :
82 0 : alpha = EvtRandom::Flat( EvtConst::twoPi );
83 0 : beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
84 :
85 0 : p4[0].applyRotateEuler( alpha, beta, -alpha );
86 0 : p4[1].applyRotateEuler( alpha, beta, -alpha );
87 :
88 0 : return 1.0;
89 : }
90 :
91 0 : if ( ndaug != 2 ) {
92 :
93 : double wtmax=0.0;
94 0 : double pm[5][30],pmin,pmax,psum,rnd[30];
95 0 : double ran,wt,pa,costh,sinth,phi,p[4][30],be[4],bep,temp;
96 : int i,il,ilr,i1,il1u,il1,il2r,ilu;
97 : int il2=0;
98 :
99 0 : for(i=0;i<ndaug;i++){
100 0 : pm[4][i]=0.0;
101 0 : rnd[i]=0.0;
102 : }
103 :
104 0 : pm[0][0]=mp;
105 0 : pm[1][0]=0.0;
106 0 : pm[2][0]=0.0;
107 0 : pm[3][0]=0.0;
108 0 : pm[4][0]=mp;
109 :
110 : psum=0.0;
111 0 : for(i=1;i<ndaug+1;i++){
112 0 : psum=psum+mass[i-1];
113 : }
114 :
115 0 : pm[4][ndaug-1]=mass[ndaug-1];
116 :
117 0 : switch (ndaug) {
118 :
119 : case 1:
120 : wtmax=1.0/16.0;
121 0 : break;
122 : case 2:
123 : wtmax=1.0/150.0;
124 0 : break;
125 : case 3:
126 : wtmax=1.0/2.0;
127 0 : break;
128 : case 4:
129 : wtmax=1.0/5.0;
130 0 : break;
131 : case 5:
132 : wtmax=1.0/15.0;
133 0 : break;
134 : case 6:
135 : wtmax=1.0/15.0;
136 0 : break;
137 : case 7:
138 : wtmax=1.0/15.0;
139 0 : break;
140 : case 8:
141 : wtmax=1.0/15.0;
142 0 : break;
143 : case 9:
144 : wtmax=1.0/15.0;
145 0 : break;
146 : case 10:
147 : wtmax=1.0/15.0;
148 0 : break;
149 : case 11:
150 : wtmax=1.0/15.0;
151 0 : break;
152 : case 12:
153 : wtmax=1.0/15.0;
154 0 : break;
155 : case 13:
156 : wtmax=1.0/15.0;
157 0 : break;
158 : case 14:
159 : wtmax=1.0/15.0;
160 0 : break;
161 : case 15:
162 : wtmax=1.0/15.0;
163 0 : break;
164 : default:
165 0 : report(Severity::Error,"EvtGen") << "too many daughters for phase space..." << ndaug << " "<< mp <<endl;;
166 0 : break;
167 : }
168 :
169 0 : pmax=mp-psum+mass[ndaug-1];
170 :
171 : pmin=0.0;
172 :
173 0 : for(ilr=2;ilr<ndaug+1;ilr++){
174 0 : il=ndaug+1-ilr;
175 0 : pmax=pmax+mass[il-1];
176 0 : pmin=pmin+mass[il+1-1];
177 0 : wtmax=wtmax*EvtPawt(pmax,pmin,mass[il-1]);
178 : }
179 :
180 : do{
181 :
182 0 : rnd[0]=1.0;
183 : il1u=ndaug-1;
184 :
185 0 : for (il1=2;il1<il1u+1;il1++){
186 0 : ran=EvtRandom::Flat();
187 0 : for (il2r=2;il2r<il1+1;il2r++){
188 0 : il2=il1+1-il2r;
189 0 : if (ran<=rnd[il2-1]) goto two39;
190 0 : rnd[il2+1-1]=rnd[il2-1];
191 : }
192 : two39:
193 0 : rnd[il2+1-1]=ran;
194 : }
195 :
196 0 : rnd[ndaug-1]=0.0;
197 : wt=1.0;
198 0 : for(ilr=2;ilr<ndaug+1;ilr++){
199 0 : il=ndaug+1-ilr;
200 0 : pm[4][il-1]=pm[4][il+1-1]+mass[il-1]+
201 0 : (rnd[il-1]-rnd[il+1-1])*(mp-psum);
202 0 : wt=wt*EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
203 : }
204 0 : if (wt>wtmax) {
205 0 : report(Severity::Error,"EvtGen") << "wtmax to small in EvtPhaseSpace with "
206 0 : << ndaug <<" daughters"<<endl;;
207 0 : }
208 0 : } while (wt<EvtRandom::Flat(wtmax));
209 :
210 : ilu=ndaug-1;
211 :
212 0 : for (il=1;il<ilu+1;il++){
213 0 : pa=EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
214 0 : costh=EvtRandom::Flat(-1.0,1.0);
215 0 : sinth=sqrt(1.0-costh*costh);
216 0 : phi=EvtRandom::Flat(EvtConst::twoPi);
217 0 : p[1][il-1]=pa*sinth*cos(phi);
218 0 : p[2][il-1]=pa*sinth*sin(phi);
219 0 : p[3][il-1]=pa*costh;
220 0 : pm[1][il+1-1]=-p[1][il-1];
221 0 : pm[2][il+1-1]=-p[2][il-1];
222 0 : pm[3][il+1-1]=-p[3][il-1];
223 0 : p[0][il-1]=sqrt(pa*pa+mass[il-1]*mass[il-1]);
224 0 : pm[0][il+1-1]=sqrt(pa*pa+pm[4][il+1-1]*pm[4][il+1-1]);
225 : }
226 :
227 0 : p[0][ndaug-1]=pm[0][ndaug-1];
228 0 : p[1][ndaug-1]=pm[1][ndaug-1];
229 0 : p[2][ndaug-1]=pm[2][ndaug-1];
230 0 : p[3][ndaug-1]=pm[3][ndaug-1];
231 :
232 0 : for (ilr=2;ilr<ndaug+1;ilr++){
233 0 : il=ndaug+1-ilr;
234 0 : be[0]=pm[0][il-1]/pm[4][il-1];
235 0 : be[1]=pm[1][il-1]/pm[4][il-1];
236 0 : be[2]=pm[2][il-1]/pm[4][il-1];
237 0 : be[3]=pm[3][il-1]/pm[4][il-1];
238 :
239 0 : for (i1=il;i1<ndaug+1;i1++){
240 0 : bep=be[1]*p[1][i1-1]+be[2]*p[2][i1-1]+
241 0 : be[3]*p[3][i1-1]+be[0]*p[0][i1-1];
242 0 : temp=(p[0][i1-1]+bep)/(be[0]+1.0);
243 0 : p[1][i1-1]=p[1][i1-1]+temp*be[1];
244 0 : p[2][i1-1]=p[2][i1-1]+temp*be[2];
245 0 : p[3][i1-1]=p[3][i1-1]+temp*be[3];
246 0 : p[0][i1-1]=bep;
247 :
248 : }
249 : }
250 :
251 0 : for (ilr=0;ilr<ndaug;ilr++){
252 0 : p4[ilr].set(p[0][ilr],p[1][ilr],p[2][ilr],p[3][ilr]);
253 : }
254 :
255 : return 1.0;
256 0 : }
257 :
258 0 : return 1.0;
259 :
260 0 : }
261 :
262 :
263 : double EvtGenKine::PhaseSpacePole(double M, double m1, double m2, double m3,
264 : double a,EvtVector4R p4[10])
265 :
266 : // generate kinematics for 3 body decays, pole for the m1,m2 mass.
267 :
268 : {
269 :
270 : //f1 = 1 (phasespace)
271 : //f2 = a*(1/m12sq)^2
272 :
273 0 : double m12sqmax=(M-m3)*(M-m3);
274 0 : double m12sqmin=(m1+m2)*(m1+m2);
275 :
276 0 : double m13sqmax=(M-m2)*(M-m2);
277 0 : double m13sqmin=(m1+m3)*(m1+m3);
278 :
279 0 : double v1=(m12sqmax-m12sqmin)*(m13sqmax-m13sqmin);
280 0 : double v2= a*(1.0/m12sqmin-1.0/m12sqmax)*(m13sqmax-m13sqmin);
281 :
282 : double m12sq,m13sq;
283 :
284 0 : double r=v1/(v1+v2);
285 :
286 : double m13min,m13max;
287 :
288 0 : do{
289 :
290 0 : m13sq=EvtRandom::Flat(m13sqmin,m13sqmax);
291 :
292 0 : if (r>EvtRandom::Flat()){
293 0 : m12sq=EvtRandom::Flat(m12sqmin,m12sqmax);
294 0 : }
295 : else{
296 0 : m12sq=1.0/(1.0/m12sqmin-EvtRandom::Flat()*(1.0/m12sqmin-1.0/m12sqmax));
297 : }
298 :
299 : //kinematically allowed?
300 0 : double E3star=(M*M-m12sq-m3*m3)/sqrt(4*m12sq);
301 0 : double E1star=(m12sq+m1*m1-m2*m2)/sqrt(4*m12sq);
302 0 : double p3star=sqrt(E3star*E3star-m3*m3);
303 0 : double p1star=sqrt(E1star*E1star-m1*m1);
304 0 : m13max=(E3star+E1star)*(E3star+E1star)-
305 0 : (p3star-p1star)*(p3star-p1star);
306 0 : m13min=(E3star+E1star)*(E3star+E1star)-
307 0 : (p3star+p1star)*(p3star+p1star);
308 :
309 0 : }while(m13sq<m13min||m13sq>m13max);
310 :
311 :
312 0 : double E2=(M*M+m2*m2-m13sq)/(2.0*M);
313 0 : double E3=(M*M+m3*m3-m12sq)/(2.0*M);
314 0 : double E1=M-E2-E3;
315 0 : double p1mom=sqrt(E1*E1-m1*m1);
316 0 : double p3mom=sqrt(E3*E3-m3*m3);
317 0 : double cost13=(2.0*E1*E3+m1*m1+m3*m3-m13sq)/(2.0*p1mom*p3mom);
318 :
319 : //report(Severity::Info,"EvtGen") << m13sq << endl;
320 : //report(Severity::Info,"EvtGen") << m12sq << endl;
321 : //report(Severity::Info,"EvtGen") << E1 << endl;
322 : //report(Severity::Info,"EvtGen") << E2 << endl;
323 : //report(Severity::Info,"EvtGen") << E3 << endl;
324 : //report(Severity::Info,"EvtGen") << p1mom << endl;
325 : //report(Severity::Info,"EvtGen") << p3mom << endl;
326 : //report(Severity::Info,"EvtGen") << cost13 << endl;
327 :
328 :
329 0 : p4[2].set(E3,0.0,0.0,p3mom);
330 0 : p4[0].set(E1,p1mom*sqrt(1.0-cost13*cost13),0.0,p1mom*cost13);
331 0 : p4[1].set(E2,-p1mom*sqrt(1.0-cost13*cost13),0.0,-p1mom*cost13-p3mom);
332 :
333 : //report(Severity::Info,"EvtGen") << "p4:"<<p4[0]<<p4[1]<<p4[2]<<endl;
334 :
335 0 : double alpha = EvtRandom::Flat( EvtConst::twoPi );
336 0 : double beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
337 0 : double gamma = EvtRandom::Flat( EvtConst::twoPi );
338 :
339 0 : p4[0].applyRotateEuler( alpha, beta, gamma );
340 0 : p4[1].applyRotateEuler( alpha, beta, gamma );
341 0 : p4[2].applyRotateEuler( alpha, beta, gamma );
342 :
343 0 : return 1.0+a/(m12sq*m12sq);
344 :
345 : }
346 :
|
