BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtGenKine Class Reference

#include <EvtGenKine.hh>

Static Public Member Functions

static double PhaseSpace (int ndaug, double mass[30], EvtVector4R p4[30], double mp)
static double PhaseSpacePole (double M, double m1, double m2, double m3, double a, EvtVector4R p4[10])

Detailed Description

Definition at line 25 of file EvtGenKine.hh.

Member Function Documentation

◆ PhaseSpace()

double EvtGenKine::PhaseSpace ( int ndaug,
double mass[30],
EvtVector4R p4[30],
double mp )
static

Definition at line 47 of file EvtGenKine.cc.

53{
54
55 double energy, p3, alpha, beta;
56
57 if ( ndaug == 1 )
58 {
59 p4[0].set( mass[0], 0.0, 0.0, 0.0 );
60 return 1.0;
61 }
62
63 if ( ndaug == 2 )
64 {
65
66 // Two body phase space
67
68 energy = ( mp * mp + mass[0] * mass[0] - mass[1] * mass[1] ) / ( 2.0 * mp );
69
70 p3 = sqrt( energy * energy - mass[0] * mass[0] );
71
72 p4[0].set( energy, 0.0, 0.0, p3 );
73
74 energy = mp - energy;
75 p3 = -1.0 * p3;
76 p4[1].set( energy, 0.0, 0.0, p3 );
77
78 // Now rotate four vectors.
79
80 alpha = EvtRandom::Flat( EvtConst::twoPi );
81 beta = acos( EvtRandom::Flat( -1.0, 1.0 ) );
82
83 p4[0].applyRotateEuler( alpha, beta, -alpha );
84 p4[1].applyRotateEuler( alpha, beta, -alpha );
85
86 return 1.0;
87 }
88
89 if ( ndaug != 2 )
90 {
91
92 double wtmax = 0.0;
93 double pm[5][30], to[4], pmin, pmax, psum, rnd[30];
94 double ran, wt, pa, costh, sinth, phi, p[4][30], be[4], bep, temp;
95 int i, il, ilr, i1, il1u, il1, il2r, ilu;
96 int il2 = 0;
97
98 for ( i = 0; i < ndaug; i++ )
99 {
100 pm[4][i] = 0.0;
101 rnd[i] = 0.0;
102 }
103
104 pm[0][0] = mp;
105 pm[1][0] = 0.0;
106 pm[2][0] = 0.0;
107 pm[3][0] = 0.0;
108 pm[4][0] = mp;
109
110 to[0] = mp;
111 to[1] = 0.0;
112 to[2] = 0.0;
113 to[3] = 0.0;
114
115 psum = 0.0;
116 for ( i = 1; i < ndaug + 1; i++ ) { psum = psum + mass[i - 1]; }
117
118 pm[4][ndaug - 1] = mass[ndaug - 1];
119
120 switch ( ndaug )
121 {
122
123 case 1: wtmax = 1.0 / 16.0; break;
124 case 2: wtmax = 1.0 / 150.0; break;
125 case 3: wtmax = 1.0 / 2.0; break;
126 case 4: wtmax = 1.0 / 5.0; break;
127 case 5: wtmax = 1.0 / 15.0; break;
128 case 6: wtmax = 1.0 / 15.0; break;
129 case 7: wtmax = 1.0 / 15.0; break;
130 case 8: wtmax = 1.0 / 15.0; break;
131 case 9: wtmax = 1.0 / 15.0; break;
132 case 10: wtmax = 1.0 / 15.0; break;
133 case 11: wtmax = 1.0 / 15.0; break;
134 case 12: wtmax = 1.0 / 15.0; break;
135 case 13: wtmax = 1.0 / 15.0; break;
136 case 14: wtmax = 1.0 / 15.0; break;
137 case 15: wtmax = 1.0 / 15.0; break;
138 default:
139 report( ERROR, "EvtGen" ) << "too many daughters for phase space..." << ndaug << " "
140 << mp << endl;
141 ;
142 break;
143 }
144
145 pmax = mp - psum + mass[ndaug - 1];
146
147 pmin = 0.0;
148
149 for ( ilr = 2; ilr < ndaug + 1; ilr++ )
150 {
151 il = ndaug + 1 - ilr;
152 pmax = pmax + mass[il - 1];
153 pmin = pmin + mass[il + 1 - 1];
154 wtmax = wtmax * EvtPawt( pmax, pmin, mass[il - 1] );
155 }
156
157 // report(INFO,"EvtGen") << "wtmax:"<<wtmax<<endl;
158
159 do {
160
161 rnd[0] = 1.0;
162 il1u = ndaug - 1;
163
164 for ( il1 = 2; il1 < il1u + 1; il1++ )
165 {
166 ran = EvtRandom::Flat();
167 for ( il2r = 2; il2r < il1 + 1; il2r++ )
168 {
169 il2 = il1 + 1 - il2r;
170 if ( ran <= rnd[il2 - 1] ) goto two39;
171 rnd[il2 + 1 - 1] = rnd[il2 - 1];
172 }
173 two39:
174 rnd[il2 + 1 - 1] = ran;
175 }
176
177 rnd[ndaug - 1] = 0.0;
178 wt = 1.0;
179 for ( ilr = 2; ilr < ndaug + 1; ilr++ )
180 {
181 il = ndaug + 1 - ilr;
182 pm[4][il - 1] = pm[4][il + 1 - 1] + mass[il - 1] +
183 ( rnd[il - 1] - rnd[il + 1 - 1] ) * ( mp - psum );
184 wt = wt * EvtPawt( pm[4][il - 1], pm[4][il + 1 - 1], mass[il - 1] );
185 }
186 if ( wt > wtmax )
187 {
188 report( ERROR, "EvtGen" )
189 << "wtmax to small in EvtPhaseSpace with " << ndaug << " daughters" << endl;
190 ;
191 }
192 } while ( wt < EvtRandom::Flat( wtmax ) );
193
194 ilu = ndaug - 1;
195
196 for ( il = 1; il < ilu + 1; il++ )
197 {
198 pa = EvtPawt( pm[4][il - 1], pm[4][il + 1 - 1], mass[il - 1] );
199 costh = EvtRandom::Flat( -1.0, 1.0 );
200 sinth = sqrt( 1.0 - costh * costh );
201 phi = EvtRandom::Flat( EvtConst::twoPi );
202 p[1][il - 1] = pa * sinth * cos( phi );
203 p[2][il - 1] = pa * sinth * sin( phi );
204 p[3][il - 1] = pa * costh;
205 pm[1][il + 1 - 1] = -p[1][il - 1];
206 pm[2][il + 1 - 1] = -p[2][il - 1];
207 pm[3][il + 1 - 1] = -p[3][il - 1];
208 p[0][il - 1] = sqrt( pa * pa + mass[il - 1] * mass[il - 1] );
209 pm[0][il + 1 - 1] = sqrt( pa * pa + pm[4][il + 1 - 1] * pm[4][il + 1 - 1] );
210 }
211
212 p[0][ndaug - 1] = pm[0][ndaug - 1];
213 p[1][ndaug - 1] = pm[1][ndaug - 1];
214 p[2][ndaug - 1] = pm[2][ndaug - 1];
215 p[3][ndaug - 1] = pm[3][ndaug - 1];
216
217 for ( ilr = 2; ilr < ndaug + 1; ilr++ )
218 {
219 il = ndaug + 1 - ilr;
220 be[0] = pm[0][il - 1] / pm[4][il - 1];
221 be[1] = pm[1][il - 1] / pm[4][il - 1];
222 be[2] = pm[2][il - 1] / pm[4][il - 1];
223 be[3] = pm[3][il - 1] / pm[4][il - 1];
224
225 for ( i1 = il; i1 < ndaug + 1; i1++ )
226 {
227 bep = be[1] * p[1][i1 - 1] + be[2] * p[2][i1 - 1] + be[3] * p[3][i1 - 1] +
228 be[0] * p[0][i1 - 1];
229 temp = ( p[0][i1 - 1] + bep ) / ( be[0] + 1.0 );
230 p[1][i1 - 1] = p[1][i1 - 1] + temp * be[1];
231 p[2][i1 - 1] = p[2][i1 - 1] + temp * be[2];
232 p[3][i1 - 1] = p[3][i1 - 1] + temp * be[3];
233 p[0][i1 - 1] = bep;
234 }
235 }
236
237 for ( ilr = 0; ilr < ndaug; ilr++ )
238 { p4[ilr].set( p[0][ilr], p[1][ilr], p[2][ilr], p[3][ilr] ); }
239
240 return 1.0;
241 }
242
243 return 1.0;
244}
mass
double mass
Definition CosmicGenerator.cxx:128
EvtPawt
double EvtPawt(double a, double b, double c)
Definition EvtGenKine.cc:32
alpha
double alpha
Definition EvtPhokharaDef.hh:17
mp
double mp
Definition EvtPhokharaDef.hh:17
report
ostream & report(Severity severity, const char *facility)
Definition EvtReport.cc:34
ERROR
@ ERROR
Definition EvtReport.hh:49
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
energy
************Class m_ypar INTEGER m_KeyWgt INTEGER m_KeyIHVP 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:28
EvtRandom::Flat
static double Flat()
Definition EvtRandom.cc:69
EvtVector4R::applyRotateEuler
void applyRotateEuler(double alpha, double beta, double gamma)
Definition EvtVector4R.cc:69
EvtVector4R::set
void set(int i, double d)
Definition EvtVector4R.hh:173

Referenced by EvtBtoXsEtap::decay(), EvtBtoXsgamma::decay(), EvtBtoXsll::decay(), and EvtParticle::initializePhaseSpace().

◆ PhaseSpacePole()

double EvtGenKine::PhaseSpacePole ( double M,
double m1,
double m2,
double m3,
double a,
EvtVector4R p4[10] )
static

Definition at line 246 of file EvtGenKine.cc.

251{
252
253 // f1 = 1 (phasespace)
254 // f2 = a*(1/m12sq)^2
255
256 double m12sqmax = ( M - m3 ) * ( M - m3 );
257 double m12sqmin = ( m1 + m2 ) * ( m1 + m2 );
258
259 double m13sqmax = ( M - m2 ) * ( M - m2 );
260 double m13sqmin = ( m1 + m3 ) * ( m1 + m3 );
261
262 double v1 = ( m12sqmax - m12sqmin ) * ( m13sqmax - m13sqmin );
263 double v2 = a * ( 1.0 / m12sqmin - 1.0 / m12sqmax ) * ( m13sqmax - m13sqmin );
264
265 double m12sq, m13sq;
266
267 double r = v1 / ( v1 + v2 );
268
269 // report(INFO,"EvtGen") << "v1,v2:"<<v1<<","<<v2<<endl;
270
271 double m13min, m13max;
272
273 do {
274
275 m13sq = EvtRandom::Flat( m13sqmin, m13sqmax );
276
277 if ( r > EvtRandom::Flat() ) { m12sq = EvtRandom::Flat( m12sqmin, m12sqmax ); }
278 else
279 {
280 m12sq =
281 1.0 / ( 1.0 / m12sqmin - EvtRandom::Flat() * ( 1.0 / m12sqmin - 1.0 / m12sqmax ) );
282 }
283
284 // kinematically allowed?
285 double E3star = ( M * M - m12sq - m3 * m3 ) / sqrt( 4 * m12sq );
286 double E1star = ( m12sq + m1 * m1 - m2 * m2 ) / sqrt( 4 * m12sq );
287 double p3star = sqrt( E3star * E3star - m3 * m3 );
288 double p1star = sqrt( E1star * E1star - m1 * m1 );
289 m13max =
290 ( E3star + E1star ) * ( E3star + E1star ) - ( p3star - p1star ) * ( p3star - p1star );
291 m13min =
292 ( E3star + E1star ) * ( E3star + E1star ) - ( p3star + p1star ) * ( p3star + p1star );
293
294 } while ( m13sq < m13min || m13sq > m13max );
295
296 double E2 = ( M * M + m2 * m2 - m13sq ) / ( 2.0 * M );
297 double E3 = ( M * M + m3 * m3 - m12sq ) / ( 2.0 * M );
298 double E1 = M - E2 - E3;
299 double p1mom = sqrt( E1 * E1 - m1 * m1 );
300 double p3mom = sqrt( E3 * E3 - m3 * m3 );
301 double cost13 = ( 2.0 * E1 * E3 + m1 * m1 + m3 * m3 - m13sq ) / ( 2.0 * p1mom * p3mom );
302
303 // report(INFO,"EvtGen") << m13sq << endl;
304 // report(INFO,"EvtGen") << m12sq << endl;
305 // report(INFO,"EvtGen") << E1 << endl;
306 // report(INFO,"EvtGen") << E2 << endl;
307 // report(INFO,"EvtGen") << E3 << endl;
308 // report(INFO,"EvtGen") << p1mom << endl;
309 // report(INFO,"EvtGen") << p3mom << endl;
310 // report(INFO,"EvtGen") << cost13 << endl;
311
312 p4[2].set( E3, 0.0, 0.0, p3mom );
313 p4[0].set( E1, p1mom * sqrt( 1.0 - cost13 * cost13 ), 0.0, p1mom * cost13 );
314 p4[1].set( E2, -p1mom * sqrt( 1.0 - cost13 * cost13 ), 0.0, -p1mom * cost13 - p3mom );
315
316 // report(INFO,"EvtGen") << "p4:"<<p4[0]<<p4[1]<<p4[2]<<endl;
317
318 double alpha = EvtRandom::Flat( EvtConst::twoPi );
319 double beta = acos( EvtRandom::Flat( -1.0, 1.0 ) );
320 double gamma = EvtRandom::Flat( EvtConst::twoPi );
321
322 p4[0].applyRotateEuler( alpha, beta, gamma );
323 p4[1].applyRotateEuler( alpha, beta, gamma );
324 p4[2].applyRotateEuler( alpha, beta, gamma );
325
326 return 1.0 + a / ( m12sq * m12sq );
327}
m2
double double * m2
Definition qcdloop1.h:83
m1
double * m1
Definition qcdloop1.h:83

Referenced by EvtParticle::initializePhaseSpace().

The documentation for this class was generated from the following files: