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 -#include - 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:"<