9 #define M_PI 3.141592653589793238463
17 #if defined __cplusplus
101 xDataTOM_element *energyElement, *linearElement, *functional, *frameElement;
102 char const *nativeData;
103 double projectileMass_MeV, targetMass_MeV;
109 energy->
e_inCOMFactor = targetMass_MeV / ( projectileMass_MeV + targetMass_MeV );
112 energy->
type = energyType;
121 if( linearElement == NULL ) {
140 frameElement = functional; }
142 char const *toUnits[3] = {
"MeV",
"MeV",
"1/MeV" };
144 frameElement = linearElement;
167 if( strcmp( child->name,
"weighted" ) )
goto err;
185 char const *toUnits[2] = {
"MeV",
"" };
189 if( strcmp( child->
name,
"weight" ) == 0 ) {
191 else if( strcmp( child->
name,
"evaporation" ) == 0 ) {
215 char const *toUnits[2] = {
"MeV",
"MeV" };
227 if( std::fabs( 1. - norm ) > 0.001 )
printf(
"bad norm = %e\n", norm );
243 char const *U, *toUnits[2] = {
"MeV",
"MeV" };
264 char const *U, *toUnits[2] = {
"MeV",
"MeV" };
285 char const *U, *toUnits[2] = {
"MeV",
"MeV" };
297 toUnits[1] =
"1/MeV";
313 double E, T_M, EFL, EFH, argList[3], xs[] = { 1
e-5, 1
e-3, 1e-1, 1
e1, 1
e3, 1
e5, 3e7 },
norm;
322 char const *EF, *TMUnits[2] = {
"MeV",
"MeV" };
324 nXs =
sizeof( xs ) /
sizeof( xs[0] );
347 if( ( dists->
Ws = (
double *)
smr_malloc2( smr, length *
sizeof(
double ), 1,
"dists->Ws" ) ) == NULL )
goto err;
350 for( iE = 0; iE <
length; iE++ ) {
353 dist = &(dists->
dist[iE]);
357 (
void *) argList, 1
e-3, 0, 12, &status ) ) == NULL )
goto err;
366 if( ( dist->
Xs = (
double *)
smr_malloc2( smr, 3 * n *
sizeof(
double ), 0,
"dist->Xs" ) ) == NULL )
goto err;
368 dist->
pdf = &(dist->
Xs[
n]);
371 for( i1 = 0; i1 <
n; i1++ ) {
373 dist->
Xs[i1] = point->
x;
374 dist->
pdf[i1] = point->
y;
384 for( i1 = 0; i1 <
n; i1++ ) dist->
pdf[i1] /=
norm;
395 if( ptwXY_TM != NULL )
ptwXY_free( ptwXY_TM );
397 if( cdfX != NULL ) cdfX =
ptwX_free( cdfX );
406 double *
parameters = (
double *) argList, EFL, EFH, T_M;
422 double u1, u2, E1, E2, gamma1, gamma2, signG = 1;
424 u1 = std::sqrt( Ep ) - std::sqrt( E_F );
426 u2 = std::sqrt( Ep ) + std::sqrt( E_F );
431 if( *status !=
nfu_Okay )
return( 0. );
440 if( *status !=
nfu_Okay )
return( 0. );
441 return( ( u2 * std::sqrt( u2 ) * E2 - u1 * std::sqrt( u1 ) * E1 + signG * ( gamma2 - gamma1 ) ) / ( 3 * std::sqrt( E_F * T_M ) ) );
450 double xs[2] = { 0.0, 1.0 }, productMass_MeV,
norm;
488 int numberOfProducts = *((
int *) argList);
489 double e = 0.5 * ( 3 * numberOfProducts - 8 );
506 switch( energy->
type ) {
514 randomEp = decaySamplingInfo->
rng( decaySamplingInfo->
rngState );
518 decaySamplingInfo->
Ep = sampled.
x;
524 decaySamplingInfo->
Ep = theta * sampled.
x;
529 decaySamplingInfo->
Ep *=
theta;
534 decaySamplingInfo->
Ep *=
theta;
543 decaySamplingInfo->
Ep = sampled.
x;
564 double a = e_in_U_theta,
b,
c,
x, norm_a, xMin = 0., xMax =
a, sqrt_x, sqrt_pi_2 = std::sqrt(
M_PI ) / 2.;
566 sqrt_x = std::sqrt( a );
567 norm_a = sqrt_pi_2 * erf( sqrt_x ) - sqrt_x *
G4Exp( -a );
568 b = norm_a * decaySamplingInfo->
rng( decaySamplingInfo->
rngState );
569 for( i1 = 0; i1 < 16; i1++ ) {
570 x = 0.5 * ( xMin + xMax );
571 sqrt_x = std::sqrt( x );
572 c = sqrt_pi_2 * erf( sqrt_x ) - sqrt_x *
G4Exp( -x );
580 decaySamplingInfo->
Ep =
x;
590 double a = e_in_U_theta,
b,
c,
x, norm_a, xMin = 0., xMax =
a;
592 norm_a = 1 - ( 1 +
a ) *
G4Exp( -a );
593 b = 1. - norm_a * decaySamplingInfo->
rng( decaySamplingInfo->
rngState );
594 for( i1 = 0; i1 < 16; i1++ ) {
595 x = 0.5 * ( xMin + xMax );
596 c = ( 1 +
x ) *
G4Exp( -x );
604 decaySamplingInfo->
Ep =
x;
615 double WattMin = 0., WattMax = e_in_U,
x,
y,
z, energyOut,
rand1,
rand2;
617 x = 1. + ( Watt_b / ( 8. * Watt_a ) );
618 y = ( x + std::sqrt( x * x - 1. ) ) / Watt_a;
621 G4int icounter_max=1024;
625 if ( icounter > icounter_max ) {
626 G4cout <<
"Loop-counter exceeded the threshold value at " << __LINE__ <<
"th line of " << __FILE__ <<
"." <<
G4endl;
631 energyOut = y *
rand1;
633 while( ( ( rand2 - z * ( rand1 + 1. ) ) * ( rand2 - z * ( rand1 + 1. ) ) > Watt_b * y * rand1 ) || ( energyOut < WattMin ) || ( energyOut > WattMax ) );
634 decaySamplingInfo->
Ep = energyOut;
647 double rW = decaySamplingInfo->
rng( decaySamplingInfo->
rngState ), cumulativeW = 0.,
weight;
654 if( cumulativeW >= rW )
break;
659 #if defined __cplusplus