EvtGen 2.2.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
Loading...
Searching...
No Matches
EvtGenKine.cpp
Go to the documentation of this file.
1
2/***********************************************************************
3* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
4* *
5* This file is part of EvtGen. *
6* *
7* EvtGen is free software: you can redistribute it and/or modify *
8* it under the terms of the GNU General Public License as published by *
9* the Free Software Foundation, either version 3 of the License, or *
10* (at your option) any later version. *
11* *
12* EvtGen is distributed in the hope that it will be useful, *
13* but WITHOUT ANY WARRANTY; without even the implied warranty of *
14* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15* GNU General Public License for more details. *
16* *
17* You should have received a copy of the GNU General Public License *
18* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
19***********************************************************************/
20
21#include "EvtGenBase/EvtGenKine.hh"
22
23#include "EvtGenBase/EvtConst.hh"
24#include "EvtGenBase/EvtParticle.hh"
25#include "EvtGenBase/EvtRandom.hh"
26#include "EvtGenBase/EvtReport.hh"
27#include "EvtGenBase/EvtVector4R.hh"
28
29#include <cmath>
30#include <iostream>
31
32using std::endl;
33
34double EvtPawt( double a, double b, double c )
35{
36 double temp = ( a * a - ( b + c ) * ( b + c ) ) *
37 ( a * a - ( b - c ) * ( b - c ) );
38
39 if ( temp <= 0 ) {
40 return 0.0;
41 }
42
43 return sqrt( temp ) / ( 2.0 * a );
44}
45
46double EvtGenKine::PhaseSpace( int ndaug, double mass[30], EvtVector4R p4[30],
47 double mp )
48
49// N body phase space routine. Send parent with
50// daughters already defined ( Number and masses )
51// Returns four vectors in parent frame.
52
53{
54 double energy, p3, alpha, beta;
55
56 if ( ndaug == 1 ) {
57 p4[0].set( mass[0], 0.0, 0.0, 0.0 );
58 return 1.0;
59 }
60
61 if ( ndaug == 2 ) {
62 //Two body phase space
63
64 energy = ( mp * mp + mass[0] * mass[0] - mass[1] * mass[1] ) /
65 ( 2.0 * mp );
66
67 p3 = 0.0;
68 if ( energy > mass[0] ) {
69 p3 = sqrt( energy * energy - mass[0] * mass[0] );
70 }
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 double wtmax = 0.0;
91 double pm[5][30], pmin, pmax, psum, rnd[30];
92 double ran, wt, pa, costh, sinth, phi, p[4][30], be[4], bep, temp;
93 int i, il, ilr, i1, il1u, il1, il2r, ilu;
94 int il2 = 0;
95
96 for ( i = 0; i < ndaug; i++ ) {
97 pm[4][i] = 0.0;
98 rnd[i] = 0.0;
99 }
100
101 pm[0][0] = mp;
102 pm[1][0] = 0.0;
103 pm[2][0] = 0.0;
104 pm[3][0] = 0.0;
105 pm[4][0] = mp;
106
107 psum = 0.0;
108 for ( i = 1; i < ndaug + 1; i++ ) {
109 psum = psum + mass[i - 1];
110 }
111
112 pm[4][ndaug - 1] = mass[ndaug - 1];
113
114 switch ( ndaug ) {
115 case 1:
116 wtmax = 1.0 / 16.0;
117 break;
118 case 2:
119 wtmax = 1.0 / 150.0;
120 break;
121 case 3:
122 wtmax = 1.0 / 2.0;
123 break;
124 case 4:
125 wtmax = 1.0 / 5.0;
126 break;
127 case 5:
128 wtmax = 1.0 / 15.0;
129 break;
130 case 6:
131 wtmax = 1.0 / 15.0;
132 break;
133 case 7:
134 wtmax = 1.0 / 15.0;
135 break;
136 case 8:
137 wtmax = 1.0 / 15.0;
138 break;
139 case 9:
140 wtmax = 1.0 / 15.0;
141 break;
142 case 10:
143 wtmax = 1.0 / 15.0;
144 break;
145 case 11:
146 wtmax = 1.0 / 15.0;
147 break;
148 case 12:
149 wtmax = 1.0 / 15.0;
150 break;
151 case 13:
152 wtmax = 1.0 / 15.0;
153 break;
154 case 14:
155 wtmax = 1.0 / 15.0;
156 break;
157 case 15:
158 wtmax = 1.0 / 15.0;
159 break;
160 default:
161 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
162 << "too many daughters for phase space..." << ndaug << " "
163 << mp << endl;
164 ;
165 break;
166 }
167
168 pmax = mp - psum + mass[ndaug - 1];
169
170 pmin = 0.0;
171
172 for ( ilr = 2; ilr < ndaug + 1; ilr++ ) {
173 il = ndaug + 1 - ilr;
174 pmax = pmax + mass[il - 1];
175 pmin = pmin + mass[il + 1 - 1];
176 wtmax = wtmax * EvtPawt( pmax, pmin, mass[il - 1] );
177 }
178
179 do {
180 rnd[0] = 1.0;
181 il1u = ndaug - 1;
182
183 for ( il1 = 2; il1 < il1u + 1; il1++ ) {
184 ran = EvtRandom::Flat();
185 for ( il2r = 2; il2r < il1 + 1; il2r++ ) {
186 il2 = il1 + 1 - il2r;
187 if ( ran <= rnd[il2 - 1] )
188 goto two39;
189 rnd[il2 + 1 - 1] = rnd[il2 - 1];
190 }
191 two39:
192 rnd[il2 + 1 - 1] = ran;
193 }
194
195 rnd[ndaug - 1] = 0.0;
196 wt = 1.0;
197 for ( ilr = 2; ilr < ndaug + 1; ilr++ ) {
198 il = ndaug + 1 - ilr;
199 pm[4][il - 1] = pm[4][il + 1 - 1] + mass[il - 1] +
200 ( rnd[il - 1] - rnd[il + 1 - 1] ) * ( mp - psum );
201 wt = wt *
202 EvtPawt( pm[4][il - 1], pm[4][il + 1 - 1], mass[il - 1] );
203 }
204 if ( wt > wtmax ) {
205 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
206 << "wtmax to small in EvtPhaseSpace with " << ndaug
207 << " daughters" << endl;
208 ;
209 }
210 } while ( wt < EvtRandom::Flat( wtmax ) );
211
212 ilu = ndaug - 1;
213
214 for ( il = 1; il < ilu + 1; il++ ) {
215 pa = EvtPawt( pm[4][il - 1], pm[4][il + 1 - 1], mass[il - 1] );
216 costh = EvtRandom::Flat( -1.0, 1.0 );
217 sinth = sqrt( 1.0 - costh * costh );
218 phi = EvtRandom::Flat( EvtConst::twoPi );
219 p[1][il - 1] = pa * sinth * cos( phi );
220 p[2][il - 1] = pa * sinth * sin( phi );
221 p[3][il - 1] = pa * costh;
222 pm[1][il + 1 - 1] = -p[1][il - 1];
223 pm[2][il + 1 - 1] = -p[2][il - 1];
224 pm[3][il + 1 - 1] = -p[3][il - 1];
225 p[0][il - 1] = sqrt( pa * pa + mass[il - 1] * mass[il - 1] );
226 pm[0][il + 1 - 1] = sqrt( pa * pa +
227 pm[4][il + 1 - 1] * pm[4][il + 1 - 1] );
228 }
229
230 p[0][ndaug - 1] = pm[0][ndaug - 1];
231 p[1][ndaug - 1] = pm[1][ndaug - 1];
232 p[2][ndaug - 1] = pm[2][ndaug - 1];
233 p[3][ndaug - 1] = pm[3][ndaug - 1];
234
235 for ( ilr = 2; ilr < ndaug + 1; ilr++ ) {
236 il = ndaug + 1 - ilr;
237 be[0] = pm[0][il - 1] / pm[4][il - 1];
238 be[1] = pm[1][il - 1] / pm[4][il - 1];
239 be[2] = pm[2][il - 1] / pm[4][il - 1];
240 be[3] = pm[3][il - 1] / pm[4][il - 1];
241
242 for ( i1 = il; i1 < ndaug + 1; i1++ ) {
243 bep = be[1] * p[1][i1 - 1] + be[2] * p[2][i1 - 1] +
244 be[3] * p[3][i1 - 1] + be[0] * p[0][i1 - 1];
245 temp = ( p[0][i1 - 1] + bep ) / ( be[0] + 1.0 );
246 p[1][i1 - 1] = p[1][i1 - 1] + temp * be[1];
247 p[2][i1 - 1] = p[2][i1 - 1] + temp * be[2];
248 p[3][i1 - 1] = p[3][i1 - 1] + temp * be[3];
249 p[0][i1 - 1] = bep;
250 }
251 }
252
253 for ( ilr = 0; ilr < ndaug; ilr++ ) {
254 p4[ilr].set( p[0][ilr], p[1][ilr], p[2][ilr], p[3][ilr] );
255 }
256
257 return 1.0;
258 }
259
260 return 1.0;
261}
262
263double 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 //f1 = 1 (phasespace)
270 //f2 = a*(1/m12sq)^2
271
272 double m12sqmax = ( M - m3 ) * ( M - m3 );
273 double m12sqmin = ( m1 + m2 ) * ( m1 + m2 );
274
275 double m13sqmax = ( M - m2 ) * ( M - m2 );
276 double m13sqmin = ( m1 + m3 ) * ( m1 + m3 );
277
278 double v1 = ( m12sqmax - m12sqmin ) * ( m13sqmax - m13sqmin );
279 double v2 = a * ( 1.0 / m12sqmin - 1.0 / m12sqmax ) * ( m13sqmax - m13sqmin );
280
281 double m12sq, m13sq;
282
283 double r = v1 / ( v1 + v2 );
284
285 double m13min, m13max;
286
287 do {
288 m13sq = EvtRandom::Flat( m13sqmin, m13sqmax );
289
290 if ( r > EvtRandom::Flat() ) {
291 m12sq = EvtRandom::Flat( m12sqmin, m12sqmax );
292 } else {
293 m12sq = 1.0 /
294 ( 1.0 / m12sqmin -
295 EvtRandom::Flat() * ( 1.0 / m12sqmin - 1.0 / m12sqmax ) );
296 }
297
298 //kinematically allowed?
299 double E3star = ( M * M - m12sq - m3 * m3 ) / sqrt( 4 * m12sq );
300 double E1star = ( m12sq + m1 * m1 - m2 * m2 ) / sqrt( 4 * m12sq );
301 double p3star = sqrt( E3star * E3star - m3 * m3 );
302 double p1star = sqrt( E1star * E1star - m1 * m1 );
303 m13max = ( E3star + E1star ) * ( E3star + E1star ) -
304 ( p3star - p1star ) * ( p3star - p1star );
305 m13min = ( E3star + E1star ) * ( E3star + E1star ) -
306 ( p3star + p1star ) * ( p3star + p1star );
307
308 } while ( m13sq < m13min || m13sq > m13max );
309
310 double E2 = ( M * M + m2 * m2 - m13sq ) / ( 2.0 * M );
311 double E3 = ( M * M + m3 * m3 - m12sq ) / ( 2.0 * M );
312 double E1 = M - E2 - E3;
313 double p1mom = sqrt( E1 * E1 - m1 * m1 );
314 double p3mom = sqrt( E3 * E3 - m3 * m3 );
315 double cost13 = ( 2.0 * E1 * E3 + m1 * m1 + m3 * m3 - m13sq ) /
316 ( 2.0 * p1mom * p3mom );
317
318 //EvtGenReport(EVTGEN_INFO,"EvtGen") << m13sq << endl;
319 //EvtGenReport(EVTGEN_INFO,"EvtGen") << m12sq << endl;
320 //EvtGenReport(EVTGEN_INFO,"EvtGen") << E1 << endl;
321 //EvtGenReport(EVTGEN_INFO,"EvtGen") << E2 << endl;
322 //EvtGenReport(EVTGEN_INFO,"EvtGen") << E3 << endl;
323 //EvtGenReport(EVTGEN_INFO,"EvtGen") << p1mom << endl;
324 //EvtGenReport(EVTGEN_INFO,"EvtGen") << p3mom << endl;
325 //EvtGenReport(EVTGEN_INFO,"EvtGen") << cost13 << endl;
326
327 p4[2].set( E3, 0.0, 0.0, p3mom );
328 p4[0].set( E1, p1mom * sqrt( 1.0 - cost13 * cost13 ), 0.0, p1mom * cost13 );
329 p4[1].set( E2, -p1mom * sqrt( 1.0 - cost13 * cost13 ), 0.0,
330 -p1mom * cost13 - p3mom );
331
332 //EvtGenReport(EVTGEN_INFO,"EvtGen") << "p4:"<<p4[0]<<p4[1]<<p4[2]<<endl;
333
334 double alpha = EvtRandom::Flat( EvtConst::twoPi );
335 double beta = acos( EvtRandom::Flat( -1.0, 1.0 ) );
336 double gamma = EvtRandom::Flat( EvtConst::twoPi );
337
338 p4[0].applyRotateEuler( alpha, beta, gamma );
339 p4[1].applyRotateEuler( alpha, beta, gamma );
340 p4[2].applyRotateEuler( alpha, beta, gamma );
341
342 return 1.0 + a / ( m12sq * m12sq );
343}
344
345/*
346 * Function which takes two invariant masses squared in 3-body decay and
347 * parent after makeDaughters() and generateMassTree() and
348 * calculates/generates momenta of daughters and sets those.
349 */
350void EvtGenKine::ThreeBodyKine( const double m12Sq, const double m23Sq,
351 EvtParticle* p )
352{
353 const double mParent = p->mass();
354 EvtParticle* daug1 = p->getDaug( 0 );
355 EvtParticle* daug2 = p->getDaug( 1 );
356 EvtParticle* daug3 = p->getDaug( 2 );
357 const double mDaug1 = daug1->mass();
358 const double mDaug2 = daug2->mass();
359 const double mDaug3 = daug3->mass();
360 const double mParentSq{ mParent * mParent };
361 const double mDaug1Sq{ mDaug1 * mDaug1 };
362 const double mDaug2Sq{ mDaug2 * mDaug2 };
363 const double mDaug3Sq{ mDaug3 * mDaug3 };
364 const double invMParent{ 1. / mParent };
365
366 const double En1 = 0.5 * ( mParentSq + mDaug1Sq - m23Sq ) * invMParent;
367 const double En3 = 0.5 * ( mParentSq + mDaug3Sq - m12Sq ) * invMParent;
368 const double En2 = mParent - En1 - En3;
369 const double p1mag = std::sqrt( En1 * En1 - mDaug1Sq );
370 const double p2mag = std::sqrt( En2 * En2 - mDaug2Sq );
371 double cosPhi = 0.5 * ( mDaug1Sq + mDaug2Sq + 2 * En1 * En2 - m12Sq ) /
372 ( p1mag * p2mag );
373
374 double sinPhi = std::sqrt( 1 - cosPhi * cosPhi );
375 if ( EvtRandom::Flat( 0., 1. ) > 0.5 ) {
376 sinPhi *= -1;
377 }
378 const double p2x = p2mag * cosPhi;
379 const double p2y = p2mag * sinPhi;
380 const double p3x = -p1mag - p2x;
381 const double p3y = -p2y;
382
383 // Construct 4-momenta and rotate them randomly in space
384 EvtVector4R p1( En1, p1mag, 0., 0. );
385 EvtVector4R p2( En2, p2x, p2y, 0. );
386 EvtVector4R p3( En3, p3x, p3y, 0. );
387 const double euler1 = EvtRandom::Flat( 0., EvtConst::twoPi );
388 const double euler2 = std::acos( EvtRandom::Flat( -1.0, 1.0 ) );
389 const double euler3 = EvtRandom::Flat( 0., EvtConst::twoPi );
390 p1.applyRotateEuler( euler1, euler2, euler3 );
391 p2.applyRotateEuler( euler1, euler2, euler3 );
392 p3.applyRotateEuler( euler1, euler2, euler3 );
393
394 // set momenta for daughters
395 daug1->init( daug1->getId(), p1 );
396 daug2->init( daug2->getId(), p2 );
397 daug3->init( daug3->getId(), p3 );
398
399 return;
400}
EvtPawt
double EvtPawt(double a, double b, double c)
Definition EvtGenKine.cpp:34
EvtGenReport
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=nullptr)
Definition EvtReport.cpp:32
EVTGEN_ERROR
@ EVTGEN_ERROR
Definition EvtReport.hh:49
EvtConst::twoPi
static const double twoPi
Definition EvtConst.hh:27
EvtGenKine::PhaseSpacePole
static double PhaseSpacePole(double M, double m1, double m2, double m3, double a, EvtVector4R p4[10])
Definition EvtGenKine.cpp:263
EvtGenKine::PhaseSpace
static double PhaseSpace(int ndaug, double mass[30], EvtVector4R p4[30], double mp)
Definition EvtGenKine.cpp:46
EvtGenKine::ThreeBodyKine
static void ThreeBodyKine(const double m12Sq, const double m23Sq, EvtParticle *p)
Definition EvtGenKine.cpp:350
EvtParticle::init
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
EvtParticle::getId
EvtId getId() const
Definition EvtParticle.cpp:124
EvtParticle::getDaug
EvtParticle * getDaug(const int i)
Definition EvtParticle.hh:173
EvtParticle::mass
double mass() const
Definition EvtParticle.cpp:159
EvtRandom::Flat
static double Flat()
Definition EvtRandom.cpp:95
EvtVector4R::applyRotateEuler
void applyRotateEuler(double alpha, double beta, double gamma)
Definition EvtVector4R.cpp:83
EvtVector4R::set
void set(int i, double d)
Definition EvtVector4R.hh:168
Generated by doxygen 1.16.1