Page MenuHomeHEPForge
Files F19302687

D110.1759179618.diff
No OneTemporary
Actions

Size
31 KB
Referenced Files
None
Subscribers
None

D110.1759179618.diff
View Options

diff --git a/src/EvtGenBase/EvtGenKine.cpp b/src/EvtGenBase/EvtGenKine.cpp
--- a/src/EvtGenBase/EvtGenKine.cpp
+++ b/src/EvtGenBase/EvtGenKine.cpp
@@ -27,320 +27,780 @@
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtVector4R.hh"
-#include <cmath>
-#include <iostream>
-
using std::endl;
-double EvtPawt( double a, double b, double c )
+inline void orderdouble( double* x, int i, int j )
{
- double temp = ( a * a - ( b + c ) * ( b + c ) ) *
- ( a * a - ( b - c ) * ( b - c ) );
+ auto t = std::min( x[i], x[j] );
+ x[j] = std::max( x[i], x[j] );
+ x[i] = t;
+}
- if ( temp <= 0 ) {
- return 0.0;
- }
+#define Od( i, j ) orderdouble( x, i, j )
- return sqrt( temp ) / ( 2.0 * a );
+void sort15( double* x )
+{
+ Od( 0, 1 );
+ Od( 2, 3 );
+ Od( 4, 5 );
+ Od( 6, 7 );
+ Od( 8, 9 );
+ Od( 10, 11 );
+ Od( 12, 13 );
+ Od( 0, 2 );
+ Od( 4, 6 );
+ Od( 8, 10 );
+ Od( 12, 14 );
+ Od( 1, 3 );
+ Od( 5, 7 );
+ Od( 9, 11 );
+ Od( 0, 4 );
+ Od( 8, 12 );
+ Od( 1, 5 );
+ Od( 9, 13 );
+ Od( 2, 6 );
+ Od( 10, 14 );
+ Od( 3, 7 );
+ Od( 0, 8 );
+ Od( 1, 9 );
+ Od( 2, 10 );
+ Od( 3, 11 );
+ Od( 4, 12 );
+ Od( 5, 13 );
+ Od( 6, 14 );
+ Od( 5, 10 );
+ Od( 6, 9 );
+ Od( 3, 12 );
+ Od( 13, 14 );
+ Od( 7, 11 );
+ Od( 1, 2 );
+ Od( 4, 8 );
+ Od( 1, 4 );
+ Od( 7, 13 );
+ Od( 2, 8 );
+ Od( 11, 14 );
+ Od( 2, 4 );
+ Od( 5, 6 );
+ Od( 9, 10 );
+ Od( 11, 13 );
+ Od( 3, 8 );
+ Od( 7, 12 );
+ Od( 6, 8 );
+ Od( 10, 12 );
+ Od( 3, 5 );
+ Od( 7, 9 );
+ Od( 3, 4 );
+ Od( 5, 6 );
+ Od( 7, 8 );
+ Od( 9, 10 );
+ Od( 11, 12 );
+ Od( 6, 7 );
+ Od( 8, 9 );
}
-double EvtGenKine::PhaseSpace( int ndaug, double mass[30], EvtVector4R p4[30],
- double mp )
+void sort14( double* x )
+{
+ Od( 0, 1 );
+ Od( 2, 3 );
+ Od( 4, 5 );
+ Od( 6, 7 );
+ Od( 8, 9 );
+ Od( 10, 11 );
+ Od( 12, 13 );
+ Od( 0, 2 );
+ Od( 4, 6 );
+ Od( 8, 10 );
+ Od( 1, 3 );
+ Od( 5, 7 );
+ Od( 9, 11 );
+ Od( 0, 4 );
+ Od( 8, 12 );
+ Od( 1, 5 );
+ Od( 9, 13 );
+ Od( 2, 6 );
+ Od( 3, 7 );
+ Od( 0, 8 );
+ Od( 1, 9 );
+ Od( 2, 10 );
+ Od( 3, 11 );
+ Od( 4, 12 );
+ Od( 5, 13 );
+ Od( 5, 10 );
+ Od( 6, 9 );
+ Od( 3, 12 );
+ Od( 7, 11 );
+ Od( 1, 2 );
+ Od( 4, 8 );
+ Od( 1, 4 );
+ Od( 7, 13 );
+ Od( 2, 8 );
+ Od( 2, 4 );
+ Od( 5, 6 );
+ Od( 9, 10 );
+ Od( 11, 13 );
+ Od( 3, 8 );
+ Od( 7, 12 );
+ Od( 6, 8 );
+ Od( 10, 12 );
+ Od( 3, 5 );
+ Od( 7, 9 );
+ Od( 3, 4 );
+ Od( 5, 6 );
+ Od( 7, 8 );
+ Od( 9, 10 );
+ Od( 11, 12 );
+ Od( 6, 7 );
+ Od( 8, 9 );
+}
-// N body phase space routine. Send parent with
-// daughters already defined ( Number and masses )
-// Returns four vectors in parent frame.
+void sort13( double* x )
+{
+ Od( 1, 3 );
+ Od( 5, 11 );
+ Od( 4, 7 );
+ Od( 8, 9 );
+ Od( 0, 12 );
+ Od( 6, 10 );
+ Od( 1, 4 );
+ Od( 6, 8 );
+ Od( 9, 10 );
+ Od( 11, 12 );
+ Od( 0, 5 );
+ Od( 3, 7 );
+ Od( 1, 2 );
+ Od( 5, 9 );
+ Od( 10, 12 );
+ Od( 0, 6 );
+ Od( 8, 11 );
+ Od( 0, 1 );
+ Od( 4, 5 );
+ Od( 3, 8 );
+ Od( 2, 6 );
+ Od( 1, 2 );
+ Od( 3, 4 );
+ Od( 6, 7 );
+ Od( 8, 9 );
+ Od( 5, 11 );
+ Od( 1, 3 );
+ Od( 6, 10 );
+ Od( 2, 4 );
+ Od( 5, 8 );
+ Od( 9, 11 );
+ Od( 7, 12 );
+ Od( 2, 3 );
+ Od( 4, 6 );
+ Od( 7, 10 );
+ Od( 4, 5 );
+ Od( 7, 9 );
+ Od( 10, 11 );
+ Od( 6, 8 );
+ Od( 3, 4 );
+ Od( 5, 6 );
+ Od( 7, 8 );
+ Od( 9, 10 );
+ Od( 6, 7 );
+ Od( 8, 9 );
+}
+void sort12( double* x )
{
- double energy, p3, alpha, beta;
+ Od( 0, 1 );
+ Od( 2, 3 );
+ Od( 4, 5 );
+ Od( 6, 7 );
+ Od( 8, 9 );
+ Od( 10, 11 );
+ Od( 1, 3 );
+ Od( 5, 7 );
+ Od( 9, 11 );
+ Od( 0, 2 );
+ Od( 4, 6 );
+ Od( 8, 10 );
+ Od( 5, 9 );
+ Od( 6, 10 );
+ Od( 1, 2 );
+ Od( 7, 11 );
+ Od( 0, 4 );
+ Od( 1, 5 );
+ Od( 9, 10 );
+ Od( 2, 6 );
+ Od( 3, 7 );
+ Od( 4, 8 );
+ Od( 5, 9 );
+ Od( 1, 2 );
+ Od( 6, 10 );
+ Od( 0, 4 );
+ Od( 7, 11 );
+ Od( 3, 8 );
+ Od( 1, 4 );
+ Od( 7, 10 );
+ Od( 2, 3 );
+ Od( 5, 6 );
+ Od( 8, 9 );
+ Od( 3, 5 );
+ Od( 6, 8 );
+ Od( 2, 4 );
+ Od( 7, 9 );
+ Od( 3, 4 );
+ Od( 5, 6 );
+ Od( 7, 8 );
+}
- if ( ndaug == 1 ) {
- p4[0].set( mass[0], 0.0, 0.0, 0.0 );
- return 1.0;
- }
+void sort11( double* x )
+{
+ Od( 0, 9 );
+ Od( 1, 8 );
+ Od( 2, 7 );
+ Od( 3, 6 );
+ Od( 4, 5 );
+ Od( 0, 3 );
+ Od( 4, 10 );
+ Od( 1, 2 );
+ Od( 6, 9 );
+ Od( 7, 8 );
+ Od( 0, 1 );
+ Od( 2, 3 );
+ Od( 5, 8 );
+ Od( 9, 10 );
+ Od( 6, 7 );
+ Od( 1, 2 );
+ Od( 4, 6 );
+ Od( 8, 10 );
+ Od( 5, 9 );
+ Od( 0, 4 );
+ Od( 7, 8 );
+ Od( 1, 5 );
+ Od( 2, 9 );
+ Od( 3, 6 );
+ Od( 1, 4 );
+ Od( 5, 7 );
+ Od( 2, 3 );
+ Od( 6, 9 );
+ Od( 2, 4 );
+ Od( 6, 7 );
+ Od( 8, 9 );
+ Od( 3, 5 );
+ Od( 3, 4 );
+ Od( 5, 6 );
+ Od( 7, 8 );
+}
- if ( ndaug == 2 ) {
- //Two body phase space
+void sort10( double* x )
+{
+ Od( 1, 8 );
+ Od( 0, 4 );
+ Od( 5, 9 );
+ Od( 3, 7 );
+ Od( 2, 6 );
+ Od( 0, 3 );
+ Od( 6, 9 );
+ Od( 4, 7 );
+ Od( 0, 1 );
+ Od( 3, 6 );
+ Od( 8, 9 );
+ Od( 2, 5 );
+ Od( 1, 5 );
+ Od( 7, 9 );
+ Od( 0, 2 );
+ Od( 4, 8 );
+ Od( 1, 2 );
+ Od( 3, 4 );
+ Od( 5, 6 );
+ Od( 7, 8 );
+ Od( 1, 3 );
+ Od( 6, 8 );
+ Od( 2, 4 );
+ Od( 5, 7 );
+ Od( 2, 3 );
+ Od( 6, 7 );
+ Od( 4, 6 );
+ Od( 3, 5 );
+ Od( 4, 5 );
+}
- energy = ( mp * mp + mass[0] * mass[0] - mass[1] * mass[1] ) /
- ( 2.0 * mp );
+void sort9( double* x )
+{
+ Od( 1, 2 );
+ Od( 4, 5 );
+ Od( 7, 8 );
+ Od( 0, 1 );
+ Od( 3, 4 );
+ Od( 6, 7 );
+ Od( 1, 2 );
+ Od( 4, 5 );
+ Od( 7, 8 );
+ Od( 3, 6 );
+ Od( 0, 3 );
+ Od( 5, 8 );
+ Od( 4, 7 );
+ Od( 3, 6 );
+ Od( 2, 5 );
+ Od( 1, 4 );
+ Od( 1, 3 );
+ Od( 5, 8 );
+ Od( 4, 7 );
+ Od( 2, 6 );
+ Od( 5, 7 );
+ Od( 2, 3 );
+ Od( 4, 6 );
+ Od( 3, 4 );
+ Od( 5, 6 );
+}
- p3 = 0.0;
- if ( energy > mass[0] ) {
- p3 = sqrt( energy * energy - mass[0] * mass[0] );
- }
+void sort8( double* x )
+{
+ Od( 0, 1 );
+ Od( 2, 3 );
+ Od( 4, 5 );
+ Od( 6, 7 );
+ Od( 1, 3 );
+ Od( 5, 7 );
+ Od( 0, 2 );
+ Od( 4, 6 );
+ Od( 1, 2 );
+ Od( 5, 6 );
+ Od( 3, 7 );
+ Od( 2, 6 );
+ Od( 1, 5 );
+ Od( 0, 4 );
+ Od( 3, 5 );
+ Od( 2, 4 );
+ Od( 1, 2 );
+ Od( 3, 4 );
+ Od( 5, 6 );
+}
- p4[0].set( energy, 0.0, 0.0, p3 );
+void sort7( double* x )
+{
+ Od( 1, 2 );
+ Od( 3, 4 );
+ Od( 5, 6 );
+ Od( 0, 2 );
+ Od( 4, 6 );
+ Od( 3, 5 );
+ Od( 2, 6 );
+ Od( 1, 5 );
+ Od( 0, 4 );
+ Od( 2, 5 );
+ Od( 0, 3 );
+ Od( 2, 4 );
+ Od( 1, 3 );
+ Od( 0, 1 );
+ Od( 2, 3 );
+ Od( 4, 5 );
+}
- energy = mp - energy;
- p3 = -1.0 * p3;
- p4[1].set( energy, 0.0, 0.0, p3 );
+void sort6( double* x )
+{
+ Od( 0, 1 );
+ Od( 2, 3 );
+ Od( 4, 5 );
+ Od( 0, 2 );
+ Od( 3, 5 );
+ Od( 1, 4 );
+ Od( 0, 1 );
+ Od( 2, 3 );
+ Od( 4, 5 );
+ Od( 1, 2 );
+ Od( 3, 4 );
+ Od( 2, 3 );
+}
- //Now rotate four vectors.
+void sort5( double* x )
+{
+ Od( 1, 2 );
+ Od( 3, 4 );
+ Od( 1, 3 );
+ Od( 0, 2 );
+ Od( 2, 4 );
+ Od( 0, 3 );
+ Od( 0, 1 );
+ Od( 2, 3 );
+ Od( 1, 2 );
+}
- alpha = EvtRandom::Flat( EvtConst::twoPi );
- beta = acos( EvtRandom::Flat( -1.0, 1.0 ) );
+void sort4( double* x )
+{
+ Od( 0, 1 );
+ Od( 2, 3 );
+ Od( 0, 2 );
+ Od( 1, 3 );
+ Od( 1, 2 );
+}
- p4[0].applyRotateEuler( alpha, beta, -alpha );
- p4[1].applyRotateEuler( alpha, beta, -alpha );
+void sort3( double* x )
+{
+ Od( 0, 1 );
+ Od( 1, 2 );
+ Od( 0, 1 );
+}
- return 1.0;
- }
+void sort2( double* x )
+{
+ Od( 0, 1 );
+}
- if ( ndaug != 2 ) {
- double wtmax = 0.0;
- double pm[5][30], pmin, pmax, psum, rnd[30];
- double ran, wt, pa, costh, sinth, phi, p[4][30], be[4], bep, temp;
- int i, il, ilr, i1, il1u, il1, il2r, ilu;
- int il2 = 0;
+void sort1( double* )
+{
+}
- for ( i = 0; i < ndaug; i++ ) {
- pm[4][i] = 0.0;
- rnd[i] = 0.0;
- }
+typedef void( fun_t )( double* );
+fun_t* sortfuns[] = { sort1, sort2, sort3, sort4, sort5,
+ sort6, sort7, sort8, sort9, sort10,
+ sort11, sort12, sort13, sort14, sort15 };
- pm[0][0] = mp;
- pm[1][0] = 0.0;
- pm[2][0] = 0.0;
- pm[3][0] = 0.0;
- pm[4][0] = mp;
+inline double sqr( double x )
+{
+ return x * x;
+}
- psum = 0.0;
- for ( i = 1; i < ndaug + 1; i++ ) {
- psum = psum + mass[i - 1];
- }
+// return s = sin and c = cos of phi = k/2^32*2*M_PI
+// tested that |s|<=1 and |c|<=1
+// max |s^2 + c^2 - 1| <= 4.440892e-16
+// max relative difference |s-s_true|/s_true < 1.716228e-15
+void __attribute__( ( noinline ) )
+usincos( unsigned long kw, double& s, double& c )
+{
+ const static double st[] = { 0,
+ 4.90676743274180142550e-2,
+ 9.80171403295606019942e-2,
+ 1.46730474455361751659e-1,
+ 1.95090322016128267848e-1,
+ 2.42980179903263889948e-1,
+ 2.90284677254462367636e-1,
+ 3.36889853392220050689e-1,
+ 3.82683432365089771728e-1,
+ 4.27555093430282094321e-1,
+ 4.71396736825997648556e-1,
+ 5.14102744193221726594e-1,
+ 5.55570233019602224743e-1,
+ 5.95699304492433343467e-1,
+ 6.34393284163645498215e-1,
+ 6.71558954847018400625e-1,
+ 7.07106781186547524401e-1,
+ 7.40951125354959091176e-1,
+ 7.73010453362736960811e-1,
+ 8.03207531480644909807e-1,
+ 8.31469612302545237079e-1,
+ 8.57728610000272069902e-1,
+ 8.81921264348355029713e-1,
+ 9.03989293123443331586e-1,
+ 9.23879532511286756128e-1,
+ 9.41544065183020778413e-1,
+ 9.56940335732208864936e-1,
+ 9.70031253194543992604e-1,
+ 9.80785280403230449126e-1,
+ 9.89176509964780973452e-1,
+ 9.95184726672196886245e-1,
+ 9.98795456205172392715e-1,
+ 1 };
+ double x = (long)( kw << ( 5 + 2 ) ), x2 = x * x;
+ static const double as[3] = { 2.661032484442617284e-21,
+ -3.140503474026838861e-63,
+ 1.111886075860967104e-105 },
+ ac[3] = { -3.540546941629423600e-42,
+ 2.089245440416171469e-84,
+ -4.931274944723895814e-127 };
+ unsigned k = (unsigned long)kw >> 32;
+ int jk = ( k + ( 1 << 24 ) ) << 1, mask = jk >> 31,
+ absj = ( ( jk >> 26 ) ^ mask ) - mask;
+ double s0 = st[absj], c0 = st[32 - absj];
+ static const double sign[2] = { 1, -1 };
+ s0 *= sign[( k + 0 ) >> 31];
+ c0 *= sign[( k + ( 1 << 30 ) ) >> 31];
+ double sn = ( as[0] + x2 * ( as[1] + x2 * ( as[2] ) ) );
+ double dcs = x * ( ac[0] + x2 * ( ac[1] + x2 * ( ac[2] ) ) );
+ s = s0 + x * ( sn * c0 + dcs * s0 );
+ c = c0 + x * ( dcs * c0 - sn * s0 );
+}
- pm[4][ndaug - 1] = mass[ndaug - 1];
+// return s = sin and c = cos of phi = r*2*M_PI
+// 0 <= r < 1
+// minimal nonzero r is supposed to be 2^-53
+void rsincos( double r, double& s, double& c )
+{
+ long kw = r * 9007199254740992ul;
+ usincos( kw << 11, s, c );
+}
- switch ( ndaug ) {
- case 1:
- wtmax = 1.0 / 16.0;
- break;
- case 2:
- wtmax = 1.0 / 150.0;
- break;
- case 3:
- wtmax = 1.0 / 2.0;
- break;
- case 4:
- wtmax = 1.0 / 5.0;
- break;
- case 5:
- wtmax = 1.0 / 15.0;
- break;
- case 6:
- wtmax = 1.0 / 15.0;
- break;
- case 7:
- wtmax = 1.0 / 15.0;
- break;
- case 8:
- wtmax = 1.0 / 15.0;
- break;
- case 9:
- wtmax = 1.0 / 15.0;
- break;
- case 10:
- wtmax = 1.0 / 15.0;
- break;
- case 11:
- wtmax = 1.0 / 15.0;
- break;
- case 12:
- wtmax = 1.0 / 15.0;
- break;
- case 13:
- wtmax = 1.0 / 15.0;
- break;
- case 14:
- wtmax = 1.0 / 15.0;
- break;
- case 15:
- wtmax = 1.0 / 15.0;
- break;
- default:
- EvtGenReport( EVTGEN_ERROR, "EvtGen" )
- << "too many daughters for phase space..." << ndaug << " "
- << mp << endl;
- ;
- break;
- }
+double EvtGenKine::PhaseSpace( int ndaug, double mass[30],
+ EvtVector4R p4[30], double mp )
- pmax = mp - psum + mass[ndaug - 1];
+// N body phase space routine. Send parent with
+// daughters already defined ( Number and masses )
+// Returns four vectors in parent frame.
- pmin = 0.0;
+{
+ // Maximal weight scale -- based on numerical experiments
+ const static double wtscale[] = {
+ 1.000000e+00, 1.000000e+00, 5.000707e-01, 1.924501e-01, 6.249988e-02,
+ 1.789312e-02, 4.632982e-03, 1.104825e-03, 2.450442e-04, 5.095439e-05,
+ 1.022054e-05, 1.834890e-06, 3.216923e-07, 5.540728e-08, 9.468573e-09 };
- for ( ilr = 2; ilr < ndaug + 1; ilr++ ) {
- il = ndaug + 1 - ilr;
- pmax = pmax + mass[il - 1];
- pmin = pmin + mass[il + 1 - 1];
- wtmax = wtmax * EvtPawt( pmax, pmin, mass[il - 1] );
+ if ( ndaug == 1 ) {
+ p4[0].set( mass[0], 0.0, 0.0, 0.0 );
+ } else if ( ndaug == 2 ) {
+ //Two body phase space
+ double en = ( mp * mp + mass[0] * mass[0] - mass[1] * mass[1] ) /
+ ( 2.0 * mp );
+ double p3 = ( en > mass[0] ) ? sqrt( en * en - mass[0] * mass[0] ) : 0;
+
+ //Now uniformly distribute over sphere
+ double s, c;
+ rsincos( EvtRandom::random(), s, c );
+ double z = EvtRandom::random() * 2 - 1, r = sqrt( 1 - z * z );
+ double pt = p3 * r, px = pt * c, py = pt * s, pz = p3 * z;
+ p4[0].set( en, px, py, pz );
+ p4[1].set( mp - en, -px, -py, -pz );
+ } else if ( ndaug == 3 ) {
+#define ipck( i, j, k ) ( i | ( j << 2 ) | ( k << 4 ) )
+ const static unsigned char indx[9] = { ipck( 2, 1, 0 ),
+ ipck( 2, 0, 1 ),
+ 0,
+ ipck( 1, 2, 0 ),
+ 0,
+ ipck( 0, 2, 1 ),
+ 0,
+ ipck( 1, 0, 2 ),
+ ipck( 0, 1, 2 ) };
+#undef ipck
+ double m0 = mass[0], m1 = mass[1], m2 = mass[2];
+ unsigned i0 = ( m0 > m1 ) + ( m0 > m2 );
+ unsigned i1 = ( m1 >= m0 ) + ( m1 > m2 );
+ unsigned i2 = ( m2 >= m0 ) + ( m2 >= m1 );
+ unsigned I = i0 + 3 * ( i1 + 3 * i2 ), J = indx[( I - 5 ) >> 1],
+ j0 = J & 3, j1 = ( J >> 2 ) & 3, j2 = ( J >> 4 ) & 3;
+ double wtmax = wtscale[ndaug - 1] * wtscale[ndaug - 1];
+ auto order = []( double& a, double& b ) -> void {
+ auto t = std::max( a, b );
+ a = std::min( a, b );
+ b = t;
+ };
+ order( m0, m1 );
+ order( m1, m2 );
+ double u0 = m0 + m1, v0 = m0 - m1, M1max = mp - m2, M02 = sqr( M1max );
+ double u02 = u0 * u0, v02 = v0 * v0;
+ order( m0, m1 );
+ double u1 = u0 + m2, v1 = u0 - m2, M12 = sqr( mp ), dE = mp - u1;
+ if ( dE <= 0 ) {
+ printf( "Not enough energy: Mtot = %f Etot = %f\n", u1, mp );
+ exit( -1 );
}
+ wtmax /= M12 * M02;
+ wtmax *= ( M02 - u02 ) * ( M02 - v02 ) * ( M12 - u1 * u1 ) *
+ ( M12 - v1 * v1 );
+ double wt, wtd, R, M1, p0, p1;
do {
- rnd[0] = 1.0;
- il1u = ndaug - 1;
-
- for ( il1 = 2; il1 < il1u + 1; il1++ ) {
- ran = EvtRandom::Flat();
- for ( il2r = 2; il2r < il1 + 1; il2r++ ) {
- il2 = il1 + 1 - il2r;
- if ( ran <= rnd[il2 - 1] )
- goto two39;
- rnd[il2 + 1 - 1] = rnd[il2 - 1];
- }
- two39:
- rnd[il2 + 1 - 1] = ran;
- }
+ M1 = u0 + dE * EvtRandom::random();
+ R = wtmax * sqr( EvtRandom::random() );
+ double M02 = M1 * M1;
+ p0 = ( M02 - u02 ) * ( M02 - v02 );
+ double u1 = M1 + m2, v1 = M1 - m2;
+ p1 = ( M12 - u1 * u1 ) * ( M12 - v1 * v1 );
+ wt = p0 * p1;
+ wtd = M02 * M12;
+ } while ( wt < wtd * R );
+ double iM1 = 1 / M1;
+ p0 = sqrt( p0 ) * 0.5 * iM1;
+ p1 = sqrt( p1 ) / ( 2 * mp );
+
+ // for tree particle decay all momenta lie in a plane so let us
+ // generate momenta in the x-y plane with the third particle
+ // momentum along x-axis and then randomly rotate them
+ double E3 = sqrt( p1 * p1 + m2 * m2 );
+ double E12 = mp - E3;
+
+ // energies of particles 1 and 2 in their rest frame where they are back-to-back
+ double E1 = sqrt( p0 * p0 + m0 * m0 ), E2 = sqrt( p0 * p0 + m1 * m1 );
+
+ // projection on the x-y plane of uniform rotations in 3D
+ double cth = EvtRandom::random() * 2 - 1, sth = sqrt( 1 - cth * cth );
+ double px = p0 * cth, py = p0 * sth;
+
+ // boost 1 and 2 particle into the 1-2-3 rest frame
+ double g = E12 * iM1, gb = p1 * iM1;
+ double gpx = g * px, gbpx = gb * px;
+ double v0x = gpx - gb * E1, v0e = g * E1 - gbpx;
+ double v1x = -gpx - gb * E2, v1e = g * E2 + gbpx;
+ double v2x = p1;
+
+ double x, y;
+ rsincos( EvtRandom::random(), x, y );
+ double c, s;
+ rsincos( EvtRandom::random(), c, s );
+ double u = EvtRandom::random() * 2, z = u - 1, r = sqrt( 1 - z * z );
+ double Rx = s * y - c * x, Ry = c * y + s * x, ux = u * x, uy = u * y;
+ double R00 = ux * Rx + c, R01 = s - ux * Ry, R10 = uy * Rx - s,
+ R11 = c - uy * Ry, R20 = -r * Rx, R21 = r * Ry;
+ double pyx = R01 * py, pyy = R11 * py, pyz = R21 * py;
+
+ p4[j0].set( v0e, R00 * v0x + pyx, R10 * v0x + pyy, R20 * v0x + pyz );
+ p4[j1].set( v1e, R00 * v1x - pyx, R10 * v1x - pyy, R20 * v1x - pyz );
+ p4[j2].set( E3, R00 * v2x, R10 * v2x, R20 * v2x );
+ } else if ( ndaug < 16 ) {
+ const int nmax = 15;
+ double M[nmax], rndf[nmax], p[nmax - 1];
+ const double* m = mass;
+
+ double E0 = 0.0;
+ for ( int i = 0; i < ndaug; i++ ) {
+ E0 += m[i];
+ M[i] = m[i];
+ }
+ double dE = mp - E0;
+ if ( dE <= 0 ) {
+ printf( "Not enough energy: Mtot = %f Etot = %f\n", E0, mp );
+ exit( -1 );
+ }
+ rndf[0] = 0.0;
+ rndf[ndaug - 1] = 1;
+
+ double wtmax = wtscale[ndaug - 1] * wtscale[ndaug - 1];
+ sortfuns[ndaug - 1]( M );
+ m = M;
+ double Mmin = 0.0, Mmax = dE + m[0], wtmaxd = 1.0;
+ for ( int i = 0; i < ndaug - 1; i++ ) {
+ Mmin += m[i];
+ Mmax += m[i + 1];
+ double u = Mmin + m[i + 1], v = Mmin - m[i + 1], M2 = Mmax * Mmax;
+ wtmax *= ( M2 - u * u ) * ( M2 - v * v );
+ wtmaxd *= M2;
+ }
+ wtmax /= wtmaxd;
+ m = mass;
+ fun_t* sortfun = sortfuns[ndaug - 3];
+ double wt, wtd, R;
+ do {
+ for ( int i = 1; i < ndaug - 1; i++ )
+ rndf[i] = EvtRandom::random();
+ sortfun( rndf + 1 );
- rnd[ndaug - 1] = 0.0;
wt = 1.0;
- for ( ilr = 2; ilr < ndaug + 1; ilr++ ) {
- il = ndaug + 1 - ilr;
- pm[4][il - 1] = pm[4][il + 1 - 1] + mass[il - 1] +
- ( rnd[il - 1] - rnd[il + 1 - 1] ) * ( mp - psum );
- wt = wt *
- EvtPawt( pm[4][il - 1], pm[4][il + 1 - 1], mass[il - 1] );
- }
- if ( wt > wtmax ) {
- EvtGenReport( EVTGEN_ERROR, "EvtGen" )
- << "wtmax to small in EvtPhaseSpace with " << ndaug
- << " daughters" << endl;
- ;
- }
- } while ( wt < EvtRandom::Flat( wtmax ) );
-
- ilu = ndaug - 1;
-
- for ( il = 1; il < ilu + 1; il++ ) {
- pa = EvtPawt( pm[4][il - 1], pm[4][il + 1 - 1], mass[il - 1] );
- costh = EvtRandom::Flat( -1.0, 1.0 );
- sinth = sqrt( 1.0 - costh * costh );
- phi = EvtRandom::Flat( EvtConst::twoPi );
- p[1][il - 1] = pa * sinth * cos( phi );
- p[2][il - 1] = pa * sinth * sin( phi );
- p[3][il - 1] = pa * costh;
- pm[1][il + 1 - 1] = -p[1][il - 1];
- pm[2][il + 1 - 1] = -p[2][il - 1];
- pm[3][il + 1 - 1] = -p[3][il - 1];
- p[0][il - 1] = sqrt( pa * pa + mass[il - 1] * mass[il - 1] );
- pm[0][il + 1 - 1] = sqrt( pa * pa +
- pm[4][il + 1 - 1] * pm[4][il + 1 - 1] );
+ wtd = 1.0;
+ R = wtmax * sqr( EvtRandom::random() );
+ double M0 = m[0];
+ int i = 1;
+ do {
+ double f = rndf[i] - rndf[i - 1];
+ double m1 = m[i], M1 = m1 + M0 + f * dE;
+ double ma = M0 + m1, ms = M0 - m1, M2 = M1 * M1;
+ double t = ( M2 - ma * ma ) * ( M2 - ms * ms );
+ wt *= t;
+ wtd *= M2;
+ p[i - 1] = t;
+ M[i] = M1;
+ M0 = M1;
+ } while ( ++i < ndaug );
+ } while ( wt < wtd * R );
+
+ if ( wt > wtd * wtmax ) {
+ printf( "Warning: current weight is higher than supposed maximum: %e > %e\n",
+ sqrt( wt / wtd ), sqrt( wtmax ) );
}
- p[0][ndaug - 1] = pm[0][ndaug - 1];
- p[1][ndaug - 1] = pm[1][ndaug - 1];
- p[2][ndaug - 1] = pm[2][ndaug - 1];
- p[3][ndaug - 1] = pm[3][ndaug - 1];
-
- for ( ilr = 2; ilr < ndaug + 1; ilr++ ) {
- il = ndaug + 1 - ilr;
- be[0] = pm[0][il - 1] / pm[4][il - 1];
- be[1] = pm[1][il - 1] / pm[4][il - 1];
- be[2] = pm[2][il - 1] / pm[4][il - 1];
- be[3] = pm[3][il - 1] / pm[4][il - 1];
-
- for ( i1 = il; i1 < ndaug + 1; i1++ ) {
- bep = be[1] * p[1][i1 - 1] + be[2] * p[2][i1 - 1] +
- be[3] * p[3][i1 - 1] + be[0] * p[0][i1 - 1];
- temp = ( p[0][i1 - 1] + bep ) / ( be[0] + 1.0 );
- p[1][i1 - 1] = p[1][i1 - 1] + temp * be[1];
- p[2][i1 - 1] = p[2][i1 - 1] + temp * be[2];
- p[3][i1 - 1] = p[3][i1 - 1] + temp * be[3];
- p[0][i1 - 1] = bep;
- }
+ double* iM = rndf;
+ for ( int i = 0; i < ndaug - 1; i++ ) {
+ iM[i + 1] = 1 / M[i + 1];
+ p[i] = sqrt( p[i] ) * ( 0.5 * iM[i + 1] );
}
- for ( ilr = 0; ilr < ndaug; ilr++ ) {
- p4[ilr].set( p[0][ilr], p[1][ilr], p[2][ilr], p[3][ilr] );
- }
+ p4[0].set( sqrt( p[0] * p[0] + m[0] * m[0] ), 0, p[0], 0 );
+ int i = 1;
+ while ( 1 ) {
+ p4[i].set( sqrt( p[i - 1] * p[i - 1] + m[i] * m[i] ), 0, -p[i - 1],
+ 0 );
+
+ double cz = EvtRandom::random() * 2 - 1, sz = sqrt( 1 - cz * cz );
+ double sy, cy;
+ rsincos( EvtRandom::random(), sy, cy );
+ for ( int j = 0; j <= i; j++ ) {
+ double x = p4[j].get( 1 ), y = p4[j].get( 2 ), z = p4[j].get( 3 );
+ double xp = cz * x - sz * y,
+ yp = sz * x + cz * y; // rotation around z
+ double zp = sy * xp + cy * z;
+ xp = cy * xp - sy * z; // rotation around y
+ p4[j].set( 1, xp );
+ p4[j].set( 2, yp );
+ p4[j].set( 3, zp );
+ }
+
+ if ( i == ( ndaug - 1 ) )
+ break;
- return 1.0;
+ double E = sqrt( p[i] * p[i] + M[i] * M[i] ), gamma = E * iM[i],
+ betagamma = p[i] * iM[i];
+ for ( int j = 0; j <= i; j++ ) {
+ double e = p4[j].get( 0 ), py = p4[j].get( 2 );
+ p4[j].set( 0, gamma * e + betagamma * py );
+ p4[j].set( 2, gamma * py + betagamma * e );
+ }
+ i++;
+ }
+ } else {
+ printf( "No more than 15 particles! Ndaughter = %d", ndaug );
+ exit( -1 );
}
return 1.0;
}
-double EvtGenKine::PhaseSpacePole( double M, double m1, double m2, double m3,
- double a, EvtVector4R p4[10] )
+double PhaseSpacePole( double M, double m1, double m2, double m3, double a,
+ EvtVector4R p4[10] )
// generate kinematics for 3 body decays, pole for the m1,m2 mass.
-
{
- //f1 = 1 (phasespace)
- //f2 = a*(1/m12sq)^2
-
- double m12sqmax = ( M - m3 ) * ( M - m3 );
- double m12sqmin = ( m1 + m2 ) * ( m1 + m2 );
-
- double m13sqmax = ( M - m2 ) * ( M - m2 );
- double m13sqmin = ( m1 + m3 ) * ( m1 + m3 );
-
- double v1 = ( m12sqmax - m12sqmin ) * ( m13sqmax - m13sqmin );
- double v2 = a * ( 1.0 / m12sqmin - 1.0 / m12sqmax ) * ( m13sqmax - m13sqmin );
-
+ double m12sqmin = ( m1 + m2 ) * ( m1 + m2 ),
+ m12sqmax = ( M - m3 ) * ( M - m3 );
+ double m13sqmin = ( m1 + m3 ) * ( m1 + m3 ),
+ m13sqmax = ( M - m2 ) * ( M - m2 );
+ double d12 = m12sqmax - m12sqmin, d13 = m13sqmax - m13sqmin;
+ double M2 = M * M, m32 = m3 * m3, m12 = m1 * m1, m22 = m2 * m2;
+ double c0 = M2 - m32, c1 = m12 - m22, c2 = m32 + m12;
+ double ab12 = log( m12sqmax / m12sqmin ); // \int_m12sqmin^m12sqmax dx/x
+ double r = d12 / ( d12 + a * ab12 );
double m12sq, m13sq;
-
- double r = v1 / ( v1 + v2 );
-
- double m13min, m13max;
-
do {
- m13sq = EvtRandom::Flat( m13sqmin, m13sqmax );
-
- if ( r > EvtRandom::Flat() ) {
- m12sq = EvtRandom::Flat( m12sqmin, m12sqmax );
- } else {
- m12sq = 1.0 /
- ( 1.0 / m12sqmin -
- EvtRandom::Flat() * ( 1.0 / m12sqmin - 1.0 / m12sqmax ) );
- }
-
- //kinematically allowed?
- double E3star = ( M * M - m12sq - m3 * m3 ) / sqrt( 4 * m12sq );
- double E1star = ( m12sq + m1 * m1 - m2 * m2 ) / sqrt( 4 * m12sq );
- double p3star = sqrt( E3star * E3star - m3 * m3 );
- double p1star = sqrt( E1star * E1star - m1 * m1 );
- m13max = ( E3star + E1star ) * ( E3star + E1star ) -
- ( p3star - p1star ) * ( p3star - p1star );
- m13min = ( E3star + E1star ) * ( E3star + E1star ) -
- ( p3star + p1star ) * ( p3star + p1star );
-
- } while ( m13sq < m13min || m13sq > m13max );
-
- double E2 = ( M * M + m2 * m2 - m13sq ) / ( 2.0 * M );
- double E3 = ( M * M + m3 * m3 - m12sq ) / ( 2.0 * M );
+ double z0 = EvtRandom::random(), z1 = EvtRandom::random(),
+ z2 = EvtRandom::random();
+ m13sq = m13sqmin + z0 * d13;
+ m12sq = ( r > z1 ) ? m12sqmin + z2 * d12 : m12sqmin * exp( z2 * ab12 );
+ double E3s = c0 - m12sq, E1s = m12sq + c1;
+ double w = 2 * m12sq, e = 4 * m12sq;
+ double A = ( m13sq - c2 ) * w - E3s * E1s;
+ double B = ( E3s * E3s - e * m32 ) * ( E1s * E1s - e * m12 );
+ if ( A * A < B )
+ break;
+ } while ( true );
+ double iM = 0.5 / M;
+ double E2 = ( M2 + m22 - m13sq ) * iM;
+ double E3 = ( M2 + m32 - m12sq ) * iM;
double E1 = M - E2 - E3;
- double p1mom = sqrt( E1 * E1 - m1 * m1 );
- double p3mom = sqrt( E3 * E3 - m3 * m3 );
- double cost13 = ( 2.0 * E1 * E3 + m1 * m1 + m3 * m3 - m13sq ) /
- ( 2.0 * p1mom * p3mom );
-
- //EvtGenReport(EVTGEN_INFO,"EvtGen") << m13sq << endl;
- //EvtGenReport(EVTGEN_INFO,"EvtGen") << m12sq << endl;
- //EvtGenReport(EVTGEN_INFO,"EvtGen") << E1 << endl;
- //EvtGenReport(EVTGEN_INFO,"EvtGen") << E2 << endl;
- //EvtGenReport(EVTGEN_INFO,"EvtGen") << E3 << endl;
- //EvtGenReport(EVTGEN_INFO,"EvtGen") << p1mom << endl;
- //EvtGenReport(EVTGEN_INFO,"EvtGen") << p3mom << endl;
- //EvtGenReport(EVTGEN_INFO,"EvtGen") << cost13 << endl;
-
- p4[2].set( E3, 0.0, 0.0, p3mom );
- p4[0].set( E1, p1mom * sqrt( 1.0 - cost13 * cost13 ), 0.0, p1mom * cost13 );
- p4[1].set( E2, -p1mom * sqrt( 1.0 - cost13 * cost13 ), 0.0,
- -p1mom * cost13 - p3mom );
-
- //EvtGenReport(EVTGEN_INFO,"EvtGen") << "p4:"<<p4[0]<<p4[1]<<p4[2]<<endl;
-
- double alpha = EvtRandom::Flat( EvtConst::twoPi );
- double beta = acos( EvtRandom::Flat( -1.0, 1.0 ) );
- double gamma = EvtRandom::Flat( EvtConst::twoPi );
-
- p4[0].applyRotateEuler( alpha, beta, gamma );
- p4[1].applyRotateEuler( alpha, beta, gamma );
- p4[2].applyRotateEuler( alpha, beta, gamma );
-
- return 1.0 + a / ( m12sq * m12sq );
+ double p1 = sqrt( E1 * E1 - m12 );
+ double p3 = sqrt( E3 * E3 - m32 );
+ double cost13 = ( 2.0 * E1 * E3 + ( m12 + m32 - m13sq ) ) / ( 2.0 * p1 * p3 );
+
+ double px = p1 * cost13;
+ double v0x = px;
+ double v1x = -px - p3;
+ double py = p1 * sqrt( 1.0 - cost13 * cost13 );
+ double v2x = p3;
+
+ double x, y;
+ rsincos( EvtRandom::random(), x, y );
+ double c, s;
+ rsincos( EvtRandom::random(), c, s );
+ double u = EvtRandom::random() * 2, z = u - 1, t = sqrt( 1 - z * z );
+ double Rx = s * y - c * x, Ry = c * y + s * x, ux = u * x, uy = u * y;
+ double R00 = ux * Rx + c, R01 = s - ux * Ry, R10 = uy * Rx - s,
+ R11 = c - uy * Ry, R20 = -t * Rx, R21 = t * Ry;
+ double pyx = R01 * py, pyy = R11 * py, pyz = R21 * py;
+
+ p4[0].set( E1, R00 * v0x + pyx, R10 * v0x + pyy, R20 * v0x + pyz );
+ p4[1].set( E2, R00 * v1x - pyx, R10 * v1x - pyy, R20 * v1x - pyz );
+ p4[2].set( E3, R00 * v2x, R10 * v2x, R20 * v2x );
+
+ return 1.0 + a / m12sq;
}
/*

File Metadata

Mime Type
text/plain
Expires
Mon, Sep 29, 10:00 PM (16 h, 31 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6571736
Default Alt Text
D110.1759179618.diff (31 KB)

Event Timeline

HEPForge.org · .