evtgen is hosted by Hepforge, IPPP Durham
EvtGen  2.0.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
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/EvtPatches.hh"
25 #include "EvtGenBase/EvtRandom.hh"
26 #include "EvtGenBase/EvtReport.hh"
27 #include "EvtGenBase/EvtVector4R.hh"
28 #include "EvtGenBase/EvtParticle.hh"
29 
30 #include <iostream>
31 #include <math.h>
32 using std::endl;
33 
34 double 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 
46 double 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 
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  //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  */
350 void 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 }
EvtGenKine::PhaseSpacePole
static double PhaseSpacePole(double M, double m1, double m2, double m3, double a, EvtVector4R p4[10])
Definition: EvtGenKine.cpp:263
EvtVector4R::applyRotateEuler
void applyRotateEuler(double alpha, double beta, double gamma)
Definition: EvtVector4R.cpp:83
EvtGenReport
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
Definition: EvtReport.cpp:33
EvtParticle::init
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
EvtPawt
double EvtPawt(double a, double b, double c)
Definition: EvtGenKine.cpp:34
EvtConst::twoPi
static const double twoPi
Definition: EvtConst.hh:27
EvtParticle::getId
EvtId getId() const
Definition: EvtParticle.cpp:130
EvtVector4R::set
void set(int i, double d)
Definition: EvtVector4R.hh:167
EvtRandom::Flat
static double Flat()
Definition: EvtRandom.cpp:72
EvtGenKine::ThreeBodyKine
static void ThreeBodyKine(const double m12Sq, const double m23Sq, EvtParticle *p)
Definition: EvtGenKine.cpp:350
EvtGenKine::PhaseSpace
static double PhaseSpace(int ndaug, double mass[30], EvtVector4R p4[30], double mp)
Definition: EvtGenKine.cpp:46
EvtParticle::mass
double mass() const
Definition: EvtParticle.cpp:165
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cpp:91
Generated by   doxygen 1.8.15