lars/geant4
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
minicalosim-spawning
geant4/source/processes/hadronic/models/lend/include/MCGIDI_headerSource.hpp
T
Gabriele Cosmo b4a16de652 Import Geant4 11.4.0 source tree
2025-12-05 08:54:02 +01:00

3449 lines
182 KiB
C++
Raw Permalink Blame History

/*
# <<BEGIN-copyright>>
# Copyright 2019, Lawrence Livermore National Security, LLC.
# This file is part of the gidiplus package (https://github.com/LLNL/gidiplus).
# gidiplus is licensed under the MIT license (see https://opensource.org/licenses/MIT).
# SPDX-License-Identifier: MIT
# <<END-copyright>>
*/
#ifndef MCGIDI_headerSource_hpp_included
#define MCGIDI_headerSource_hpp_included 1
#include <climits>
#ifdef MCGIDI_USE_DOUBLES
#define crossSectionSumError 1e-8
#else
#define crossSectionSumError 1e-6
#endif
// From file: MCGIDI_URR.cpp
/* *********************************************************************************************************//**
* Updates *this* if *a_protare* has a non-negative *URR_index*.
*
* @param a_protare [in] The protare whose *URR_index* is used to see if *this* needs updating.
* @param a_energy [in] The energy of the projectile.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::URR_protareInfos::updateProtare( MCGIDI::Protare const *a_protare, double a_energy, RNG && a_rng ) {
for( std::size_t i1 = 0; i1 < a_protare->numberOfProtares( ); ++i1 ) {
ProtareSingle *protareSingle = const_cast<ProtareSingle *>( a_protare->protare( i1 ) );
if( protareSingle->URR_index( ) >= 0 ) {
URR_protareInfo &URR_protare_info = m_URR_protareInfos[static_cast<std::size_t>(protareSingle->URR_index())];
URR_protare_info.m_inURR = protareSingle->inURR( a_energy );
if( URR_protare_info.inURR( ) ) URR_protare_info.m_rng_Value = a_rng( );
}
}
}
/* *********************************************************************************************************//**
* This function samples an energy and cosine of the angle for a photon for Klein Nishina scattering (i.e, incoherent photo-atomic scattering).
*
* @param a_energyIn [in] The energy of the incoming photon.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energyOut [in] The energy of the scattered photon.
* @param a_mu [in] The cosine of the angle of the scattered photon's z-axis and the incoming photon's z-axis.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI_sampleKleinNishina( double a_energyIn, RNG && a_rng, double *a_energyOut, double *a_mu ) {
/*
Description
Sample the Klein-Nishina distribution.
The unit of energy is the rest mass of the electron.
Reference: R. N. Blomquist and E. N. Gelbard, Nuclear Science
and Engineering, 83, 380-384 (1983)
This routine was taken from MCAPM which was from MCNP with only cosmetic changes.
Input
a_energyIn - incident photon energy ( in electron rest mass units )
*rng - user supplied random number generator
Output
*a_energyOut - exiting photon energy ( in electron rest mass units )
*a_mu - exiting photon cosine
*/
double a1, b1, t1, s1, r1, mu, energyOut;
a1 = 1.0 / a_energyIn;
b1 = 1.0 / ( 1.0 + 2.0 * a_energyIn );
if( a_energyIn < 3.0 ) { // Kahn''s method ( e < 1.5 MeV ) AECU-3259.
bool reject = true;
t1 = 1.0 / ( 1.0 + 8.0 * b1 );
do {
if( a_rng( ) <= t1 ) {
r1 = 2.0 * a_rng( );
s1 = 1.0 / ( 1.0 + a_energyIn * r1 );
mu = 1.0 - r1;
reject = a_rng( ) > 4.0 * s1 * ( 1.0 - s1 ); }
else {
s1 = ( 1.0 + 2.0 * a_energyIn * a_rng( ) ) * b1;
mu = 1.0 + a1 * ( 1.0 - 1.0 / s1 );
reject = a_rng( ) > 0.5 * ( mu * mu + s1 );
}
} while( reject );
energyOut = a_energyIn / ( 1 + a_energyIn * ( 1 - mu ) ); }
else { // Koblinger''s method ( e > 1.5 MeV ) NSE 56, 218 ( 1975 ).
t1 = a_rng( ) * ( 4.0 * a1 + 0.5 * ( 1.0 - b1 * b1 ) - ( 1.0 - 2.0 * ( 1.0 + a_energyIn ) * ( a1 * a1 ) ) * log( b1 ) );
if( t1 > 2.0 * a1 ) {
if( t1 > 4.0 * a1 ) {
if( t1 > 4.0 * a1 + 0.5 * ( 1.0 - b1 * b1 ) ) {
energyOut = a_energyIn * pow( b1, a_rng( ) );
mu = 1.0 + a1 - 1.0 / energyOut; }
else {
energyOut = a_energyIn * sqrt( 1.0 - a_rng( ) * ( 1.0 - b1 * b1 ) );
mu = 1.0 + a1 - 1.0 / energyOut;
} }
else {
energyOut = a_energyIn * ( 1.0 + a_rng( ) * ( b1 - 1.0 ) );
mu = 1.0 + a1 - 1.0 / energyOut; } }
else {
r1 = 2.0 * a_rng( );
mu = 1.0 - r1;
energyOut = 1.0 / ( a1 + r1 );
}
}
*a_mu = mu;
*a_energyOut = energyOut;
return;
}
/* *********************************************************************************************************//**
* This method samples the outgoing product data for the two outgoing particles in a two-body outgoing channel.
* First, is samples *mu*, the cosine of the product's outgoing angle, since this is for two-body interactions, *mu*
* is in the center-of-mass frame. It then calls kinetics_COMKineticEnergy2LabEnergyAndMomentum.
*
* @param a_X [in] The energy of the projectile in the lab frame.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::Distribution::sample( double a_X, MCGIDI::Sampling::Input &a_input, RNG && a_rng ) const {
switch( type( ) ) {
case Distributions::Type::none:
break;
case Distributions::Type::unspecified:
static_cast<Distributions::Unspecified const *>( this )->sample( a_X, a_input, a_rng );
break;
case Distributions::Type::angularTwoBody:
static_cast<Distributions::AngularTwoBody const *>( this )->sample( a_X, a_input, a_rng );
break;
case Distributions::Type::KalbachMann:
static_cast<Distributions::KalbachMann const *>( this )->sample( a_X, a_input, a_rng );
break;
case Distributions::Type::uncorrelated:
static_cast<Distributions::Uncorrelated const *>( this )->sample( a_X, a_input, a_rng );
break;
case Distributions::Type::branching3d:
static_cast<Distributions::Branching3d const *>( this )->sample( a_X, a_input, a_rng );
break;
case Distributions::Type::energyAngularMC:
static_cast<Distributions::EnergyAngularMC const *>( this )->sample( a_X, a_input, a_rng );
break;
case Distributions::Type::angularEnergyMC:
static_cast<Distributions::AngularEnergyMC const *>( this )->sample( a_X, a_input, a_rng );
break;
case Distributions::Type::coherentPhotoAtomicScattering:
static_cast<Distributions::CoherentPhotoAtomicScattering const *>( this )->sample( a_X, a_input, a_rng );
break;
case Distributions::Type::incoherentPhotoAtomicScattering:
static_cast<Distributions::IncoherentPhotoAtomicScattering const *>( this )->sample( a_X, a_input, a_rng );
break;
case Distributions::Type::incoherentBoundToFreePhotoAtomicScattering:
static_cast<Distributions::IncoherentBoundToFreePhotoAtomicScattering const *>( this )->sample( a_X, a_input, a_rng );
break;
case Distributions::Type::incoherentPhotoAtomicScatteringElectron:
static_cast<Distributions::IncoherentPhotoAtomicScatteringElectron const *>( this )->sample( a_X, a_input, a_rng );
break;
case Distributions::Type::pairProductionGamma:
static_cast<Distributions::PairProductionGamma const *>( this )->sample( a_X, a_input, a_rng );
break;
case Distributions::Type::coherentElasticTNSL:
static_cast<Distributions::CoherentElasticTNSL const *>( this )->sample( a_X, a_input, a_rng );
break;
case Distributions::Type::incoherentElasticTNSL:
static_cast<Distributions::IncoherentElasticTNSL const *>( this )->sample( a_X, a_input, a_rng );
break;
}
}
// From file: MCGIDI_distributions.cpp
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [out] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::Distribution::angleBiasing( Reaction const *a_reaction, double a_temperature, double a_energy_in, double a_mu_lab,
RNG && a_rng, double &a_energy_out ) const {
double probability = 0.0;
a_energy_out = 0.0;
switch( type( ) ) {
case Distributions::Type::none:
break;
case Distributions::Type::unspecified:
probability = static_cast<Distributions::Unspecified const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
case Distributions::Type::angularTwoBody:
probability = static_cast<Distributions::AngularTwoBody const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
case Distributions::Type::KalbachMann:
probability = static_cast<Distributions::KalbachMann const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
case Distributions::Type::uncorrelated:
probability = static_cast<Distributions::Uncorrelated const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
case Distributions::Type::branching3d:
probability = static_cast<Distributions::Branching3d const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
case Distributions::Type::energyAngularMC:
probability = static_cast<Distributions::EnergyAngularMC const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
case Distributions::Type::angularEnergyMC:
probability = static_cast<Distributions::AngularEnergyMC const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
case Distributions::Type::coherentPhotoAtomicScattering:
probability = static_cast<Distributions::CoherentPhotoAtomicScattering const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
case Distributions::Type::incoherentPhotoAtomicScattering:
probability = static_cast<Distributions::IncoherentPhotoAtomicScattering const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
case Distributions::Type::incoherentBoundToFreePhotoAtomicScattering:
probability = static_cast<Distributions::IncoherentBoundToFreePhotoAtomicScattering const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
case Distributions::Type::incoherentPhotoAtomicScatteringElectron:
probability = static_cast<Distributions::IncoherentPhotoAtomicScatteringElectron const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
case Distributions::Type::pairProductionGamma:
probability = static_cast<Distributions::PairProductionGamma const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
case Distributions::Type::coherentElasticTNSL:
probability = static_cast<Distributions::CoherentElasticTNSL const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
case Distributions::Type::incoherentElasticTNSL:
probability = static_cast<Distributions::IncoherentElasticTNSL const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
a_rng, a_energy_out );
break;
}
return( probability );
}
/* *********************************************************************************************************//**
* This function calculates the products outgoing data (i.e., energy, velocity/momentum) for the two products of
* a two-body interaction give the cosine of the first product's outgoing angle.
*
* @param a_beta [in] The velocity/speedOflight of the com frame relative to the lab frame.
* @param a_kinetic_com [in] Total kinetic energy (K1 + K2) in the COM frame.
* @param a_m3cc [in] The mass of the first product.
* @param a_m4cc [in] The mass of the second product.
* @param a_input [in] Sample options requested by user and where the products' outgoing data are returned.
***********************************************************************************************************/
inline LUPI_HOST_DEVICE void kinetics_COMKineticEnergy2LabEnergyAndMomentum( double a_beta, double a_kinetic_com,
double a_m3cc, double a_m4cc, MCGIDI::Sampling::Input &a_input ) {
/*
Relativity:
E = K + m, E^2 = K^2 + 2 K m + m^2, E^2 - m^2 = p^2 = K^2 + 2 K m
pc p v
---- = v, --- = --- = beta = b
E E c
K ( K + 2 m )
b^2 = ---------------
( K + m )^2
*/
double x, v_p, p, pp3, pp4, px3, py3, pz3, pz4, pz, p_perp2, E3, E4, gamma, m3cc2 = a_m3cc * a_m3cc, m4cc2 = a_m4cc * a_m4cc;
p = sqrt( a_kinetic_com * ( a_kinetic_com + 2. * a_m3cc ) * ( a_kinetic_com + 2. * a_m4cc ) *
( a_kinetic_com + 2. * ( a_m3cc + a_m4cc ) ) ) / ( 2. * ( a_kinetic_com + a_m3cc + a_m4cc ) );
py3 = p * sqrt( 1 - a_input.m_mu * a_input.m_mu );
px3 = py3 * cos( a_input.m_phi );
py3 *= sin( a_input.m_phi );
pz = p * a_input.m_mu;
if( 1 ) { // FIXME Assuming the answer is wanted in the lab frame for now.
a_input.m_frame = GIDI::Frame::lab;
E3 = sqrt( p * p + m3cc2 );
E4 = sqrt( p * p + m4cc2 );
gamma = sqrt( 1. / ( 1. - a_beta * a_beta ) );
pz3 = gamma * ( pz + a_beta * E3 );
pz4 = gamma * ( -pz + a_beta * E4 ); }
else { // COM frame.
a_input.m_frame = GIDI::Frame::centerOfMass;
pz3 = pz;
pz4 = -pz;
}
p_perp2 = px3 * px3 + py3 * py3;
a_input.m_px_vx1 = px3;
a_input.m_py_vy1 = py3;
a_input.m_pz_vz1 = pz3;
pp3 = p_perp2 + pz3 * pz3;
x = ( a_m3cc > 0 ) ? pp3 / ( 2 * m3cc2 ) : 1.;
if( x < 1e-5 ) {
a_input.m_energyOut1 = a_m3cc * x * ( 1 - 0.5 * x * ( 1 - x ) ); }
else {
a_input.m_energyOut1 = sqrt( m3cc2 + pp3 ) - a_m3cc;
}
a_input.m_px_vx2 = -px3;
a_input.m_py_vy2 = -py3;
a_input.m_pz_vz2 = pz4;
pp4 = p_perp2 + pz4 * pz4;
x = ( a_m4cc > 0 ) ? pp4 / ( 2 * m4cc2 ) : 1.;
if( x < 1e-5 ) {
a_input.m_energyOut2 = a_m4cc * x * ( 1 - 0.5 * x * ( 1 - x ) ); }
else {
a_input.m_energyOut2 = sqrt( m4cc2 + pp4 ) - a_m4cc;
}
if( a_input.wantVelocity( ) ) {
v_p = MCGIDI_speedOfLight_cm_sec / sqrt( pp3 + m3cc2 );
a_input.m_px_vx1 *= v_p;
a_input.m_py_vy1 *= v_p;
a_input.m_pz_vz1 *= v_p;
v_p = MCGIDI_speedOfLight_cm_sec / sqrt( pp4 + m4cc2 );
a_input.m_px_vx2 *= v_p;
a_input.m_py_vy2 *= v_p;
a_input.m_pz_vz2 *= v_p;
}
}
/* *********************************************************************************************************//**
* This method samples the outgoing product data for the two outgoing particles in a two-body outgoing channel.
* First, is samples *mu*, the cosine of the product's outgoing angle, since this is for two-body interactions, *mu*
* is in the center-of-mass frame. It then calls kinetics_COMKineticEnergy2LabEnergyAndMomentum.
*
* @param a_X [in] The energy of the projectile in the lab frame.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::AngularTwoBody::sample( double a_X, MCGIDI::Sampling::Input &a_input, RNG && a_rng ) const {
double initialMass = projectileMass( ) + targetMass( ), finalMass = productMass( ) + m_residualMass;
double beta = sqrt( a_X * ( a_X + 2. * projectileMass( ) ) ) / ( a_X + initialMass ); // beta = v/c.
double _x = targetMass( ) * ( a_X - m_twoBodyThreshold ) / ( finalMass * finalMass );
double Kp; // Kp is the total kinetic energy for m3 and m4 in the COM frame.
a_input.setSampledType( Sampling::SampledType::firstTwoBody );
if( m_Upscatter ) {
if( ( a_input.m_upscatterModel == Sampling::Upscatter::Model::B ) || ( a_input.m_upscatterModel == Sampling::Upscatter::Model::BSnLimits )
|| ( a_input.m_upscatterModel == Sampling::Upscatter::Model::DBRC ) ) {
if( upscatterModelB( a_X, a_input, a_rng ) ) return;
}
}
if( _x < 2e-5 ) {
Kp = finalMass * _x * ( 1 - 0.5 * _x * ( 1 - _x ) ); }
else { // This is the relativistic formula derived from E^2 - (pc)^2 is frame independent.
Kp = sqrt( finalMass * finalMass + 2 * targetMass( ) * ( a_X - m_twoBodyThreshold ) ) - finalMass;
}
if( Kp < 0 ) Kp = 0.; // FIXME There needs to be a better test here.
a_input.m_mu = m_angular->sample( a_X, a_rng( ), a_rng );
a_input.m_phi = 2. * M_PI * a_rng( );
kinetics_COMKineticEnergy2LabEnergyAndMomentum( beta, Kp, productMass( ), m_residualMass, a_input );
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [out] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::AngularTwoBody::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
double a_energy_in, double a_mu_lab, RNG && a_rng, double &a_energy_out ) const {
a_energy_out = 0.0;
double initialMass = projectileMass( ) + targetMass( ), finalMass = productMass( ) + m_residualMass;
double _x = targetMass( ) * ( a_energy_in - m_twoBodyThreshold ) / ( finalMass * finalMass );
double Kp; // Total kinetic energy of products in the center-of-mass.
if( _x < 2e-5 ) {
Kp = finalMass * _x * ( 1 - 0.5 * _x * ( 1 - _x ) ); }
else { // This is the relativistic formula derived from E^2 - (pc)^2 which is frame independent (i.e., an invariant).
Kp = sqrt( finalMass * finalMass + 2.0 * targetMass( ) * ( a_energy_in - m_twoBodyThreshold ) ) - finalMass;
}
if( Kp < 0 ) Kp = 0.; // FIXME There needs to be a better test here.
double energy_product_com = 0.5 * Kp * ( Kp + 2.0 * m_residualMass ) / ( Kp + productMass( ) + m_residualMass );
if( productMass( ) == 0.0 ) {
double boostBeta = sqrt( a_energy_in * ( a_energy_in + 2. * projectileMass( ) ) ) / ( a_energy_in + initialMass ); // Good, even for projectileMass = 0.
double one_mu_beta = 1.0 - a_mu_lab * boostBeta;
double mu_com = ( a_mu_lab - boostBeta ) / one_mu_beta;
double Jacobian = ( 1.0 - boostBeta * boostBeta ) / ( one_mu_beta * one_mu_beta );
a_energy_out = sqrt( 1.0 - boostBeta * boostBeta ) * energy_product_com * ( 1.0 + mu_com * boostBeta );
return( Jacobian * m_angular->evaluate( a_energy_in, mu_com ) );
}
double productBeta = MCGIDI_particleBeta( productMass( ), energy_product_com );
double boostBeta = sqrt( a_energy_in * ( a_energy_in + 2. * projectileMass( ) ) ) / ( a_energy_in + initialMass ); // beta = v/c.
double muPlus = 0.0, JacobianPlus = 0.0, muMinus = 0.0, JacobianMinus = 0.0;
int numberOfMus = muCOM_From_muLab( a_mu_lab, boostBeta, productBeta, muPlus, JacobianPlus, muMinus, JacobianMinus );
if( numberOfMus == 0 ) return( 0.0 );
double probability = JacobianPlus * m_angular->evaluate( a_energy_in, muPlus );
if( numberOfMus == 2 ) {
double probabilityMinus = JacobianMinus * m_angular->evaluate( a_energy_in, muMinus );
probability += probabilityMinus;
if( probabilityMinus > a_rng( ) * probability ) {
muPlus = muMinus;
}
}
double productBeta2 = productBeta * productBeta;
double productBetaLab2 = productBeta2 + boostBeta * boostBeta * ( 1.0 - productBeta2 * ( 1.0 - muPlus * muPlus ) ) + 2.0 * muPlus * productBeta * boostBeta;
productBetaLab2 /= 1.0 - muPlus * productBeta * boostBeta;
a_energy_out = MCGIDI::particleKineticEnergyFromBeta2( productMass( ), productBetaLab2 );
return( probability );
}
/* *********************************************************************************************************//**
* This method samples a targets velocity for elastic upscattering for upscatter model B and then calculates the outgoing
* product data for the projectile and target.
*
* @param a_kineticLab [in] The kinetic energy of the projectile in the lab frame.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE bool MCGIDI::Distributions::AngularTwoBody::upscatterModelB( double a_kineticLab, Sampling::Input &a_input, RNG && a_rng ) const {
const double Two_sqrtPi = 1.1283791670955125739;
const double C0 = 1.0410423479, C1 = 3.9626339162e-4, C2 =-1.8654539193e-3, C3 = 1.0264818153e-4;
double neutronMass = projectileMass( ); // Mass are in incident energy unit / c**2.
double _targetMass = targetMass( );
double temperature = a_input.temperature( );
double kineticLabMax = 1e4 * temperature;
if( a_input.m_upscatterModel == Sampling::Upscatter::Model::BSnLimits ) {
double kineticLabMax200 = 200.0 * temperature;
kineticLabMax = 1e3 * temperature * neutronMass / _targetMass;
if( kineticLabMax < kineticLabMax200 ) kineticLabMax = kineticLabMax200;
if( a_kineticLab >= 0.1 ) kineticLabMax = 0.9 * a_kineticLab; }
else {
if( kineticLabMax > 1e-2 ) {
kineticLabMax = 1e-2;
if( kineticLabMax < 100.0 * temperature ) {
kineticLabMax = 100.0 * temperature;
if( kineticLabMax > 10.0 ) kineticLabMax = 10.0; // Assumes energy is in MeV.
}
}
}
if( a_kineticLab > kineticLabMax ) return( false ); // Only for low neutron energy.
a_input.m_frame = GIDI::Frame::lab;
double muProjectileTarget, relativeBeta, targetBeta;
double targetThermalBeta = MCGIDI_particleBeta( _targetMass, temperature );
double neutronBeta = MCGIDI_particleBeta( neutronMass, a_kineticLab );
a_input.m_numberOfDBRC_rejections = 0;
bool continueLoop = false;
double crossSectionMax = 0.0;
if( a_input.m_upscatterModel == Sampling::Upscatter::Model::DBRC ) {
double targetThermalSpeed = m_modelDBRC_data->targetThermalSpeed( temperature ); // Non-relativistic calculations, unlike targetThermalBeta.
crossSectionMax = m_modelDBRC_data->crossSectionMax( a_kineticLab, targetThermalSpeed );
}
do {
continueLoop = false;
++a_input.m_numberOfDBRC_rejections;
do {
int MethodP1orP2 = 0; /* Assume P2 */
if( a_rng( ) * ( neutronBeta + Two_sqrtPi * targetThermalBeta ) < neutronBeta ) MethodP1orP2 = 1;
muProjectileTarget = 1.0 - 2.0 * a_rng( );
if( MethodP1orP2 == 0 ) { // x Exp( -x ) term.
targetBeta = targetThermalBeta * sqrt( -log( ( 1.0 - a_rng( ) ) * ( 1.0 - a_rng( ) ) ) ); }
else { // x^2 Exp( -x^2 ) term.
double x1;
do {
x1 = a_rng( );
x1 = sqrt( -log( ( 1.0 - a_rng( ) ) * ( 1.0 - x1 * x1 ) ) );
x1 = x1 / ( ( ( C3 * x1 + C2 ) * x1 + C1 ) * x1 + C0 );
} while( x1 > 4.0 );
targetBeta = targetThermalBeta * x1;
}
relativeBeta = sqrt( targetBeta * targetBeta + neutronBeta * neutronBeta - 2 * muProjectileTarget * targetBeta * neutronBeta );
} while( relativeBeta < ( targetBeta + neutronBeta ) * a_rng( ) );
if( a_input.m_upscatterModel == Sampling::Upscatter::Model::DBRC ) {
double relativeNeutronEnergy = 0.5 * productMass( ) * relativeBeta * relativeBeta;
if( m_modelDBRC_data->evaluate( relativeNeutronEnergy ) < a_rng( ) * crossSectionMax ) continueLoop = true;
}
} while( continueLoop );
double m1_12 = neutronMass / ( neutronMass + _targetMass );
double m2_12 = _targetMass / ( neutronMass + _targetMass );
double cosRelative = 0.0; // Cosine of angle between projectile velocity and relative velocity.
if( relativeBeta != 0.0 ) cosRelative = ( neutronBeta - muProjectileTarget * targetBeta ) / relativeBeta;
if( cosRelative > 1.0 ) {
cosRelative = 1.0; }
else if( cosRelative < -1.0 ) {
cosRelative = -1.0;
}
double sinRelative = sqrt( 1.0 - cosRelative * cosRelative ); // Sine of angle between projectile velocity and relative velocity.
a_input.m_muLab = muProjectileTarget;
a_input.m_targetBeta = targetBeta;
a_input.m_relativeBeta = relativeBeta;
double betaNeutronOut = m2_12 * relativeBeta;
double kineticEnergyRelative = particleKineticEnergy( neutronMass, betaNeutronOut );
double muCOM = m_angular->sample( kineticEnergyRelative, a_rng( ), a_rng );
double phiCOM = 2.0 * M_PI * a_rng( );
double SCcom = sqrt( 1.0 - muCOM * muCOM );
double SScom = SCcom * sin( phiCOM );
SCcom *= cos( phiCOM );
a_input.m_pz_vz1 = betaNeutronOut * ( muCOM * cosRelative - SCcom * sinRelative );
a_input.m_px_vx1 = betaNeutronOut * ( muCOM * sinRelative + SCcom * cosRelative );
a_input.m_py_vy1 = betaNeutronOut * SScom;
double massRatio = -neutronMass / _targetMass;
a_input.m_pz_vz2 = massRatio * a_input.m_pz_vz1;
a_input.m_px_vx2 = massRatio * a_input.m_px_vx1;
a_input.m_py_vy2 = massRatio * a_input.m_py_vy1;
double vCOMz = m1_12 * neutronBeta + m2_12 * muProjectileTarget * targetBeta; // Boost from center-of-mass to lab frame.
double vCOMx = m2_12 * sqrt( 1.0 - muProjectileTarget * muProjectileTarget ) * targetBeta;
a_input.m_pz_vz1 += vCOMz;
a_input.m_px_vx1 += vCOMx;
a_input.m_pz_vz2 += vCOMz;
a_input.m_px_vx2 += vCOMx;
double vx2_vy2 = a_input.m_px_vx1 * a_input.m_px_vx1 + a_input.m_py_vy1 * a_input.m_py_vy1;
double v2 = a_input.m_pz_vz1 * a_input.m_pz_vz1 + vx2_vy2;
a_input.m_mu = 0.0;
if( v2 != 0.0 ) a_input.m_mu = a_input.m_pz_vz1 / sqrt( v2 );
a_input.m_phi = atan2( a_input.m_py_vy1, a_input.m_px_vx1 );
a_input.m_energyOut1 = MCGIDI::particleKineticEnergyFromBeta2( neutronMass, v2 );
a_input.m_energyOut2 = MCGIDI::particleKineticEnergyFromBeta2( _targetMass, a_input.m_px_vx2 * a_input.m_px_vx2 + a_input.m_py_vy2 * a_input.m_py_vy2 + a_input.m_pz_vz2 * a_input.m_pz_vz2 );
a_input.m_px_vx1 *= MCGIDI_speedOfLight_cm_sec;
a_input.m_py_vy1 *= MCGIDI_speedOfLight_cm_sec;
a_input.m_pz_vz1 *= MCGIDI_speedOfLight_cm_sec;
a_input.m_px_vx2 *= MCGIDI_speedOfLight_cm_sec;
a_input.m_py_vy2 *= MCGIDI_speedOfLight_cm_sec;
a_input.m_pz_vz2 *= MCGIDI_speedOfLight_cm_sec;
if( !a_input.wantVelocity( ) ) { // Return momenta.
a_input.m_px_vx1 *= neutronMass; // Non-relativistic.
a_input.m_py_vy1 *= neutronMass;
a_input.m_pz_vz1 *= neutronMass;
a_input.m_px_vx2 *= _targetMass;
a_input.m_py_vy2 *= _targetMass;
a_input.m_pz_vz2 *= _targetMass;
}
double phi = 2.0 * M_PI * a_rng( );
double sine = sin( phi );
double cosine = cos( phi );
double saved = a_input.m_px_vx1;
a_input.m_px_vx1 = cosine * a_input.m_px_vx1 - sine * a_input.m_py_vy1;
a_input.m_py_vy1 = sine * saved + cosine * a_input.m_py_vy1;
return( true );
}
/* *********************************************************************************************************//**
* This method samples the outgoing product data by sampling the outgoing energy E' and mu from the uncorrelated
* E and mu probabilities. It also samples the outgoing phi uniformly between 0 and 2 pi.
*
* @param a_X [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::Uncorrelated::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
a_input.m_mu = m_angular->sample( a_X, a_rng( ), a_rng );
a_input.m_energyOut1 = m_energy->sample( a_X, a_rng( ), a_rng );
a_input.m_phi = 2. * M_PI * a_rng( );
a_input.m_frame = productFrame( );
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [in] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::Uncorrelated::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
double a_energy_in, double a_mu_lab, RNG && a_rng, double &a_energy_out ) const {
if( productFrame( ) != GIDI::Frame::lab ) {
a_energy_out = 0.0;
double initialMass = projectileMass( ) + targetMass( );
double boostBeta = sqrt( a_energy_in * ( a_energy_in + 2. * projectileMass( ) ) ) / ( a_energy_in + initialMass ); // Good, even for projectileMass = 0.
double energy_out_com = m_energy->sample( a_energy_in, a_rng( ), a_rng );
if( productMass( ) == 0.0 ) {
double one_mu_beta = 1.0 - a_mu_lab * boostBeta;
double mu_com = ( a_mu_lab - boostBeta ) / one_mu_beta;
double Jacobian = ( 1.0 - boostBeta * boostBeta ) / ( one_mu_beta * one_mu_beta );
a_energy_out = sqrt( 1.0 - boostBeta * boostBeta ) * energy_out_com * ( 1.0 + mu_com * boostBeta );
return( Jacobian * m_angular->evaluate( a_energy_in, mu_com ) );
}
double productBeta = MCGIDI_particleBeta( productMass( ), energy_out_com );
double muPlus = 0.0, JacobianPlus = 0.0, muMinus = 0.0, JacobianMinus = 0.0;
int numberOfMus = muCOM_From_muLab( a_mu_lab, boostBeta, productBeta, muPlus, JacobianPlus, muMinus, JacobianMinus );
if( numberOfMus == 0 ) return( 0.0 );
double probability = JacobianPlus * m_angular->evaluate( a_energy_in, muPlus );
if( numberOfMus == 2 ) {
double probabilityMinus = JacobianMinus * m_angular->evaluate( a_energy_in, muMinus );
probability += probabilityMinus;
if( probabilityMinus > a_rng( ) * probability ) muPlus = muMinus;
}
double productBeta2 = productBeta * productBeta;
double productBetaLab2 = productBeta2 + boostBeta * boostBeta * ( 1.0 - productBeta2 * ( 1.0 - muPlus * muPlus ) ) + 2.0 * muPlus * productBeta * boostBeta;
productBetaLab2 /= 1.0 - muPlus * productBeta * boostBeta;
a_energy_out = MCGIDI::particleKineticEnergyFromBeta2( productMass( ), productBetaLab2 );
return( probability );
}
a_energy_out = m_energy->sample( a_energy_in, a_rng( ), a_rng );
return( m_angular->evaluate( a_energy_in, a_mu_lab ) );
}
/* *********************************************************************************************************//**
* This method samples the outgoing branching photons.
*
* @param a_X [in] The energy of the projectile in the lab frame.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::Branching3d::sample( LUPI_maybeUnused double a_X, LUPI_maybeUnused Sampling::Input &a_input, LUPI_maybeUnused RNG && a_rng ) const {
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [in] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::Branching3d::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
LUPI_maybeUnused double a_energy_in, LUPI_maybeUnused double a_mu_lab, LUPI_maybeUnused RNG && a_rng, LUPI_maybeUnused double &a_energy_out ) const {
double probability = 0.0;
return( probability );
}
/* *********************************************************************************************************//**
* This method samples the outgoing product data by sampling the outgoing energy E' from the probability P(E'|E) and then samples mu from
* the probability P(mu|E,E'). It also samples the outgoing phi uniformly between 0 and 2 pi.
*
* @param a_X [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::EnergyAngularMC::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
double energyOut_1, energyOut_2;
a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
a_input.m_energyOut1 = m_energy->sample2dOf3d( a_X, a_rng( ), a_rng, &energyOut_1, &energyOut_2 );
a_input.m_mu = m_angularGivenEnergy->sample( a_X, energyOut_1, energyOut_2, a_rng( ), a_rng );
a_input.m_phi = 2. * M_PI * a_rng( );
a_input.m_frame = productFrame( );
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [in] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::EnergyAngularMC::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
double a_energy_in, double a_mu_lab, RNG && a_rng, double &a_energy_out ) const {
double probability = 0.0;
if( productFrame( ) == GIDI::Frame::centerOfMass ) {
a_energy_out = m_energy->sample( a_energy_in, a_rng( ), a_rng );
double initialMass = projectileMass( ) + targetMass( );
double boostBeta = sqrt( a_energy_in * ( a_energy_in + 2. * projectileMass( ) ) ) / ( a_energy_in + initialMass ); // Good, even for projectileMass = 0.
double energy_out_com = m_energy->sample( a_energy_in, a_rng( ), a_rng );
if( productMass( ) == 0.0 ) {
double one_mu_beta = 1.0 - a_mu_lab * boostBeta;
double mu_com = ( a_mu_lab - boostBeta ) / one_mu_beta;
double Jacobian = ( 1.0 - boostBeta * boostBeta ) / ( one_mu_beta * one_mu_beta );
a_energy_out = sqrt( 1.0 - boostBeta * boostBeta ) * energy_out_com * ( 1.0 + mu_com * boostBeta );
return( Jacobian * m_angularGivenEnergy->evaluate( a_energy_in, energy_out_com, mu_com ) );
}
double productBeta = MCGIDI_particleBeta( productMass( ), energy_out_com );
double muPlus = 0.0, JacobianPlus = 0.0, muMinus = 0.0, JacobianMinus = 0.0;
int numberOfMus = muCOM_From_muLab( a_mu_lab, boostBeta, productBeta, muPlus, JacobianPlus, muMinus, JacobianMinus );
if( numberOfMus == 0 ) return( 0.0 );
probability = JacobianPlus * m_angularGivenEnergy->evaluate( a_energy_in, energy_out_com, muPlus );
if( numberOfMus == 2 ) {
double probabilityMinus = JacobianMinus * m_angularGivenEnergy->evaluate( a_energy_in, energy_out_com, muMinus );
probability += probabilityMinus;
if( probabilityMinus > a_rng( ) * probability ) muPlus = muMinus;
}
double productBeta2 = productBeta * productBeta;
double productBetaLab2 = productBeta2 + boostBeta * boostBeta * ( 1.0 - productBeta2 * ( 1.0 - muPlus * muPlus ) ) + 2.0 * muPlus * productBeta * boostBeta;
productBetaLab2 /= 1.0 - muPlus * productBeta * boostBeta;
a_energy_out = MCGIDI::particleKineticEnergyFromBeta2( productMass( ), productBetaLab2 ); }
else {
a_energy_out = m_energy->sample( a_energy_in, a_rng( ), a_rng );
probability = m_angularGivenEnergy->evaluate( a_energy_in, a_energy_out, a_mu_lab );
}
return( probability );
}
/* *********************************************************************************************************//**
* This method samples the outgoing product data by sampling the outgoing mu from the probability P(mu|E) and then samples E' from
* the probability P(E'|E,mu). It also samples the outgoing phi uniformly between 0 and 2 pi.
*
* @param a_X [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::AngularEnergyMC::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
double mu_1, mu_2;
a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
a_input.m_mu = m_angular->sample2dOf3d( a_X, a_rng( ), a_rng, &mu_1, &mu_2 );
a_input.m_energyOut1 = m_energyGivenAngular->sample( a_X, mu_1, mu_2, a_rng( ), a_rng );
a_input.m_phi = 2. * M_PI * a_rng( );
a_input.m_frame = productFrame( );
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [in] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::AngularEnergyMC::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
double a_energy_in, double a_mu_lab, RNG && a_rng, double &a_energy_out ) const {
if( productFrame( ) != GIDI::Frame::lab ) LUPI_THROW( "AngularEnergyMC::angleBiasing: center-of-mass not supported." );
a_energy_out = m_energyGivenAngular->sample( a_energy_in, a_mu_lab, a_mu_lab, a_rng( ), a_rng );
return( m_angular->evaluate( a_energy_in, a_mu_lab ) );
}
/* *********************************************************************************************************//**
* This method samples the outgoing product data using the Kalbach-Mann formalism.
*
* @param a_X [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::KalbachMann::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
a_input.m_energyOut1 = m_f->sample( a_X, a_rng( ), a_rng );
double rValue = m_r->evaluate( a_X, a_input.m_energyOut1 );
double aValue = m_a->evaluate( a_X, a_input.m_energyOut1 );
// In the following: Cosh[ a mu ] + r Sinh[ a mu ] = ( 1 - r ) Cosh[ a mu ] + r ( Cosh[ a mu ] + Sinh[ a mu ] ).
if( a_rng( ) >= rValue ) { // Sample the '( 1 - r ) Cosh[ a mu ]' term.
double T = ( 2. * a_rng( ) - 1. ) * sinh( aValue );
a_input.m_mu = log( T + sqrt( T * T + 1. ) ) / aValue; }
else { // Sample the 'r ( Cosh[ a mu ] + Sinh[ a mu ] )' term.
double rng1 = a_rng( ), exp_a = exp( aValue );
a_input.m_mu = log( rng1 * exp_a + ( 1. - rng1 ) / exp_a ) / aValue;
}
if( a_input.m_mu < -1 ) a_input.m_mu = -1;
if( a_input.m_mu > 1 ) a_input.m_mu = 1;
a_input.m_phi = 2. * M_PI * a_rng( );
a_input.m_frame = productFrame( );
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [in] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::KalbachMann::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
double a_energy_in, double a_mu_lab, RNG && a_rng, double &a_energy_out ) const {
a_energy_out = 0.0;
double initialMass = projectileMass( ) + targetMass( );
double energy_out_com = m_f->sample( a_energy_in, a_rng( ), a_rng );
double productBeta = MCGIDI_particleBeta( productMass( ), energy_out_com );
double boostBeta = sqrt( a_energy_in * ( a_energy_in + 2. * projectileMass( ) ) ) / ( a_energy_in + initialMass ); // beta = v/c.
double muPlus = 0.0, JacobianPlus = 0.0, muMinus = 0.0, JacobianMinus = 0.0;
int numberOfMus = muCOM_From_muLab( a_mu_lab, boostBeta, productBeta, muPlus, JacobianPlus, muMinus, JacobianMinus );
if( numberOfMus == 0 ) return( 0.0 );
double rAtEnergyEnergyPrime = m_r->evaluate( a_energy_in, energy_out_com );
double aAtEnergyEnergyPrime = m_a->evaluate( a_energy_in, energy_out_com );
double aMu = aAtEnergyEnergyPrime * muPlus;
double probability = 0.5 * JacobianPlus;
if( productMass( ) == 0.0 ) {
probability *= 1.0 - rAtEnergyEnergyPrime + rAtEnergyEnergyPrime * aAtEnergyEnergyPrime * exp( aMu ) / sinh( aAtEnergyEnergyPrime ); }
else {
probability *= aAtEnergyEnergyPrime * ( cosh( aMu ) + rAtEnergyEnergyPrime * cosh( aMu ) ) / sinh( aAtEnergyEnergyPrime );
}
if( numberOfMus == 2 ) {
aMu = aAtEnergyEnergyPrime * muMinus;
double probabilityMinus = 0.5 * JacobianMinus;
if( productMass( ) == 0.0 ) {
probabilityMinus *= 1.0 - rAtEnergyEnergyPrime + rAtEnergyEnergyPrime * aAtEnergyEnergyPrime * exp( aMu ) / sinh( aAtEnergyEnergyPrime ); }
else {
probabilityMinus *= aAtEnergyEnergyPrime * ( cosh( aMu ) + rAtEnergyEnergyPrime * cosh( aMu ) ) / sinh( aAtEnergyEnergyPrime );
}
probability += probabilityMinus;
if( probabilityMinus > a_rng( ) * probability ) muPlus = muMinus;
}
double productBeta2 = productBeta * productBeta;
double productBetaLab2 = productBeta2 + boostBeta * boostBeta * ( 1.0 - productBeta2 * ( 1.0 - muPlus * muPlus ) ) + 2.0 * muPlus * productBeta * boostBeta;
productBetaLab2 /= 1.0 - muPlus * productBeta * boostBeta;
a_energy_out = MCGIDI::particleKineticEnergyFromBeta2( productMass( ), productBetaLab2 );
return( probability );
}
/* *********************************************************************************************************//**
* This method samples the outgoing product data from the coherent photo-atomic scattering law.
* It also samples the outgoing phi uniformly between 0 and 2 pi.
*
* @param a_X [in] The energy of the projectile in the lab frame.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::CoherentPhotoAtomicScattering::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
a_input.m_energyOut1 = a_X;
int intLowerIndex = binarySearchVector( a_X, m_energies );
if( intLowerIndex < 1 ) {
do {
a_input.m_mu = 1.0 - 2.0 * a_rng( );
} while( ( 1.0 + a_input.m_mu * a_input.m_mu ) < 2.0 * a_rng( ) ); }
else {
std::size_t lowerIndex = static_cast<std::size_t>( intLowerIndex );
double _a = m_a[lowerIndex];
double X_i = m_energies[lowerIndex];
double formFactor_i = m_formFactor[lowerIndex];
double formFactor_X_i = formFactor_i * X_i;
double Z = a_X / X_i;
double realAnomalousFactor = 0.0;
double imaginaryAnomalousFactor = 0.0;
if( m_anomalousDataPresent ) {
realAnomalousFactor = m_realAnomalousFactor->evaluate( a_X );
imaginaryAnomalousFactor = m_imaginaryAnomalousFactor->evaluate( a_X );
}
double anomalousFactorSquared = realAnomalousFactor * realAnomalousFactor + imaginaryAnomalousFactor * imaginaryAnomalousFactor;
double normalization = m_integratedFormFactorSquared[lowerIndex] + formFactor_X_i * formFactor_X_i * Z_a( Z, 2.0 * _a + 2.0 );
anomalousFactorSquared = 0.0;
if( anomalousFactorSquared != 0.0 ) {
double integratedFormFactor_i = m_integratedFormFactor[lowerIndex] + formFactor_X_i * Z_a( Z, _a + 2.0 );
normalization += 2.0 * integratedFormFactor_i * realAnomalousFactor + 0.5 * anomalousFactorSquared * a_X * a_X;
}
do {
double partialIntegral = a_rng( ) * normalization;
double X;
if( anomalousFactorSquared == 0.0 ) {
intLowerIndex = binarySearchVector( partialIntegral, m_integratedFormFactorSquared );
lowerIndex = static_cast<std::size_t>( intLowerIndex );
if( lowerIndex == 0 ) {
X = sqrt( 2.0 * partialIntegral ) / m_formFactor[0]; }
else {
double remainer = partialIntegral - m_integratedFormFactorSquared[lowerIndex];
double epsilon = 2.0 * m_a[lowerIndex] + 2.0;
X_i = m_energies[lowerIndex];
formFactor_i = m_formFactor[lowerIndex];
formFactor_X_i = formFactor_i * X_i;
remainer /= formFactor_X_i * formFactor_X_i;
if( fabs( epsilon ) < 1e-6 ) {
X = X_i * exp( remainer ); }
else {
X = X_i * pow( 1.0 + epsilon * remainer, 1.0 / epsilon );
}
} }
else { // Currently not implemented.
X = 0.5 * a_X;
}
double X_E = X / a_X;
a_input.m_mu = 1.0 - 2.0 * X_E * X_E;
} while( ( 1.0 + a_input.m_mu * a_input.m_mu ) < 2.0 * a_rng( ) );
}
a_input.m_phi = 2.0 * M_PI * a_rng( );
a_input.m_frame = productFrame( );
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [in] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::CoherentPhotoAtomicScattering::angleBiasing( Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
double a_energy_in, double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
a_energy_out = a_energy_in;
URR_protareInfos URR_protareInfos1;
double sigma = a_reaction->protareSingle( )->reactionCrossSection( a_reaction->reactionIndex( ), URR_protareInfos1, 0.0, a_energy_in );
double formFactor = evaluateFormFactor( a_energy_in, a_mu_lab );
double imaginaryAnomalousFactor = 0.0;
if( m_anomalousDataPresent ) {
formFactor += m_realAnomalousFactor->evaluate( a_energy_in );
imaginaryAnomalousFactor = m_imaginaryAnomalousFactor->evaluate( a_energy_in );
}
double probability = M_PI * MCGIDI_classicalElectronRadius * MCGIDI_classicalElectronRadius * ( 1.0 + a_mu_lab * a_mu_lab )
* ( formFactor * formFactor + imaginaryAnomalousFactor * imaginaryAnomalousFactor ) / sigma;
return( probability );
}
/* *********************************************************************************************************//**
* This method samples the outgoing product data by sampling the outgoing energy E' from the probability P(E'|E) and then samples mu from
* the probability P(mu|E,E'). It also samples the outgoing phi uniformly between 0 and 2 pi.
*
* @param a_X [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::IncoherentPhotoAtomicScattering::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
double k1 = a_X / PoPI_electronMass_MeV_c2;
double energyOut, mu, scatteringFactor;
if( a_X >= m_energies.back( ) ) {
MCGIDI_sampleKleinNishina( k1, a_rng, &energyOut, &mu ); }
else {
double scatteringFactorMax = evaluateScatteringFactor( a_X );
do {
MCGIDI_sampleKleinNishina( k1, a_rng, &energyOut, &mu );
scatteringFactor = evaluateScatteringFactor( a_X * sqrt( 0.5 * ( 1.0 - mu ) ) );
} while( scatteringFactor < a_rng( ) * scatteringFactorMax );
}
a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
a_input.m_energyOut1 = energyOut * PoPI_electronMass_MeV_c2;
a_input.m_mu = mu;
a_input.m_phi = 2.0 * M_PI * a_rng( );
a_input.m_frame = productFrame( );
}
/* *********************************************************************************************************//**
* This method samples the outgoing product data by sampling the outgoing energy E' from the probability P(E'|E) and then samples mu from
* the probability P(mu|E,E'). It also samples the outgoing phi uniformly between 0 and 2 pi.
*
* @param a_X [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_userrng [in] A random number generator that takes the state *a_rngState* and returns a double in the range [0.0, 1.0).
* @param a_rngState [in] The current state for the random number generator.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
double energyOut, mu, occupationNumber;
// Convert incident photon energy [MeV] to units of rest mass energy of the electron
const double alpha_in = a_X / PoPI_electronMass_MeV_c2;
double alpha_ratio, occupation_pz, occupationNumberMax;
double quad_a = 0, quad_b = 0, quad_c = 0, pz = 0; // Initialize with dummy values to silence compiler warnings
bool energetically_possible = false;
int ep_it = 0;
while( energetically_possible == false && ep_it < 1000 ){
// Sample outgoing angle
occupationNumberMax = evaluateOccupationNumber( a_X, -1.0 );
if( a_X >= 10.0 ) { // This condition is not yet correct
MCGIDI_sampleKleinNishina( alpha_in, a_rng, &energyOut, &mu ); }
else {
do {
MCGIDI_sampleKleinNishina( alpha_in, a_rng, &energyOut, &mu );
occupationNumber = evaluateOccupationNumber( a_X, mu );
} while( occupationNumber < occupationNumberMax * a_rng( ) );
}
// Sample electron momentum projection, pz
occupation_pz = occupationNumberMax*a_rng();
int intLowerIndex = binarySearchVector( occupation_pz, m_occupationNumber );
if( intLowerIndex == -1 ){
pz = m_pz.back();
}
else{
std::size_t lowerIndex = static_cast<std::size_t>( intLowerIndex );
pz = m_pz[lowerIndex] + (occupation_pz-m_occupationNumber[lowerIndex])*(m_pz[lowerIndex+1]-m_pz[lowerIndex])/(m_occupationNumber[lowerIndex+1]-m_occupationNumber[lowerIndex]);
}
// Convert pz to outgoing photon energy
alpha_ratio = energyRatio(a_X, mu);
quad_a = pz*pz - (1/alpha_ratio)*(1/alpha_ratio);
quad_b = -2*alpha_in*( pz*pz * mu - (1/alpha_ratio));
quad_c = alpha_in*alpha_in*( pz*pz - 1 );
if(quad_b*quad_b - 4*quad_a*quad_c > 0){
energetically_possible = true;
}
ep_it = ep_it + 1;
}
const double quad_1 = -quad_b/(2*quad_a) + sqrt( quad_b*quad_b - 4*quad_a*quad_c )/( 2*quad_a );
const double quad_2 = -quad_b/(2*quad_a) - sqrt( quad_b*quad_b - 4*quad_a*quad_c )/( 2*quad_a );
// Select the correct outgoing energy based on the pz value
if(pz >= 0){
if(quad_1 >= quad_2){
energyOut = quad_1;
}
else{
energyOut = quad_2;
}
}
else{
if(quad_1 >= quad_2){
energyOut = quad_2;
}
else{
energyOut = quad_1;
}
}
a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
a_input.m_energyOut1 = energyOut * PoPI_electronMass_MeV_c2;
a_input.m_mu = mu;
a_input.m_phi = 2.0 * M_PI * a_rng( );
a_input.m_frame = productFrame( );
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [in] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::IncoherentPhotoAtomicScattering::angleBiasing( Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
double a_energy_in, double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
URR_protareInfos URR_protareInfos1;
double sigma = a_reaction->protareSingle( )->reactionCrossSection( a_reaction->reactionIndex( ), URR_protareInfos1, 0.0, a_energy_in );
double norm = M_PI * MCGIDI_classicalElectronRadius * MCGIDI_classicalElectronRadius / sigma;
double one_minus_mu = 1.0 - a_mu_lab;
double k_in = a_energy_in / PoPI_electronMass_MeV_c2;
a_energy_out = a_energy_in / ( 1.0 + k_in * one_minus_mu );
double k_out = a_energy_out / PoPI_electronMass_MeV_c2;
double k_ratio = k_out / k_in;
double probability = evaluateScatteringFactor( a_energy_in * sqrt( 0.5 * one_minus_mu ) );
probability *= k_ratio * k_ratio * ( 1.0 + a_mu_lab * a_mu_lab + k_in * k_out * one_minus_mu * one_minus_mu ) * norm;
return( probability );
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_userrng [in] A random number generator that takes the state *a_rngState* and returns a double in the range [0.0, 1.0).
* @param a_rngState [in] The current state for the random number generator.
* @param a_energy_out [in] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::angleBiasing( Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
double a_energy_in, double a_mu_lab, RNG && a_rng, double &a_energy_out ) const {
URR_protareInfos URR_protareInfos1;
double sigma = a_reaction->protareSingle( )->reactionCrossSection( a_reaction->reactionIndex( ), URR_protareInfos1, 0.0, a_energy_in );
double norm = M_PI * MCGIDI_classicalElectronRadius * MCGIDI_classicalElectronRadius / sigma;
double one_minus_mu = 1.0 - a_mu_lab;
double alpha_in = a_energy_in / PoPI_electronMass_MeV_c2;
double quad_a, quad_b, quad_c, alpha_ratio, pz, occupation_pz, occupationNumberMax;
bool energetically_possible = false;
int ep_it = 0;
while( energetically_possible == false && ep_it < 1000 ){
// Sample electron momentum projection, pz
occupationNumberMax = evaluateOccupationNumber( a_energy_in, -1.0 );
occupation_pz = occupationNumberMax*a_rng();
int intLowerIndex = binarySearchVector( occupation_pz, m_occupationNumber );
pz = 0;
if( intLowerIndex == -1 ){
pz = m_pz.back();
}
else{
std::size_t lowerIndex = static_cast<std::size_t>( intLowerIndex );
pz = m_pz[lowerIndex] + (occupation_pz-m_occupationNumber[lowerIndex])*(m_pz[lowerIndex+1]-m_pz[lowerIndex])/(m_occupationNumber[lowerIndex+1]-m_occupationNumber[lowerIndex]);
}
// Convert pz to outgoing photon energy
alpha_ratio = energyRatio(a_energy_in, a_mu_lab);
quad_a = pz*pz - (1/alpha_ratio)*(1/alpha_ratio);
quad_b = -2*alpha_in*( pz*pz * a_mu_lab - (1/alpha_ratio));
quad_c = alpha_in*alpha_in*( pz*pz - 1 );
if(quad_b*quad_b - 4*quad_a*quad_c > 0){
energetically_possible = true;
}
ep_it = ep_it + 1;
}
const double quad_1 = -quad_b/(2*quad_a) + sqrt( quad_b*quad_b - 4*quad_a*quad_c )/( 2*quad_a );
const double quad_2 = -quad_b/(2*quad_a) - sqrt( quad_b*quad_b - 4*quad_a*quad_c )/( 2*quad_a );
// Select the correct outgoing energy based on the pz value
double alpha_out = 0;
if(pz >= 0){
if(quad_1 >= quad_2){
alpha_out = quad_1;
}
else{
alpha_out = quad_2;
}
}
else{
if(quad_1 >= quad_2){
alpha_out = quad_2;
}
else{
alpha_out = quad_1;
}
}
alpha_ratio = alpha_out / alpha_in;
a_energy_out = alpha_out * PoPI_electronMass_MeV_c2;
double probability = evaluateOccupationNumber( a_energy_in, a_mu_lab );
probability *= alpha_ratio * alpha_ratio * ( 1.0 + a_mu_lab * a_mu_lab + alpha_in * alpha_out * one_minus_mu * one_minus_mu ) * norm;
return( probability );
}
/* *********************************************************************************************************//**
* This method returns the outgoing electron energy and angle given that the photon when out at an angle of *a_input.m_mu*.
* Ergo, this method must be called directly after the photon has been sampled.
*
* @param a_energy [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::IncoherentPhotoAtomicScatteringElectron::sample( double a_energy, Sampling::Input &a_input, RNG && a_rng ) const {
double halfTheta = 0.5 * acos( a_input.m_mu );
double cot_psi = ( 1.0 + a_energy / PoPI_electronMass_MeV_c2 ) * tan( halfTheta );
double psi = atan( 1.0 / cot_psi );
double deltaE_photon = a_energy - a_input.m_energyOut1;
double electronMomentum2 = deltaE_photon * ( deltaE_photon + 2.0 * PoPI_electronMass_MeV_c2 ); // Square of the electron outlgoing momentum.
a_input.m_energyOut1 = electronMomentum2 / ( sqrt( electronMomentum2 + PoPI_electronMass_MeV_c2 * PoPI_electronMass_MeV_c2 ) + PoPI_electronMass_MeV_c2 );
a_input.m_mu = cos( psi );
a_input.m_phi = 2.0 * M_PI * a_rng( );
a_input.m_frame = productFrame( );
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* causing a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
* Currently, this method only returns 0.0 for the probability and outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [in] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::IncoherentPhotoAtomicScatteringElectron::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
LUPI_maybeUnused double a_energy_in, LUPI_maybeUnused double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
a_energy_out = 0;
return( 0.0 );
}
/* *********************************************************************************************************//**
* This method samples the outgoing photon by assigning the electron rest mass energy as the photon's energy and,
* if m_firstSampled is true, randomly picking mu and phi. If m_firstSampled is false, the previous sampled particle
* that filled in a_input must be the other sampled photon, then, the mu and phi for the second-sampled photon is such that
* it is back-to-back with the other photon.
*
* @param a_X [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::PairProductionGamma::sample( LUPI_maybeUnused double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
if( m_firstSampled ) {
a_input.m_mu = 1.0 - 2.0 * a_rng( );
a_input.m_phi = M_PI * a_rng( ); }
else {
a_input.m_mu *= -1.0;
a_input.m_phi += M_PI;
}
a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
a_input.m_energyOut1 = PoPI_electronMass_MeV_c2;
a_input.m_frame = productFrame( );
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [in] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::PairProductionGamma::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
LUPI_maybeUnused double a_energy_in, LUPI_maybeUnused double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
a_energy_out = PoPI_electronMass_MeV_c2;
return( 1.0 ); // 1.0 as there are two photons, each with 1/2 probability.
}
/* *********************************************************************************************************//**
* This method samples the outgoing neutron data for coherent elastic TSNL from the Debye/Waller function.
*
* @param a_energy [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::CoherentElasticTNSL::sample( double a_energy, Sampling::Input &a_input,
RNG && a_rng ) const {
if( a_energy <= m_energies[0] ) {
a_input.m_mu = 1.0; }
else {
double temperature = a_input.temperature( );
if( temperature < m_temperatures[0] ) temperature = m_temperatures[0];
if( temperature > m_temperatures.back( ) ) temperature = m_temperatures.back( );
std::size_t temperatureIndex = (std::size_t) MCGIDI::binarySearchVector( temperature, m_temperatures, true );
double const *pointer1 = &m_S_table[temperatureIndex * m_energies.size( )];
double const *pointer2 = pointer1;
double fractionFirstTemperature = 1.0;
if( temperatureIndex != ( m_temperatures.size( ) - 1 ) ) {
fractionFirstTemperature = ( m_temperatures[temperatureIndex+1] - temperature ) / ( m_temperatures[temperatureIndex+1] - m_temperatures[temperatureIndex] );
pointer2 += m_energies.size( );
}
double fractionSecondTemperature = 1.0 - fractionFirstTemperature;
int intEnergyIndexMax = MCGIDI::binarySearchVector( a_energy, m_energies, true );
std::size_t energyIndexMax = static_cast<std::size_t>( intEnergyIndexMax );
if( a_energy == m_energies[energyIndexMax] ) --energyIndexMax;
double randomTotal = a_rng( ) * ( fractionFirstTemperature * pointer1[energyIndexMax] + fractionSecondTemperature * pointer2[energyIndexMax] );
std::size_t energyIndex = 0;
for( ; energyIndex < energyIndexMax; ++energyIndex ) {
if( randomTotal <= fractionFirstTemperature * pointer1[energyIndex] + fractionSecondTemperature * pointer2[energyIndex] ) break;
}
a_input.m_mu = 1.0 - 2.0 * m_energies[energyIndex] / a_energy;
}
a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
a_input.m_energyOut1 = a_energy;
a_input.m_phi = 2.0 * M_PI * a_rng( );
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [in] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::CoherentElasticTNSL::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
double a_energy_in, LUPI_maybeUnused double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
double probability = 0.0;
a_energy_out = a_energy_in;
return( probability );
}
/* *********************************************************************************************************//**
* This method samples the outgoing neutron data for incoherent elastic TSNL from the Debye/Waller function.
*
* @param a_energy [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::IncoherentElasticTNSL::sample( double a_energy, Sampling::Input &a_input,
RNG && a_rng ) const {
double temperature = a_input.temperature( ) / m_temperatureToMeV_K;
double W_prime = m_DebyeWallerIntegral->evaluate( temperature );
double twoEW = 2 * a_energy * W_prime;
double expOfTwice_twoEW = exp( -2 * twoEW );
double sampled_cdf = a_rng( );
if( sampled_cdf > ( 1 - 1e-5 ) ) {
double Variable = ( 1.0 - sampled_cdf ) * ( 1.0 - expOfTwice_twoEW );
a_input.m_mu = 1.0 - Variable * ( 1.0 + 0.5 * Variable ) / twoEW; }
else if( sampled_cdf < expOfTwice_twoEW ) {
a_input.m_mu = -1.0 + log( sampled_cdf / expOfTwice_twoEW * ( 1.0 - expOfTwice_twoEW ) + 1.0 ) / twoEW; }
else {
a_input.m_mu = 1.0 + log( expOfTwice_twoEW + sampled_cdf * ( 1.0 - expOfTwice_twoEW ) ) / twoEW;
}
a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
a_input.m_energyOut1 = a_energy;
a_input.m_phi = 2.0 * M_PI * a_rng( );
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] The temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [in] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::IncoherentElasticTNSL::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, double a_temperature,
double a_energy_in, double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
double temperature = a_temperature / m_temperatureToMeV_K;
double W_prime = m_DebyeWallerIntegral->evaluate( temperature );
double twoEW = 2 * a_energy_in * W_prime;
double probability = exp( -twoEW * ( 1.0 - a_mu_lab ) ) * twoEW / ( 1.0 - exp( -2 * twoEW ) );
a_energy_out = a_energy_in;
return( probability );
}
/* *********************************************************************************************************//**
* The method sets all outgoing product data to 0.0 and set the sampledType to Sampling::unspecified.
*
* @param a_X [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::Unspecified::sample( LUPI_maybeUnused double a_X, Sampling::Input &a_input, LUPI_maybeUnused RNG && a_rng ) const {
a_input.setSampledType( Sampling::SampledType::unspecified );
a_input.m_energyOut1 = 0.;
a_input.m_mu = 0.;
a_input.m_phi = 0.;
a_input.m_frame = productFrame( );
}
/* *********************************************************************************************************//**
* Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
* at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy. This one should never
* be called. If called, returns 0.0 for a probability.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_temperature [in] Specifies the temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_energy_out [in] The energy of the emitted outgoing particle.
*
* @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::Unspecified::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
LUPI_maybeUnused double a_energy_in, LUPI_maybeUnused double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
a_energy_out = 0.0;
return( 0.0 );
}
// From file: MCGIDI_functions.cpp
/*
============================================================
*/
template <typename RNG >
LUPI_HOST_DEVICE int MCGIDI::Functions::Function1d::sampleBoundingInteger( double a_x1, RNG && a_rng ) const {
if( type( ) == Function1dType::TerrellFissionNeutronMultiplicityModel )
return( static_cast<TerrellFissionNeutronMultiplicityModel const *>( this )->sampleBoundingInteger( a_x1, a_rng ) );
double d_value = evaluate( a_x1 );
int iValue = (int) d_value;
if( iValue == d_value ) return( iValue );
if( d_value - iValue > a_rng( ) ) ++iValue;
return( iValue );
}
/* *********************************************************************************************************//**
* Sample the number of fission prompt neutrons using Terrell's modified Gaussian distribution.
* Method uses Red Cullen's algoritm (see UCRL-TR-222526).
*
* @param a_energy [in] The energy of the projectile.
* @param a_rng [in] The random number generator function the uses *a_rngState* to generator a double in the range [0, 1.0).
* @param a_rngState [in/out] The random number generator state.
*
* @return The sampled number of emitted, prompt neutrons for fission.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE int MCGIDI::Functions::TerrellFissionNeutronMultiplicityModel::sampleBoundingInteger( double a_energy, RNG && a_rng ) const {
const double Terrell_BSHIFT = -0.43287;
double width = M_SQRT2 * m_width;
double temp1 = m_multiplicity->evaluate( a_energy ) + 0.5;
double temp2 = temp1 / width;
double expo = exp( -temp2 * temp2 );
double cshift = temp1 + Terrell_BSHIFT * m_width * expo / ( 1.0 - expo );
double multiplicity = 1.0;
do {
double rw = sqrt( -log( a_rng( ) ) );
double theta = ( 2.0 * M_PI ) * a_rng();
multiplicity = width * rw * cos( theta ) + cshift;
} while ( multiplicity < 0.0 );
return( static_cast<int>( floor( multiplicity ) ) );
}
/* *********************************************************************************************************//**
* Returns the x-value corresponding cumulative probability *a_rngValue*.
*
* @param a_rngValue [in] The x-value to evaluate the function at.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
*
* @return The value of the function at *a_x1*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::ProbabilityBase1d::sample( double a_rngValue, RNG && a_rng ) const {
return( static_cast<Xs_pdf_cdf1d const *>( this )->sample( a_rngValue, a_rng ) );
}
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::Xs_pdf_cdf1d::sample( double a_rngValue, LUPI_maybeUnused RNG && a_rng ) const {
int intLower = binarySearchVector( a_rngValue, m_cdf );
double domainValue = 0;
if( intLower < 0 ) { // This should never happen.
LUPI_THROW( "Xs_pdf_cdf1d::sample: intLower < 0." );
}
std::size_t lower = static_cast<std::size_t>( intLower );
if( interpolation( ) == Interpolation::FLAT ) {
double fraction = ( m_cdf[lower+1] - a_rngValue ) / ( m_cdf[lower+1] - m_cdf[lower] );
domainValue = fraction * m_Xs[lower] + ( 1 - fraction ) * m_Xs[lower+1]; }
else { // Assumes lin-lin interpolation.
double slope = m_pdf[lower+1] - m_pdf[lower];
if( slope == 0.0 ) {
if( m_pdf[lower] == 0.0 ) {
domainValue = m_Xs[lower];
if( lower == 0 ) domainValue = m_Xs[1]; }
else {
double fraction = ( m_cdf[lower+1] - a_rngValue ) / ( m_cdf[lower+1] - m_cdf[lower] );
domainValue = fraction * m_Xs[lower] + ( 1 - fraction ) * m_Xs[lower+1];
} }
else {
double d1, d2;
slope = slope / ( m_Xs[lower+1] - m_Xs[lower] );
d1 = a_rngValue - m_cdf[lower];
d2 = m_cdf[lower+1] - a_rngValue;
if( d2 > d1 ) { // Closer to lower.
domainValue = m_Xs[lower] + ( sqrt( m_pdf[lower] * m_pdf[lower] + 2. * slope * d1 ) - m_pdf[lower] ) / slope; }
else { // Closer to lower + 1.
domainValue = m_Xs[lower+1] - ( m_pdf[lower+1] - sqrt( m_pdf[lower+1] * m_pdf[lower+1] - 2. * slope * d2 ) ) / slope;
}
}
}
return( domainValue );
}
/* *********************************************************************************************************//**
* This method samples an x1 from a pdf(x1|x2) given x2 and the cumulative value of the pdf as *a_rngValue*.
*
* @param a_x2 [in] The value of x2.
* @param a_rngValue [in] The value of the cumulative used to determine the x1 value.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::ProbabilityBase2d::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
double value = 0.0;
switch( type( ) ) {
case ProbabilityBase2dType::none:
break;
case ProbabilityBase2dType::weightedFunctionals:
value = static_cast<WeightedFunctionals2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
default:
value = static_cast<ProbabilityBase2d_d1 const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
}
return( value );
}
/* *********************************************************************************************************//**
* Returns the value of x1, given x2 and the cumulative probability *a_rngValue*.
*
* @param a_x2 [in] Value of the outer most independent variable (i.e., *x2*).
* @param a_rngValue [in] The value of the cumulative probability used to determine the x1 value.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
*
* @return The value *x1* where the cumulative probability is *a_rngValue* for x2 = *a_x2*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::ProbabilityBase2d_d1::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
double value = 0.0;
switch( type( ) ) {
case ProbabilityBase2dType::XYs:
case ProbabilityBase2dType::isotropic:
case ProbabilityBase2dType::discreteGamma:
case ProbabilityBase2dType::primaryGamma:
case ProbabilityBase2dType::recoil:
case ProbabilityBase2dType::NBodyPhaseSpace:
case ProbabilityBase2dType::evaporation:
case ProbabilityBase2dType::generalEvaporation:
case ProbabilityBase2dType::simpleMaxwellianFission:
case ProbabilityBase2dType::Watt:
value = static_cast<ProbabilityBase2d_d2 const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
case ProbabilityBase2dType::regions:
value = static_cast<Regions2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
case ProbabilityBase2dType::none:
case ProbabilityBase2dType::weightedFunctionals:
LUPI_THROW( "ProbabilityBase2d_d1::sample: This should never happen." );
}
return( value );
}
/* *********************************************************************************************************//**
* This method returns two x1 values for use with ProbabilityBase3d functions.
*
* @param a_x2 [in] The value of x2.
* @param a_rngValue [in] The value of the cumulative value used to determine the x1 value.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_x1_1 [in] The lower value of the x1 value.
* @param a_x1_2 [in] The upper value of the x1 value.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::ProbabilityBase2d_d1::sample2dOf3d( double a_x2, double a_rngValue, RNG && a_rng,
double *a_x1_1, double *a_x1_2 ) const {
double value = 0.0;
switch( type( ) ) {
case ProbabilityBase2dType::XYs:
value = static_cast<XYs2d const *>( this )->sample2dOf3d( a_x2, a_rngValue, a_rng, a_x1_1, a_x1_2 );
break;
default:
LUPI_THROW( "ProbabilityBase2d_d1::sample2dOf3d: not implemented." );
}
return( value );
}
/* *********************************************************************************************************//**
* Returns the value of x1, given x2 and the cumulative probability *a_rngValue*.
*
* @param a_x2 [in] Value of the outer most independent variable (i.e., *x2*).
* @param a_rngValue [in] The value of the cumulative probability used to determine the x1 value.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
*
* @return The value *x1* where the cumulative probability is *a_rngValue* for x2 = *a_x2*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::ProbabilityBase2d_d2::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
double value = 0.0;
switch( type( ) ) {
case ProbabilityBase2dType::XYs:
value = static_cast<XYs2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
case ProbabilityBase2dType::isotropic:
value = static_cast<Isotropic2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
case ProbabilityBase2dType::discreteGamma:
value = static_cast<DiscreteGamma2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
case ProbabilityBase2dType::primaryGamma:
value = static_cast<PrimaryGamma2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
case ProbabilityBase2dType::recoil:
value = static_cast<Recoil2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
case ProbabilityBase2dType::NBodyPhaseSpace:
value = static_cast<NBodyPhaseSpace2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
case ProbabilityBase2dType::evaporation:
value = static_cast<Evaporation2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
case ProbabilityBase2dType::generalEvaporation:
value = static_cast<GeneralEvaporation2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
case ProbabilityBase2dType::simpleMaxwellianFission:
value = static_cast<SimpleMaxwellianFission2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
case ProbabilityBase2dType::Watt:
value = static_cast<Watt2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
break;
case ProbabilityBase2dType::none:
case ProbabilityBase2dType::weightedFunctionals:
case ProbabilityBase2dType::regions:
LUPI_THROW( "ProbabilityBase2d_d2::sample: This should never happen." );
}
return( value );
}
/* *********************************************************************************************************//**
* This method returns two x1 values for use with ProbabilityBase3d functions.
*
* @param a_x2 [in] The value of x2.
* @param a_rngValue [in] The value of the cumulative value used to determine the x1 value.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_x1_1 [in] The lower value of the x1 value.
* @param a_x1_2 [in] The upper value of the x1 value.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::ProbabilityBase2d_d2::sample2dOf3d( double a_x2, double a_rngValue, RNG && a_rng,
double *a_x1_1, double *a_x1_2 ) const {
double value = 0.0;
switch( type( ) ) {
case ProbabilityBase2dType::XYs:
value = static_cast<XYs2d const *>( this )->sample2dOf3d( a_x2, a_rngValue, a_rng, a_x1_1, a_x1_2 );
break;
default:
LUPI_THROW( "ProbabilityBase2d_d2::sample2dOf3d: not implemented." );
}
return( value );
}
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::XYs2d::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
/*
C Samples from a pdf(x1|x2). First determine which pdf(s) to sample from given x2.
C Then use rngValue to sample from pdf1(x1) and maybe pdf2(x1) and interpolate to
C determine x1.
*/
double sampledValue = 0;
int intLower = binarySearchVector( a_x2, m_Xs );
if( intLower == -2 ) {
sampledValue = m_probabilities[0]->sample( a_rngValue, a_rng ); }
else if( intLower == -1 ) {
sampledValue = m_probabilities.back( )->sample( a_rngValue, a_rng ); }
else {
auto lower = static_cast<std::size_t>( intLower );
double sampled1 = m_probabilities[lower]->sample( a_rngValue, a_rng );
if( interpolation( ) == Interpolation::FLAT ) {
sampledValue = sampled1; }
else {
double sampled2 = m_probabilities[lower+1]->sample( a_rngValue, a_rng );
if( interpolation( ) == Interpolation::LINLIN ) {
double fraction = ( m_Xs[lower+1] - a_x2 ) / ( m_Xs[lower+1] - m_Xs[lower] );
sampledValue = fraction * sampled1 + ( 1 - fraction ) * sampled2; }
else if( interpolation( ) == Interpolation::LOGLIN ) {
double fraction = ( m_Xs[lower+1] - a_x2 ) / ( m_Xs[lower+1] - m_Xs[lower] );
sampledValue = sampled2 * pow( sampled2 / sampled1, fraction ); }
else if( interpolation( ) == Interpolation::LINLOG ) {
double fraction = log( m_Xs[lower+1] / a_x2 ) / log( m_Xs[lower+1] / m_Xs[lower] );
sampledValue = fraction * sampled1 + ( 1 - fraction ) * sampled2; }
else if( interpolation( ) == Interpolation::LOGLOG ) {
double fraction = log( m_Xs[lower+1] / a_x2 ) / log( m_Xs[lower+1] / m_Xs[lower] );
sampledValue = sampled2 * pow( sampled2 / sampled1, fraction ); }
else { // This should never happen.
LUPI_THROW( "XYs2d::sample: unsupported interpolation." );
}
}
}
return( sampledValue );
}
/* *********************************************************************************************************//**
* This method returns two x1 values for use with ProbabilityBase3d functions.
*
* @param a_x2 [in] The value of x2.
* @param a_rngValue [in] The value of the cumulative value used to determine the x1 value.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_x1_1 [in] The lower value of the x1 value.
* @param a_x1_2 [in] The upper value of the x1 value.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::XYs2d::sample2dOf3d( double a_x2, double a_rngValue, RNG && a_rng,
double *a_x1_1, double *a_x1_2 ) const {
/*
C Samples from a pdf(x1|x2). First determine which pdf(s) to sample from given x2. Then use rngValue to sample from pdf1(x1)
C and maybe pdf2(x1) and interpolate to determine x1.
*/
double sampledValue = 0;
int intLower = binarySearchVector( a_x2, m_Xs );
if( intLower == -2 ) {
sampledValue = m_probabilities[0]->sample( a_rngValue, a_rng );
*a_x1_2 = *a_x1_1 = sampledValue; }
else if( intLower == -1 ) {
sampledValue = m_probabilities.back( )->sample( a_rngValue, a_rng );
*a_x1_2 = *a_x1_1 = sampledValue; }
else {
std::size_t lower = static_cast<std::size_t>( intLower );
*a_x1_1 = m_probabilities[lower]->sample( a_rngValue, a_rng );
if( interpolation( ) == Interpolation::FLAT ) {
sampledValue = *a_x1_2 = *a_x1_1; }
else {
*a_x1_2 = m_probabilities[lower+1]->sample( a_rngValue, a_rng );
if( interpolation( ) == Interpolation::LINLIN ) {
double fraction = ( m_Xs[lower+1] - a_x2 ) / ( m_Xs[lower+1] - m_Xs[lower] );
sampledValue = fraction * *a_x1_1 + ( 1 - fraction ) * *a_x1_2; }
else if( interpolation( ) == Interpolation::LOGLIN ) {
double fraction = ( m_Xs[lower+1] - a_x2 ) / ( m_Xs[lower+1] - m_Xs[lower] );
sampledValue = *a_x1_2 * pow( *a_x1_2 / *a_x1_1, fraction ); }
else if( interpolation( ) == Interpolation::LINLOG ) {
double fraction = log( m_Xs[lower+1] / a_x2 ) / log( m_Xs[lower+1] / m_Xs[lower] );
sampledValue = fraction * *a_x1_1 + ( 1 - fraction ) * *a_x1_2; }
else if( interpolation( ) == Interpolation::LOGLOG ) {
double fraction = log( m_Xs[lower+1] / a_x2 ) / log( m_Xs[lower+1] / m_Xs[lower] );
sampledValue = *a_x1_2 * pow( *a_x1_2 / *a_x1_1 , fraction ); }
else { // This should never happen.
LUPI_THROW( "XYs2d::sample: unsupported interpolation." );
}
}
}
return( sampledValue );
}
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::Regions2d::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
int intLower = binarySearchVector( a_x2, m_Xs );
if( intLower < 0 ) {
if( intLower == -1 ) { // a_x2 > last value of m_Xs.
return( m_probabilities.back( )->sample( a_x2, a_rngValue, a_rng ) );
}
intLower = 0; // a_x2 < first value of m_Xs.
}
std::size_t lower = static_cast<std::size_t>( intLower );
return( m_probabilities[lower]->sample( a_x2, a_rngValue, a_rng ) );
}
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::Recoil2d::sample( LUPI_maybeUnused double a_x2, LUPI_maybeUnused double a_rngValue, LUPI_maybeUnused RNG && a_rng ) const {
#if !defined(__NVCC__) && !defined(__HIP__)
LUPI_THROW( "Recoil2d::sample: not implemented." );
#endif
return( 0.0 );
}
/* *********************************************************************************************************
* Sampling for NBody phase space.
*
* @param a_x2 [in] Incident energy of the projectile.
* @param a_rngValue [in] The GIDI::Protare whose data is to be used to construct *this*.
* @param a_rng [in] This argument is not used by this method but needed to match ProbabilityBase1d::sample's definition.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::NBodyPhaseSpace2d::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
double energyMax = m_energy_in_COMFactor * a_x2 + m_Q;
if( energyMax < 0.0 ) return( 0.0 ); // Kludge for now until for upscatter model A until sampling below threshold is fixed.
return( energyMax * m_massFactor * m_dist->sample( a_rngValue, a_rng ) );
}
inline LUPI_HOST_DEVICE static double MCGIDI_sampleEvaporation( double a_xMax, double a_rngValue ) {
double b1, c1, xMid, norm, xMin = 0.;
norm = 1 - ( 1 + a_xMax ) * exp( -a_xMax );
b1 = 1. - norm * a_rngValue;
for( int i1 = 0; i1 < 16; i1++ ) {
xMid = 0.5 * ( xMin + a_xMax );
c1 = ( 1 + xMid ) * exp( -xMid );
if( b1 > c1 ) {
a_xMax = xMid; }
else {
xMin = xMid;
}
}
return( xMid );
}
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::Evaporation2d::sample( double a_x2, double a_rngValue, LUPI_maybeUnused RNG && a_rng ) const {
double theta = m_theta->evaluate( a_x2 );
return( theta * MCGIDI_sampleEvaporation( ( a_x2 - m_U ) / theta, a_rngValue ) );
}
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::GeneralEvaporation2d::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
return( m_theta->evaluate( a_x2 ) * m_g->sample( a_rngValue, a_rng ) );
}
inline LUPI_HOST_DEVICE static double MCGIDI_sampleSimpleMaxwellianFission( double a_xMax, double a_rngValue ) {
double b1, c1, xMid, norm, xMin = 0., sqrt_xMid, sqrt_pi_2 = 0.5 * sqrt( M_PI );
sqrt_xMid = sqrt( a_xMax );
norm = sqrt_pi_2 * erf( sqrt_xMid ) - sqrt_xMid * exp( -a_xMax );
b1 = norm * a_rngValue;
for( int i1 = 0; i1 < 16; i1++ ) {
xMid = 0.5 * ( xMin + a_xMax );
sqrt_xMid = sqrt( xMid );
c1 = sqrt_pi_2 * erf( sqrt_xMid ) - sqrt_xMid * exp( -xMid );
if( b1 < c1 ) {
a_xMax = xMid; }
else {
xMin = xMid;
}
}
return( xMid );
}
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::SimpleMaxwellianFission2d::sample( double a_x2, double a_rngValue, LUPI_maybeUnused RNG && a_rng ) const {
double theta = m_theta->evaluate( a_x2 );
return( theta * MCGIDI_sampleSimpleMaxwellianFission( ( a_x2 - m_U ) / theta, a_rngValue ) );
}
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::Watt2d::sample( double a_x2, LUPI_maybeUnused double a_rngValue, RNG && a_rng ) const {
/*
* From MCAPM via Sample Watt Spectrum as in TART ( Kalos algorithm ).
*/
double WattMin = 0., WattMax = a_x2 - m_U, x, y, z, energyOut, rand1, rand2;
double Watt_a = 1./m_a->evaluate( a_x2 ); // Kalos algorithm uses the inverse of the 'a' parameter stored in GNDS
double Watt_b = m_b->evaluate( a_x2 );
x = 1. + ( Watt_b / ( 8. * Watt_a ) );
y = ( x + sqrt( x * x - 1. ) ) / Watt_a;
z = Watt_a * y - 1.;
do {
rand1 = -log( a_rng( ) );
rand2 = -log( a_rng( ) );
energyOut = y * rand1;
} while( ( ( rand2 - z * ( rand1 + 1. ) ) * ( rand2 - z * ( rand1 + 1. ) ) > Watt_b * y * rand1 ) ||
( energyOut < WattMin ) || ( energyOut > WattMax ) );
return( energyOut );
}
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::WeightedFunctionals2d::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
/*
c This routine assumes that the weights sum to 1.
*/
std::size_t i1;
std::size_t n1 = m_weight.size( ) - 1; // Take last point if others do not add to randomWeight.
double randomWeight = a_rng( ), cumulativeWeight = 0.;
for( i1 = 0; i1 < n1; ++i1 ) {
cumulativeWeight += m_weight[i1]->evaluate( a_x2 );
if( cumulativeWeight >= randomWeight ) break;
}
return( m_energy[i1]->sample( a_x2, a_rngValue, a_rng) );
}
/* *********************************************************************************************************//**
* This method samples an x1 from a pdf(x1|x2) given x2 and the cumulative value of the pdf as *a_rngValue*.
*
* @param a_x3 [in] The value of x3.
* @param a_x2_1 [in] The value of ?.
* @param a_x2_2 [in] The value of ?.
* @param a_rngValue [in] The value of the cumulative used to determine the x1 value.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::ProbabilityBase3d::sample( double a_x3, double a_x2_1, double a_x2_2, double a_rngValue, RNG && a_rng ) const {
return( static_cast<XYs3d const *>( this )->sample( a_x3, a_x2_1, a_x2_2, a_rngValue, a_rng ) );
}
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Probabilities::XYs3d::sample( double a_x3, double a_x2_1, double a_x2_2, double a_rngValue, RNG && a_rng ) const {
/*
C Samples from a pdf(x1|x3,x2). First determine which pdf(s) to sample from given x3
C Then use rngValue to sample from pdf2_1(x2) and maybe pdf2_2(x2) and interpolate to
C determine x1.
*/
double sampledValue = 0;
int intLower = binarySearchVector( a_x3, m_Xs );
if( intLower == -2 ) { // x3 < first value of Xs.
sampledValue = m_probabilities[0]->sample( a_x2_1, a_rngValue, a_rng ); }
else if( intLower == -1 ) { // x3 > last value of Xs.
sampledValue = m_probabilities.back( )->sample( a_x2_1, a_rngValue, a_rng ); }
else {
std::size_t lower = static_cast<std::size_t>( intLower );
double sampled1 = m_probabilities[lower]->sample( a_x2_1, a_rngValue, a_rng );
if( interpolation( ) == Interpolation::FLAT ) {
sampledValue = sampled1; }
else {
double sampled2 = m_probabilities[lower+1]->sample( a_x2_2, a_rngValue, a_rng );
if( interpolation( ) == Interpolation::LINLIN ) {
double fraction = ( m_Xs[lower+1] - a_x3 ) / ( m_Xs[lower+1] - m_Xs[lower] );
sampledValue = fraction * sampled1 + ( 1 - fraction ) * sampled2; }
else if( interpolation( ) == Interpolation::LOGLIN ) {
double fraction = ( m_Xs[lower+1] - a_x3 ) / ( m_Xs[lower+1] - m_Xs[lower] );
sampledValue = sampled2 * pow( sampled2 / sampled1, fraction ); }
else if( interpolation( ) == Interpolation::LINLOG ) {
double fraction = log( m_Xs[lower+1] / a_x3 ) / log( m_Xs[lower+1] / m_Xs[lower] );
sampledValue = fraction * sampled1 + ( 1 - fraction ) * sampled2; }
else if( interpolation( ) == Interpolation::LOGLOG ) {
double fraction = log( m_Xs[lower+1] / a_x3 ) / log( m_Xs[lower+1] / m_Xs[lower] );
sampledValue = sampled2 * pow( sampled2 / sampled1, fraction ); }
else { // This should never happen.
LUPI_THROW( "XYs3d::sample: unsupported interpolation." );
}
}
}
return( sampledValue );
}
// From file: MCGIDI_heatedCrossSections.cpp
/* *********************************************************************************************************//**
* Returns the requested reaction's multi-group cross section for target temperature *a_temperature* and projectile multi-group *a_hashIndex*.
*
* @param a_URR_protareInfos [in] URR information.
* @param a_URR_index [in] If not negative, specifies the index in *a_URR_protareInfos*.
* @param a_hashIndex [in] Specifies projectile energy hash index.
* @param a_temperature [in] The temperature of the target.
* @param a_energy [in] The energy of the projectile.
* @param a_crossSection [in] The total cross section for the protare at *a_temperature* and *a_energy*.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE std::size_t MCGIDI::HeatedCrossSectionsContinuousEnergy::sampleReaction( URR_protareInfos const &a_URR_protareInfos,
int a_URR_index, std::size_t a_hashIndex, double a_temperature, double a_energy, double a_crossSection, RNG && a_rng ) const {
std::size_t sampled_reaction_index, temperatureIndex1, temperatureIndex2, number_of_temperatures = m_temperatures.size( );
double sampleCrossSection = a_crossSection * a_rng( );
if( a_temperature <= m_temperatures[0] ) {
temperatureIndex1 = 0;
temperatureIndex2 = temperatureIndex1; }
else if( a_temperature >= m_temperatures.back( ) ) {
temperatureIndex1 = m_temperatures.size( ) - 1;
temperatureIndex2 = temperatureIndex1; }
else {
std::size_t i1 = 0;
for( ; i1 < number_of_temperatures; ++i1 ) if( a_temperature < m_temperatures[i1] ) break;
temperatureIndex1 = i1 - 1;
temperatureIndex2 = i1;
}
std::size_t numberOfReactions = m_heatedCrossSections[0]->numberOfReactions( );
double energyFraction1, energyFraction2, crossSectionSum = 0.0;
HeatedCrossSectionContinuousEnergy &heatedCrossSection1 = *m_heatedCrossSections[temperatureIndex1];
std::size_t energyIndex1 = heatedCrossSection1.evaluationInfo( a_hashIndex, a_energy, &energyFraction1 );
if( temperatureIndex1 == temperatureIndex2 ) {
for( sampled_reaction_index = 0; sampled_reaction_index < numberOfReactions; ++sampled_reaction_index ) {
crossSectionSum += heatedCrossSection1.reactionCrossSection2( sampled_reaction_index, a_URR_protareInfos, a_URR_index, a_energy,
energyIndex1, energyFraction1 );
if( crossSectionSum >= sampleCrossSection ) break;
} }
else {
double temperatureFraction2 = ( a_temperature - m_temperatures[temperatureIndex1] )
/ ( m_temperatures[temperatureIndex2] - m_temperatures[temperatureIndex1] );
double temperatureFraction1 = 1.0 - temperatureFraction2;
HeatedCrossSectionContinuousEnergy &heatedCrossSection2 = *m_heatedCrossSections[temperatureIndex2];
std::size_t energyIndex2 = heatedCrossSection2.evaluationInfo( a_hashIndex, a_energy, &energyFraction2 );
for( sampled_reaction_index = 0; sampled_reaction_index < numberOfReactions; ++sampled_reaction_index ) {
if( m_thresholds[sampled_reaction_index] >= a_energy ) continue;
crossSectionSum += temperatureFraction1 * heatedCrossSection1.reactionCrossSection2( sampled_reaction_index, a_URR_protareInfos,
a_URR_index, a_energy, energyIndex1, energyFraction1 );
crossSectionSum += temperatureFraction2 * heatedCrossSection2.reactionCrossSection2( sampled_reaction_index, a_URR_protareInfos,
a_URR_index, a_energy, energyIndex2, energyFraction2 );
if( crossSectionSum >= sampleCrossSection ) break;
}
}
if( sampled_reaction_index == numberOfReactions ) {
if( crossSectionSum < ( 1.0 - crossSectionSumError ) * a_crossSection ) {
#if LUPI_ON_GPU
MCGIDI_PRINTF( "HeatedCrossSectionsContinuousEnergy::sampleReaction: crossSectionSum %.17e less than a_crossSection = %.17e.",
crossSectionSum, a_crossSection );
#else
std::string errorString = "HeatedCrossSectionsContinuousEnergy::sampleReaction: crossSectionSum "
+ LUPI::Misc::doubleToString3( "%.17e", crossSectionSum ) + " less than a_crossSection = "
+ LUPI::Misc::doubleToString3( "%.17e", a_crossSection ) + ".";
LUPI_THROW( errorString.c_str( ) );
#endif
}
for( sampled_reaction_index = 0; sampled_reaction_index < numberOfReactions; ++sampled_reaction_index ) { // This should rarely happen so just pick the first reaction with non-zero cross section.
if( heatedCrossSection1.reactionCrossSection2( sampled_reaction_index, a_URR_protareInfos, a_URR_index, a_energy, energyIndex1,
energyFraction1, true ) > 0 ) break;
}
}
return( sampled_reaction_index );
}
/* *********************************************************************************************************//**
* Returns the requested reaction's multi-group cross section for target temperature *a_temperature* and projectile multi-group *a_hashIndex*.
*
* @param a_hashIndex [in] The multi-group index.
* @param a_temperature [in] The temperature of the target.
* @param a_energy [in] The energy of the projectile.
* @param a_crossSection [in] The total cross section for the protare at *a_temperature* and *a_energy*.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE std::size_t MCGIDI::HeatedCrossSectionsMultiGroup::sampleReaction( std::size_t a_hashIndex, double a_temperature, double a_energy, double a_crossSection,
RNG && a_rng ) const {
std::size_t i1, sampled_reaction_index, temperatureIndex1, temperatureIndex2, numberOfTemperatures = m_temperatures.size( );
double sampleCrossSection = a_crossSection * a_rng( );
if( a_temperature <= m_temperatures[0] ) {
temperatureIndex1 = 0;
temperatureIndex2 = temperatureIndex1; }
else if( a_temperature >= m_temperatures.back( ) ) {
temperatureIndex1 = m_temperatures.size( ) - 1;
temperatureIndex2 = temperatureIndex1; }
else {
for( i1 = 0; i1 < numberOfTemperatures; ++i1 ) if( a_temperature < m_temperatures[i1] ) break;
temperatureIndex1 = i1 - 1;
temperatureIndex2 = i1;
}
std::size_t numberOfReactions = m_heatedCrossSections[0]->numberOfReactions( );
double crossSectionSum = 0;
HeatedCrossSectionMultiGroup &heatedCrossSection1 = *m_heatedCrossSections[temperatureIndex1];
if( temperatureIndex1 == temperatureIndex2 ) {
for( sampled_reaction_index = 0; sampled_reaction_index < numberOfReactions; ++sampled_reaction_index ) {
crossSectionSum += heatedCrossSection1.reactionCrossSection( sampled_reaction_index, a_hashIndex, true );
if( crossSectionSum >= sampleCrossSection ) break;
} }
else {
double temperatureFraction2 = ( a_temperature - m_temperatures[temperatureIndex1] ) / ( m_temperatures[temperatureIndex2] - m_temperatures[temperatureIndex1] );
double temperatureFraction1 = 1.0 - temperatureFraction2;
HeatedCrossSectionMultiGroup &heatedCrossSection2 = *m_heatedCrossSections[temperatureIndex2];
for( sampled_reaction_index = 0; sampled_reaction_index < numberOfReactions; ++sampled_reaction_index ) {
if( m_thresholds[sampled_reaction_index] >= a_energy ) continue;
crossSectionSum += temperatureFraction1 * heatedCrossSection1.reactionCrossSection( sampled_reaction_index, a_hashIndex, true );
crossSectionSum += temperatureFraction2 * heatedCrossSection2.reactionCrossSection( sampled_reaction_index, a_hashIndex, true );
if( crossSectionSum >= sampleCrossSection ) break;
}
}
if( sampled_reaction_index == numberOfReactions ) return( MCGIDI_nullReaction );
if( m_multiGroupThresholdIndex[sampled_reaction_index] == static_cast<int>( a_hashIndex ) ) {
double energyAboveThreshold = a_energy - m_thresholds[sampled_reaction_index];
if( energyAboveThreshold <= ( a_rng( ) * ( m_projectileMultiGroupBoundariesCollapsed[a_hashIndex+1] - m_thresholds[sampled_reaction_index] ) ) )
return( MCGIDI_nullReaction );
}
return( sampled_reaction_index );
}
// From file: MCGIDI_misc.cpp
/* *********************************************************************************************************//**
* This function returns a normalized Maxwellian speed (i.e., v = |velocity|) in 3d (i.e., x^2 Exp( -x^2 ))
* where v = sqrt(2 * T / m) * x.
* Using formula in https://link.springer.com/content/pdf/10.1007%2Fs10955-011-0364-y.pdf.
* Author Nader M.A. Mohamed, title "Efficient Algorithm for Generating Maxwell Random Variables".
*
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
*
* @return The sampled normalized Maxwellian speed.
***********************************************************************************************************/
template <typename RNG>
inline LUPI_HOST_DEVICE double sampleBetaFromMaxwellian( RNG && a_rng ) {
double _g = 2.0 / ( 1.37 * 0.5 * 1.772453850905516 ); // 1.772453850905516 = sqrt( pi ).
double beta, r1;
do {
r1 = a_rng( );
beta = sqrt( -2.0 * log( r1 ) );
} while( _g * r1 * beta < a_rng( ) );
return( beta );
}
namespace MCGIDI {
/* *********************************************************************************************************//**
* This function boost a particle from one frame to another frame. The frames have a relative speed *a_boostSpeed*
* and cosine of angle *a_boostMu* between their z-axes. BRB FIXME, currently it is the x-axis.
*
* @param a_input [in] Instance containing a random number generator that returns a double in the range [0, 1).
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_product [in] The particle to boost.
***********************************************************************************************************/
template <typename RNG>
inline LUPI_HOST_DEVICE void upScatterModelABoostParticle( Sampling::Input &a_input, RNG && a_rng, Sampling::Product &a_product ) {
double C_rel = 1.0;
if( a_input.m_relativeBeta != 0.0 ) {
C_rel = ( a_input.m_projectileBeta - a_input.m_muLab * a_input.m_targetBeta ) / a_input.m_relativeBeta;
if( C_rel > 1.0 ) C_rel = 1.0; // Handle round-off issue. Probably should check how big the issue is.
if( C_rel < -1.0 ) C_rel = -1.0; // Handle round-off issue. Probably should check how big the issue is.
}
double S_rel = sqrt( 1.0 - C_rel * C_rel );
double pz_vz = a_product.m_pz_vz;
a_product.m_pz_vz = C_rel * a_product.m_pz_vz + S_rel * a_product.m_px_vx;
a_product.m_px_vx = -S_rel * pz_vz + C_rel * a_product.m_px_vx;
double targetSpeed = MCGIDI_speedOfLight_cm_sec * a_input.m_targetBeta;
a_product.m_pz_vz += a_input.m_muLab * targetSpeed;
a_product.m_px_vx += sqrt( 1.0 - a_input.m_muLab * a_input.m_muLab ) * targetSpeed;
double phi = 2.0 * M_PI * a_rng( );
double sine = sin( phi );
double cosine = cos( phi );
double px_vx = a_product.m_px_vx;
a_product.m_px_vx = cosine * a_product.m_px_vx - sine * a_product.m_py_vy;
a_product.m_py_vy = sine * px_vx + cosine * a_product.m_py_vy;
double speed2 = a_product.m_px_vx * a_product.m_px_vx + a_product.m_py_vy * a_product.m_py_vy + a_product.m_pz_vz * a_product.m_pz_vz;
speed2 /= MCGIDI_speedOfLight_cm_sec * MCGIDI_speedOfLight_cm_sec;
a_product.m_kineticEnergy = particleKineticEnergyFromBeta2( a_product.m_productMass, speed2 );
}
}
// From file: MCGIDI_outputChannel.cpp
/* *********************************************************************************************************//**
* This method adds sampled products to *a_products*.
*
* @param a_protare [in] The Protare this Reaction belongs to.
* @param a_projectileEnergy [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_products [in] The object to add all sampled products to.
***********************************************************************************************************/
template <typename RNG, typename PUSHBACK>
LUPI_HOST_DEVICE void MCGIDI::OutputChannel::sampleProducts( ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input,
RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) const {
if( m_hasFinalStatePhotons ) {
double random = a_rng( );
double cumulative = 0.0;
bool sampled = false;
for( auto productIter = m_products.begin( ); productIter != m_products.end( ); ++productIter ) {
cumulative += (*productIter)->multiplicity( )->evaluate( a_projectileEnergy );
if( cumulative >= random ) {
(*productIter)->sampleFinalState( a_protare, a_projectileEnergy, a_input, a_rng, a_push_back, a_products );
sampled = true;
break;
}
}
if( !sampled ) { // BRB: FIXME: still need to code for continuum photon.
} }
else {
for( Vector<Product *>::const_iterator iter = m_products.begin( ); iter != m_products.end( ); ++iter )
(*iter)->sampleProducts( a_protare, a_projectileEnergy, a_input, a_rng, a_push_back, a_products );
}
if( m_totalDelayedNeutronMultiplicity != nullptr ) {
double totalDelayedNeutronMultiplicity = m_totalDelayedNeutronMultiplicity->evaluate( a_projectileEnergy );
if( a_rng( ) < totalDelayedNeutronMultiplicity ) { // Assumes that totalDelayedNeutronMultiplicity < 1.0, which it is.
double sum = 0.0;
totalDelayedNeutronMultiplicity *= a_rng( );
for( std::size_t i1 = 0; i1 < (std::size_t) m_delayedNeutrons.size( ); ++i1 ) {
DelayedNeutron const *delayedNeutron1( delayedNeutron( i1 ) );
Product const &product = delayedNeutron1->product( );
sum += product.multiplicity( )->evaluate( a_projectileEnergy );
if( sum >= totalDelayedNeutronMultiplicity ) {
product.distribution( )->sample( a_projectileEnergy, a_input, a_rng );
a_input.m_delayedNeutronIndex = delayedNeutron1->delayedNeutronIndex( );
a_input.m_delayedNeutronDecayRate = delayedNeutron1->rate( );
a_products.add( a_projectileEnergy, product.intid( ), product.index( ), product.userParticleIndex( ), product.mass( ),
a_input, a_rng, a_push_back, false );
break;
}
}
}
}
}
/* *********************************************************************************************************//**
* Returns the probability for a project with energy *a_energy_in* to cause this channel to emitted a particle of index
* *a_index* at angle *a_mu_lab* as seen in the lab frame. If a particle is emitted, *a_energy_out* is its sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_index [in] The index of the particle to emit.
* @param a_temperature [in] Specifies the temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_weight [in] The probability of emitting outgoing particle into lab angle *a_mu_lab*.
* @param a_energy_out [in] The energy of the emitted outgoing particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_cumulative_weight [in] The sum of the multiplicity for other outgoing particles with index *a_index*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::OutputChannel::angleBiasing( Reaction const *a_reaction, int a_index, double a_temperature, double a_energy_in, double a_mu_lab,
double &a_weight, double &a_energy_out, RNG && a_rng, double &a_cumulative_weight ) const {
for( Vector<Product *>::const_iterator iter = m_products.begin( ); iter != m_products.end( ); ++iter )
(*iter)->angleBiasing( a_reaction, a_index, a_temperature, a_energy_in, a_mu_lab, a_weight, a_energy_out, a_rng, a_cumulative_weight );
if( ( m_totalDelayedNeutronMultiplicity != nullptr ) && ( a_index == m_neutronIndex ) ) {
for( std::size_t i1 = 0; i1 < (std::size_t) m_delayedNeutrons.size( ); ++i1 ) {
DelayedNeutron const *delayedNeutron1( delayedNeutron( i1 ) );
Product const &product = delayedNeutron1->product( );
product.angleBiasing( a_reaction, a_index, a_temperature, a_energy_in, a_mu_lab, a_weight, a_energy_out, a_rng, a_cumulative_weight );
}
}
}
/* *********************************************************************************************************//**
* Returns the probability for a project with energy *a_energy_in* to cause this channel to emitted a particle of intid
* *a_intid* at angle *a_mu_lab* as seen in the lab frame. If a particle is emitted, *a_energy_out* is its sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_intid [in] The intid of the particle to emit.
* @param a_temperature [in] Specifies the temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_weight [in] The probability of emitting outgoing particle into lab angle *a_mu_lab*.
* @param a_energy_out [in] The energy of the emitted outgoing particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_cumulative_weight [in] The sum of the multiplicity for other outgoing particles with intid *a_intid*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::OutputChannel::angleBiasingViaIntid( Reaction const *a_reaction, int a_intid, double a_temperature, double a_energy_in, double a_mu_lab,
double &a_weight, double &a_energy_out, RNG && a_rng, double &a_cumulative_weight ) const {
for( Vector<Product *>::const_iterator iter = m_products.begin( ); iter != m_products.end( ); ++iter )
(*iter)->angleBiasingViaIntid( a_reaction, a_intid, a_temperature, a_energy_in, a_mu_lab, a_weight, a_energy_out, a_rng, a_cumulative_weight );
if( ( m_totalDelayedNeutronMultiplicity != nullptr ) && ( a_intid == PoPI::Intids::neutron ) ) {
for( std::size_t i1 = 0; i1 < (std::size_t) m_delayedNeutrons.size( ); ++i1 ) {
DelayedNeutron const *delayedNeutron1( delayedNeutron( i1 ) );
Product const &product = delayedNeutron1->product( );
product.angleBiasingViaIntid( a_reaction, a_intid, a_temperature, a_energy_in, a_mu_lab, a_weight, a_energy_out, a_rng, a_cumulative_weight );
}
}
}
// From file: MCGIDI_product.cpp
/* *********************************************************************************************************//**
* This method adds sampled products to *a_products*.
*
* @param a_protare [in] The Protare this Reaction belongs to.
* @param a_projectileEnergy [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_products [in] The object to add all sampled products to.
***********************************************************************************************************/
template <typename RNG, typename PUSHBACK>
LUPI_HOST_DEVICE void MCGIDI::Product::sampleProducts( ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input,
RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) const {
#ifdef MCGIDI_USE_OUTPUT_CHANNEL
if( m_outputChannel != nullptr ) {
m_outputChannel->sampleProducts( a_protare, a_projectileEnergy, a_input, a_rng, a_push_back, a_products ); }
else {
#endif
if( m_twoBodyOrder == TwoBodyOrder::secondParticle ) {
a_products.add( a_projectileEnergy, intid( ), index( ), userParticleIndex( ), mass( ), a_input, a_rng, a_push_back, m_intid == PoPI::Intids::photon ); }
else {
int _multiplicity = m_multiplicity->sampleBoundingInteger( a_projectileEnergy, a_rng );
int __multiplicity = _multiplicity;
for( ; _multiplicity > 0; --_multiplicity ) {
m_distribution->sample( a_projectileEnergy, a_input, a_rng );
a_input.m_delayedNeutronIndex = -1;
a_input.m_delayedNeutronDecayRate = 0.0;
a_products.add( a_projectileEnergy, intid( ), index( ), userParticleIndex( ), mass( ), a_input, a_rng, a_push_back, m_intid == PoPI::Intids::photon );
}
if( m_initialStateIndex >= 0 ) {
if( __multiplicity == 0 ) {
a_protare->sampleBranchingGammas( a_input, a_projectileEnergy, m_initialStateIndex, a_rng, a_push_back, a_products );
}
}
}
#ifdef MCGIDI_USE_OUTPUT_CHANNEL
}
#endif
}
/* *********************************************************************************************************//**
* This method adds sampled products to *a_products*. In particular, the product is a capture reaction
* primary gamma what has a finalState attribute. This gamma is added as well as the gammas from the
* gamma cascade.
*
* @param a_protare [in] The Protare this Reaction belongs to.
* @param a_projectileEnergy [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_products [in] The object to add all sampled products to.
***********************************************************************************************************/
template <typename RNG, typename PUSHBACK>
LUPI_HOST_DEVICE void MCGIDI::Product::sampleFinalState( ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input,
RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) const {
m_distribution->sample( a_projectileEnergy, a_input, a_rng );
a_input.m_delayedNeutronIndex = -1;
a_input.m_delayedNeutronDecayRate = 0.0;
a_products.add( a_projectileEnergy, m_intid, m_index, m_userParticleIndex, mass( ), a_input, a_rng, a_push_back, m_intid == PoPI::Intids::photon );
if( m_initialStateIndex >= 0 ) {
a_protare->sampleBranchingGammas( a_input, a_projectileEnergy, m_initialStateIndex, a_rng, a_push_back, a_products );
}
}
/* *********************************************************************************************************//**
* Returns the weight for a projectile with energy *a_energy_in* to cause this channel to emitted a particle of index
* *a_pid* at angle *a_mu_lab* as seen in the lab frame. If a particle is emitted, *a_energy_out* is its sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_index [in] The index of the particle to emit.
* @param a_temperature [in] Specifies the temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_weight [in] The weight of emitting outgoing particle into lab angle *a_mu_lab*.
* @param a_energy_out [in] The energy of the emitted outgoing particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_cumulative_weight [in] The sum of the multiplicity for other outgoing particles with index *a_index*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Product::angleBiasing( Reaction const *a_reaction, int a_index, double a_temperature, double a_energy_in, double a_mu_lab,
double &a_weight, double &a_energy_out, RNG && a_rng, double &a_cumulative_weight ) const {
#ifdef MCGIDI_USE_OUTPUT_CHANNEL
if( m_outputChannel != nullptr ) {
m_outputChannel->angleBiasing( a_reaction, a_index, a_temperature, a_energy_in, a_mu_lab, a_weight, a_energy_out, a_rng, a_cumulative_weight ); }
else {
#endif
if( m_index != a_index ) return;
double probability = 0.0;
double energy_out = 0.0;
if( a_cumulative_weight == 0.0 ) a_energy_out = 0.0;
if( m_multiplicity->type( ) == Function1dType::branching ) { // Needs to handle F1_Branching.
}
else {
probability = m_distribution->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab, a_rng, energy_out );
}
double weight = m_multiplicity->evaluate( a_energy_in ) * probability;
a_cumulative_weight += weight;
if( weight > a_rng( ) * a_cumulative_weight ) {
a_weight = weight;
a_energy_out = energy_out;
}
#ifdef MCGIDI_USE_OUTPUT_CHANNEL
}
#endif
}
/* *********************************************************************************************************//**
* Returns the weight for a projectile with energy *a_energy_in* to cause this channel to emitted a particle of intid
* *a_intid* at angle *a_mu_lab* as seen in the lab frame. If a particle is emitted, *a_energy_out* is its sampled outgoing energy.
*
* @param a_reaction [in] The reaction containing the particle which this distribution describes.
* @param a_intid [in] The intid of the particle to emit.
* @param a_temperature [in] Specifies the temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_weight [in] The weight of emitting outgoing particle into lab angle *a_mu_lab*.
* @param a_energy_out [in] The energy of the emitted outgoing particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_cumulative_weight [in] The sum of the multiplicity for other outgoing particles with intid *a_intid*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Product::angleBiasingViaIntid( Reaction const *a_reaction, int a_intid, double a_temperature, double a_energy_in, double a_mu_lab,
double &a_weight, double &a_energy_out, RNG && a_rng , double &a_cumulative_weight ) const {
#ifdef MCGIDI_USE_OUTPUT_CHANNEL
if( m_outputChannel != nullptr ) {
m_outputChannel->angleBiasingViaIntid( a_reaction, a_intid, a_temperature, a_energy_in, a_mu_lab, a_weight, a_energy_out, a_rng, a_cumulative_weight ); }
else {
#endif
if( m_intid != a_intid ) return;
angleBiasing( a_reaction, m_index, a_temperature, a_energy_in, a_mu_lab, a_weight, a_energy_out, a_rng, a_cumulative_weight );
#ifdef MCGIDI_USE_OUTPUT_CHANNEL
}
#endif
}
// From file: MCGIDI_protare.cpp
/* *********************************************************************************************************//**
* Samples a reaction of *this* and returns its index.
*
* @param a_input [in/out] Sample options requested by user.
* @param a_URR_protareInfos [in] URR information.
* @param a_hashIndex [in] Specifies the continuous energy hash index or multi-group index.
* @param a_crossSection [in] The total cross section for the protare at *a_input.temperature()* and *a_input.energy()*.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
*
* @return The index of the sampled reaction.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE std::size_t MCGIDI::Protare::sampleReaction( Sampling::Input &a_input, URR_protareInfos const &a_URR_protareInfos,
std::size_t a_hashIndex, double a_crossSection, RNG && a_rng ) const {
std::size_t reactionIndex = MCGIDI_nullReaction;
switch( protareType( ) ) {
case ProtareType::single:
reactionIndex = static_cast<ProtareSingle const *>( this )->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex,
a_crossSection, a_rng );
break;
case ProtareType::composite:
reactionIndex = static_cast<ProtareComposite const *>( this )->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex,
a_crossSection, a_rng );
break;
case ProtareType::TNSL:
reactionIndex = static_cast<ProtareTNSL const *>( this )->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex,
a_crossSection, a_rng );
break;
}
return( reactionIndex );
}
/* *********************************************************************************************************//**
* Samples gammas from a nuclide electro-magnetic decay.
*
* @param a_input [in] Sample options requested by user.
* @param a_projectileEnergy [in] The energy of the projectile.
* @param a_initialStateIndex [in] The index in *m_nuclideGammaBranchStateInfos* whose nuclide data are used for sampling.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_products [in] The object to add all sampled gammas to.
***********************************************************************************************************/
template <typename RNG, typename PUSHBACK>
LUPI_HOST_DEVICE void MCGIDI::ProtareSingle::sampleBranchingGammas( Sampling::Input &a_input, double a_projectileEnergy, int a_initialStateIndex,
RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) const {
int initialStateIndex = a_initialStateIndex;
double energyLevelSampleWidthUpper = 0.0; // Used for GRIN continuum levels to add variaction to outgoing photons.
NuclideGammaBranchStateInfo *nuclideGammaBranchStateInfo = nullptr;
if( initialStateIndex >= 0 ) nuclideGammaBranchStateInfo = m_nuclideGammaBranchStateInfos[static_cast<std::size_t>(initialStateIndex)];
while( initialStateIndex >= 0 ) {
auto const &branchIndices = nuclideGammaBranchStateInfo->branchIndices( );
double random = a_rng( );
double sum = 0.0;
initialStateIndex = -1; // Just in case the for loop never has "sum >= random".
for( std::size_t i1 = 0; i1 < branchIndices.size( ); ++i1 ) {
NuclideGammaBranchInfo *nuclideGammaBranchInfo = m_branches[branchIndices[i1]];
sum += nuclideGammaBranchInfo->probability( );
if( sum >= random ) {
double energyLevelSampleWidthLower = 0.0;
initialStateIndex = nuclideGammaBranchInfo->residualStateIndex( );
if( initialStateIndex >= 0 ) {
nuclideGammaBranchStateInfo = m_nuclideGammaBranchStateInfos[static_cast<std::size_t>(initialStateIndex)];
energyLevelSampleWidthLower = a_rng( ) * nuclideGammaBranchStateInfo->nuclearLevelEnergyWidth( );
}
if( nuclideGammaBranchInfo->photonEmissionProbability( ) > a_rng( ) ) {
a_input.setSampledType( Sampling::SampledType::photon );
a_input.m_dataInTargetFrame = false;
a_input.m_frame = GIDI::Frame::lab;
a_input.m_energyOut1 = nuclideGammaBranchInfo->gammaEnergy( ) + energyLevelSampleWidthUpper - energyLevelSampleWidthLower;
a_input.m_mu = 1.0 - 2.0 * a_rng( );
a_input.m_phi = 2.0 * M_PI * a_rng( );
a_products.add( a_projectileEnergy, PoPI::Intids::photon, m_photonIndex, userPhotonIndex( ), 0.0, a_input, a_rng, a_push_back, true );
}
energyLevelSampleWidthUpper = energyLevelSampleWidthLower;
break;
}
}
}
}
/* *********************************************************************************************************//**
* This function is used internally to sample a target's velocity (speed and cosine of angle relative to projectile)
* for a heated target using zero temperature, multi-grouped cross sections.
*
* @param a_input [in/out] Contains needed input like the targets temperature. Also will have the target sampled velocity on return if return value is *true*.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
*
* @return Returns *true* if target velocity is sampled and false otherwise.
***********************************************************************************************************/
template <typename RNG>
inline LUPI_HOST_DEVICE bool MCGIDI::ProtareSingle::sampleTargetBetaForUpscatterModelA( Sampling::Input &a_input, RNG && a_rng ) const {
double projectileBeta = MCGIDI_particleBeta( m_projectileMass, a_input.energy( ) );
double targetThermalBeta = MCGIDI_particleBeta( m_targetMass, a_input.temperature( ) );
a_input.m_projectileBeta = projectileBeta;
a_input.m_relativeBeta = projectileBeta;
a_input.m_muLab = 0.0;
a_input.m_targetBeta = 0.0;
if( targetThermalBeta < 1e-4 * projectileBeta ) return( false );
a_input.m_modelTemperature = 0.0;
double relativeBetaMin = projectileBeta - 2.0 * targetThermalBeta;
double relativeBetaMax = projectileBeta + 2.0 * targetThermalBeta;
std::size_t maxIndex = m_upscatterModelAGroupVelocities.size( ) - 2;
int intRelativeBetaMinIndex = binarySearchVector( relativeBetaMin, m_upscatterModelAGroupVelocities, true );
std::size_t relativeBetaMinIndex = static_cast<std::size_t>( intRelativeBetaMinIndex );
int intRelativeBetaMaxIndex = binarySearchVector( relativeBetaMax, m_upscatterModelAGroupVelocities, true );
std::size_t relativeBetaMaxIndex = static_cast<std::size_t>( intRelativeBetaMaxIndex );
double targetBeta, relativeBeta, mu;
if( relativeBetaMinIndex >= maxIndex ) relativeBetaMinIndex = maxIndex;
if( relativeBetaMaxIndex >= maxIndex ) relativeBetaMaxIndex = maxIndex;
if( relativeBetaMinIndex == relativeBetaMaxIndex ) {
targetBeta = targetThermalBeta * sampleBetaFromMaxwellian( a_rng );
mu = 1.0 - 2.0 * a_rng( );
relativeBeta = sqrt( targetBeta * targetBeta + projectileBeta * projectileBeta - 2.0 * mu * targetBeta * projectileBeta ); }
else {
double reactionRate;
double reactionRateMax = 0;
for( std::size_t i1 = relativeBetaMinIndex; i1 <= relativeBetaMaxIndex; ++i1 ) {
reactionRate = m_upscatterModelACrossSection[i1] * m_upscatterModelAGroupVelocities[i1+1];
if( reactionRate > reactionRateMax ) reactionRateMax = reactionRate;
}
do {
targetBeta = targetThermalBeta * sampleBetaFromMaxwellian( a_rng );
mu = 1.0 - 2.0 * a_rng( );
relativeBeta = sqrt( targetBeta * targetBeta + projectileBeta * projectileBeta - 2.0 * mu * targetBeta * projectileBeta );
std::size_t index = static_cast<std::size_t>( binarySearchVector( relativeBeta, m_upscatterModelAGroupVelocities, true ) );
if( index > maxIndex ) index = maxIndex;
reactionRate = m_upscatterModelACrossSection[index] * relativeBeta;
} while( reactionRate < a_rng( ) * reactionRateMax );
}
a_input.m_modelEnergy = particleKineticEnergy( m_projectileMass, relativeBeta );
a_input.m_relativeBeta = relativeBeta;
a_input.m_muLab = mu;
a_input.m_targetBeta = targetBeta;
return( true );
}
/* *********************************************************************************************************//**
* Returns the index of a sampled reaction for target temperature, projectile energy and total cross section
* as specified via argument *a_input*. Random numbers are obtained via *a_rng*.
*
* @param a_input [in/out] Sample options requested by user. The values m_modelTemperature and m_modelEnergy are set by this method.
* @param a_URR_protareInfos [in] URR information.
* @param a_hashIndex [in] Specifies the continuous energy hash index or multi-group index.
* @param a_crossSection [in] The total cross section for the protare at *a_input.temperature()* and *a_input.energy()*.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE std::size_t MCGIDI::ProtareSingle::sampleReaction( Sampling::Input &a_input, URR_protareInfos const &a_URR_protareInfos,
std::size_t a_hashIndex, double a_crossSection, RNG && a_rng ) const {
std::size_t hashIndex = a_hashIndex;
double crossSection1 = a_crossSection;
a_input.m_dataInTargetFrame = false;
a_input.m_modelTemperature = a_input.m_temperature;
a_input.m_modelEnergy = a_input.m_energy;
if( upscatterModelASupported( ) && ( a_input.m_upscatterModel == Sampling::Upscatter::Model::A ) ) {
a_input.m_dataInTargetFrame = sampleTargetBetaForUpscatterModelA( a_input, a_rng );
if( a_input.m_dataInTargetFrame ) {
if( m_continuousEnergy ) {
hashIndex = m_domainHash.index( a_input.m_modelEnergy ); }
else {
hashIndex = m_multiGroupHash.index( a_input.m_modelEnergy );
}
crossSection1 = crossSection( a_URR_protareInfos, hashIndex, a_input.m_modelTemperature, a_input.m_modelEnergy, true );
}
}
if( m_continuousEnergy ) return( m_heatedCrossSections.sampleReaction( a_URR_protareInfos, m_URR_index, hashIndex,
a_input.m_modelTemperature, a_input.m_modelEnergy, crossSection1, a_rng ) );
return( m_heatedMultigroupCrossSections.sampleReaction( hashIndex, a_input.m_modelTemperature, a_input.m_modelEnergy,
crossSection1, a_rng ) );
}
// From file: MCGIDI_protareComposite.cpp
/* *********************************************************************************************************//**
* Samples a reaction of *this* and returns its index.
*
* @param a_input [in] Sample options requested by user.
* @param a_URR_protareInfos [in] URR information.
* @param a_hashIndex [in] Specifies the continuous energy hash index or multi-group index.
* @param a_crossSection [in] The total cross section for the protare at *a_input.temperature()* and *a_input.energy()*.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
*
* @return The index of the sampled reaction.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE std::size_t MCGIDI::ProtareComposite::sampleReaction( Sampling::Input &a_input, URR_protareInfos const &a_URR_protareInfos,
std::size_t a_hashIndex, double a_crossSection, RNG && a_rng ) const {
std::size_t length = static_cast<std::size_t>( m_protares.size( ) );
std::size_t reaction_index = 0;
double cross_section_sum = 0.0;
double cross_section_rng = a_rng( ) * a_crossSection;
for( std::size_t i1 = 0; i1 < length; ++i1 ) {
double cross_section = m_protares[i1]->crossSection( a_URR_protareInfos, a_hashIndex, a_input.temperature( ), a_input.energy( ), true );
cross_section_sum += cross_section;
if( cross_section_sum > cross_section_rng ) {
std::size_t reaction_index2 = m_protares[i1]->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex, cross_section, a_rng );
reaction_index += reaction_index2;
if( reaction_index2 == MCGIDI_nullReaction ) reaction_index = MCGIDI_nullReaction;
break;
}
reaction_index += m_protares[i1]->numberOfReactions( );
}
return( reaction_index );
}
// From file: MCGIDI_protareTNSL.cpp
/* *********************************************************************************************************//**
* Returns the total cross section.
*
* @param a_input [in] Sample options requested by user.
* @param a_URR_protareInfos [in] URR information.
* @param a_hashIndex [in] Specifies the continuous energy hash index or multi-group index.
* @param a_crossSection [in] The total cross section for the protare at *a_input.temperature()* and *a_input.energy()*.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
*
* @return The index of the sampled reaction.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE std::size_t MCGIDI::ProtareTNSL::sampleReaction( Sampling::Input &a_input, URR_protareInfos const &a_URR_protareInfos,
std::size_t a_hashIndex, double a_crossSection, RNG && a_rng ) const {
std::size_t reactionIndex = 0;
if( ( a_input.energy( ) < m_TNSL_maximumEnergy ) && ( a_input.temperature( ) <= m_TNSL_maximumTemperature ) ) {
double TNSL_crossSection = m_TNSL->crossSection( a_URR_protareInfos, a_hashIndex, a_input.temperature( ), a_input.energy( ), true );
if( TNSL_crossSection > a_rng( ) * a_crossSection ) {
reactionIndex = m_TNSL->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex, TNSL_crossSection, a_rng ); }
else {
reactionIndex = m_protareWithoutElastic->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex,
a_crossSection - TNSL_crossSection, a_rng );
if( reactionIndex != MCGIDI_nullReaction ) reactionIndex += m_numberOfTNSLReactions + 1;
} }
else {
reactionIndex = m_protareWithElastic->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex, a_crossSection, a_rng );
if( reactionIndex != MCGIDI_nullReaction ) reactionIndex += m_numberOfTNSLReactions;
}
return( reactionIndex );
}
// From file: MCGIDI_reaction.cpp
/* *********************************************************************************************************//**
* This method adds sampled products to *a_products*.
*
* @param a_protare [in] The Protare this Reaction belongs to.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_products [in] The object to add all sampled products to.
* @param a_checkOrphanProducts [in] If true, associated orphan products are also sampled.
***********************************************************************************************************/
template <typename RNG, typename PUSHBACK>
LUPI_HOST_DEVICE void MCGIDI::Reaction::sampleProducts( Protare const *a_protare, Sampling::Input &a_input,
RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products, bool a_checkOrphanProducts ) const {
double projectileEnergy = a_input.modelEnergy( );
a_input.m_GRIN_intermediateResidual = -1;
a_input.m_reaction = this;
a_input.m_projectileMass = m_projectileMass;
a_input.m_targetMass = m_targetMass;
a_input.m_relativeBeta = MCGIDI_particleBeta( m_projectileMass, projectileEnergy );
if( m_GRIN_specialSampleProducts ) {
if( m_GRIN_capture != nullptr ) {
if( projectileEnergy < m_GRIN_maximumCaptureIncidentEnergy ) {
if( m_GRIN_capture->sampleProducts( (ProtareSingle const *) a_protare, projectileEnergy, a_input, a_rng, a_push_back, a_products ) ) {
return;
}
} }
else if( m_GRIN_inelastic != nullptr ) {
if( m_GRIN_inelastic->sampleProducts( (ProtareSingle const *) a_protare, projectileEnergy, a_input, a_rng, a_push_back, a_products ) ) {
return;
}
}
}
#ifdef MCGIDI_USE_OUTPUT_CHANNEL
m_outputChannel->sampleProducts( m_protareSingle, projectileEnergy, a_input, a_rng, a_push_back, a_products );
#else
if( m_hasFinalStatePhotons ) {
double random = a_rng( );
double cumulative = 0.0;
bool sampled = false;
for( auto productIter = m_products.begin( ); productIter != m_products.end( ); ++productIter ) {
cumulative += (*productIter)->multiplicity( )->evaluate( projectileEnergy );
if( cumulative >= random ) {
(*productIter)->sampleFinalState( m_protareSingle, projectileEnergy, a_input, a_rng, a_push_back, a_products );
sampled = true;
break;
}
}
if( !sampled ) { // BRB: FIXME: still need to code for continuum photon.
} }
else {
for( auto productIter = m_products.begin( ); productIter != m_products.end( ); ++productIter ) {
(*productIter)->sampleProducts( m_protareSingle, projectileEnergy, a_input, a_rng, a_push_back, a_products );
}
}
if( m_totalDelayedNeutronMultiplicity != nullptr ) {
double totalDelayedNeutronMultiplicity = m_totalDelayedNeutronMultiplicity->evaluate( projectileEnergy );
if( a_rng( ) < totalDelayedNeutronMultiplicity ) { // Assumes that totalDelayedNeutronMultiplicity < 1.0, which it is.
double sum = 0.0;
totalDelayedNeutronMultiplicity *= a_rng( );
for( std::size_t i1 = 0; i1 < (std::size_t) m_delayedNeutrons.size( ); ++i1 ) {
DelayedNeutron const *delayedNeutron1 = m_delayedNeutrons[i1];
Product const &product = delayedNeutron1->product( );
sum += product.multiplicity( )->evaluate( projectileEnergy );
if( sum >= totalDelayedNeutronMultiplicity ) {
product.distribution( )->sample( projectileEnergy, a_input, a_rng );
a_input.m_delayedNeutronIndex = delayedNeutron1->delayedNeutronIndex( );
a_input.m_delayedNeutronDecayRate = delayedNeutron1->rate( );
a_products.add( a_input.energy( ), product.intid( ), product.index( ), product.userParticleIndex( ), product.mass( ), a_input, a_rng, a_push_back, false );
break;
}
}
}
}
if( m_fissionResiduaIntid != -1 ) { // Special treatment to add 2 ENDL 99120 or 99125 products.
a_input.setSampledType( MCGIDI::Sampling::SampledType::unspecified );
a_input.m_frame = GIDI::Frame::lab;
a_input.m_energyOut1 = 0.0;
a_input.m_mu = 0.0;
a_input.m_phi = 0.0;
a_input.m_delayedNeutronIndex = -1;
a_input.m_delayedNeutronDecayRate = 0.0;
a_products.add( 0.0, m_fissionResiduaIntid, m_fissionResiduaIndex, m_fissionResiduaUserIndex, m_fissionResidualMass, a_input, a_rng, a_push_back, false );
a_products.add( 0.0, m_fissionResiduaIntid, m_fissionResiduaIndex, m_fissionResiduaUserIndex, m_fissionResidualMass, a_input, a_rng, a_push_back, false );
}
#endif
if( a_checkOrphanProducts ) {
for( auto productIter = m_associatedOrphanProducts.begin( ); productIter != m_associatedOrphanProducts.end( ); ++productIter ) {
(*productIter)->sampleProducts( m_protareSingle, projectileEnergy, a_input, a_rng, a_push_back, a_products );
}
}
}
/* *********************************************************************************************************//**
* This method adds sampled products to *a_products*.
*
* @param a_protare [in] The ProtareSingle this Reaction belongs to.
* @param a_projectileEnergy [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_products [in] The object to add all sampled products to.
***********************************************************************************************************/
template <typename RNG, typename PUSHBACK>
LUPI_HOST_DEVICE bool MCGIDI::GRIN_capture::sampleProducts( ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input,
RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) const {
std::size_t index = 0;
double random = a_rng( );
for( ; index < m_summedProbabilities.size( ) - 1; ++index ) {
if( random < m_summedProbabilities[index] ) break;
}
GRIN_captureLevelProbability *GRIN_captureLevelProbability1 = m_captureLevelProbabilities[index];
double availableEnergy = m_captureNeutronSeparationEnergy + a_projectileEnergy;
int primaryCaptureLevelIndex = GRIN_captureLevelProbability1->sampleCaptureLevel( a_protare, availableEnergy, a_rng );
NuclideGammaBranchStateInfo const *nuclideGammaBranchStateInfo = a_protare->nuclideGammaBranchStateInfos( )[static_cast<std::size_t>(primaryCaptureLevelIndex)];
a_input.m_GRIN_intermediateResidual = nuclideGammaBranchStateInfo->intid( );
a_input.setSampledType( Sampling::SampledType::photon );
a_input.m_dataInTargetFrame = false;
a_input.m_frame = GIDI::Frame::lab;
a_input.m_energyOut1 = availableEnergy - nuclideGammaBranchStateInfo->nuclearLevelEnergy( );
a_input.m_mu = 2 * a_rng( ) - 1.0;
a_input.m_phi = 2.0 * M_PI * a_rng( );
a_products.add( a_projectileEnergy, PoPI::Intids::photon, a_protare->photonIndex( ), a_protare->userPhotonIndex( ),
0.0, a_input, a_rng, a_push_back, true );
a_protare->sampleBranchingGammas( a_input, a_projectileEnergy, primaryCaptureLevelIndex, a_rng, a_push_back, a_products );
if( m_residualIntid != -1 ) {
a_input.m_energyOut1 = 0.0;
a_input.m_mu = 0.0;
a_input.m_phi = 0.0;
a_products.add( a_projectileEnergy, m_residualIntid, m_residualIndex, m_residualUserIndex, m_residualMass, a_input, a_rng, a_push_back, false );
}
return( true );
}
/* *********************************************************************************************************//**
* This method adds sampled products to *a_products*.
*
* @param a_protare [in] The ProtareSingle this Reaction belongs to.
* @param a_projectileEnergy [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_products [in] The object to add all sampled products to.
***********************************************************************************************************/
template <typename RNG, typename PUSHBACK>
LUPI_HOST_DEVICE bool MCGIDI::GRIN_inelastic::sampleProducts( ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input,
RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) const {
std::size_t index = 1;
for( ; index < m_energies.size( ); ++index ) {
if( m_energies[index] > a_projectileEnergy ) break;
}
--index;
GRIN_inelasticForEnergy *inelasticForEnergy = m_inelasticForEnergy[index];
int levelIndex = inelasticForEnergy->sampleLevelIndex( a_projectileEnergy, a_rng( ) );
if( levelIndex < 0 ) return( false );
NuclideGammaBranchStateInfo const *nuclideGammaBranchStateInfo = a_protare->nuclideGammaBranchStateInfos( )[static_cast<std::size_t>(levelIndex)];
a_input.m_GRIN_intermediateResidual = nuclideGammaBranchStateInfo->intid( );
double residualMass = m_targetMass + nuclideGammaBranchStateInfo->nuclearLevelEnergy( );
double initialMass = m_neutronMass + m_targetMass;
double finalMass = m_neutronMass + residualMass;
double twoBodyThreshold = 0.5 * ( finalMass * finalMass - initialMass * initialMass ) / m_targetMass;
double betaBoast = sqrt( a_projectileEnergy * ( a_projectileEnergy + 2. * m_neutronMass ) )
/ ( a_projectileEnergy + m_neutronMass + m_targetMass ); // betaBoast = v/c.
double _x = m_targetMass * ( a_projectileEnergy - twoBodyThreshold ) / ( finalMass * finalMass );
if( _x < 0 ) _x = 0.; // FIXME There needs to be a better test here.
double Kp;
if( _x < 2e-5 ) {
Kp = finalMass * _x * ( 1 - 0.5 * _x * ( 1 - _x ) ); }
else { // This is the relativistic formula derived from E^2 - (pc)^2 is frame independent.
Kp = sqrt( finalMass * finalMass + 2 * m_targetMass * ( a_projectileEnergy - twoBodyThreshold ) ) - finalMass;
}
if( Kp < 0 ) Kp = 0.; // FIXME There needs to be a better test here.
a_input.setSampledType( Sampling::SampledType::firstTwoBody );
a_input.m_mu = 1.0 - 2.0 * a_rng( );
a_input.m_phi = 2. * M_PI * a_rng( );
kinetics_COMKineticEnergy2LabEnergyAndMomentum( betaBoast, Kp, m_neutronMass, residualMass, a_input );
a_input.m_delayedNeutronIndex = -1;
a_input.m_delayedNeutronDecayRate = 0.0;
a_products.add( a_projectileEnergy, PoPI::Intids::neutron, m_neutronIndex, m_neutronUserParticleIndex, m_neutronMass, a_input, a_rng,
a_push_back, false );
a_products.add( a_projectileEnergy, m_targetIntid, m_targetIndex, m_targetUserParticleIndex,
m_targetMass, a_input, a_rng, a_push_back, false );
a_protare->sampleBranchingGammas( a_input, a_projectileEnergy, levelIndex, a_rng, a_push_back, a_products );
return( true );
}
/* *********************************************************************************************************//**
* This method samples a capture state level and returns an index into the a_protare->m_nuclideGammaBranchStateInfos vector
* of the sampled state level.
*
* @param a_energy [in] The neutron separation energy plus the projectile energy.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
*
* @return An integer of the sampled state in a_protare->m_nuclideGammaBranchStateInfos.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE int MCGIDI::GRIN_captureToCompound::sampleCaptureLevel( ProtareSingle const *a_protare, double a_energy,
RNG && a_rng, bool a_checkEnergy ) const {
if( a_checkEnergy ) {
NuclideGammaBranchStateInfo *nuclideGammaBranchStateInfo = a_protare->nuclideGammaBranchStateInfos( )[m_index];
if( nuclideGammaBranchStateInfo->nuclearLevelEnergy( ) < a_energy ) return( -1 );
}
double random = a_rng( );
std::size_t index = 0;
for( ; index < m_continuumIndices.m_levels.size( ) - 1; ++index ) {
if( m_continuumIndices.m_summedProbabilities[index] >= random ) break;
}
return( m_continuumIndices.m_levels[index] );
}
/* *********************************************************************************************************//**
* This method samples a capture state level and returns an index into the a_protare->m_nuclideGammaBranchStateInfos vector
* of the sampled state level.
*
* @param a_protare [in] The ProtareSingle this Reaction belongs to.
* @param a_energy [in] The neutron separation energy plus the projectile energy.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
*
* @return An integer of the sampled state in a_protare->m_nuclideGammaBranchStateInfos.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE int MCGIDI::GRIN_captureLevelProbability::sampleCaptureLevel( ProtareSingle const *a_protare, double a_energy, RNG && a_rng ) {
double random = a_rng( );
for( std::size_t index = 0; index < m_knownLevelsAndProbabilities.m_levels.size( ); ++index ) {
if( m_knownLevelsAndProbabilities.m_summedProbabilities[index] >= random ) {
return( m_knownLevelsAndProbabilities.m_levels[index] );
}
}
for( std::size_t i1 = 0; i1 < m_captureToCompounds.size( ) - 1; ++i1 ) {
GRIN_captureToCompound const *GRIN_captureToCompound1 = m_captureToCompounds[i1];
int index = GRIN_captureToCompound1->sampleCaptureLevel( a_protare, a_energy, a_rng, true );
if( index > -1 ) return( index );
}
return( m_captureToCompounds.back( )->sampleCaptureLevel( a_protare, a_energy, a_rng, false ) );
}
/* *********************************************************************************************************//**
* This method adds a null product to *a_products*. When running in multi-group mode, a sampled reaction may be rejected if the threshold
* is in the multi-group that the projectile is in. If this happens, only null products should be returned. This type of behavior was need
* in MCAPM but is probably not needed for MCGIDI.
*
* @param a_protare [in] The Protare this Reaction belongs to.
* @param a_projectileEnergy [in] The energy of the projectile.
* @param a_input [in] Sample options requested by user.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_products [in] The object to add all sampled products to.
***********************************************************************************************************/
template <typename RNG, typename PUSHBACK>
LUPI_HOST_DEVICE void MCGIDI::Reaction::sampleNullProducts( Protare const &a_protare, double a_projectileEnergy, Sampling::Input &a_input,
RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) {
a_input.m_GRIN_intermediateResidual = -1;
a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
a_input.m_dataInTargetFrame = false;
a_input.m_frame = GIDI::Frame::lab;
a_input.m_delayedNeutronIndex = -1;
a_input.m_delayedNeutronDecayRate = 0.0;
a_input.m_energyOut1 = a_projectileEnergy;
a_input.m_mu = 1.0;
a_input.m_phi = 0.0;
a_products.add( a_projectileEnergy, a_protare.projectileIntid( ), a_protare.projectileIndex( ), a_protare.projectileUserIndex( ), a_protare.projectileMass( ), a_input, a_rng, a_push_back, false );
}
/* *********************************************************************************************************//**
* Returns the weight for a project with energy *a_energy_in* to cause this reaction to emitted a particle of index
* *a_index* at angle *a_mu_lab* as seen in the lab frame. If a particle is emitted, *a_energy_out* is its sampled outgoing energy.
*
* @param a_index [in] The index of the particle to emit.
* @param a_temperature [in] Specifies the temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_energy_out [in] The energy of the emitted outgoing particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_cumulative_weight [in] The cumulative multiplicity.
* @param a_checkOrphanProducts [in] If true, associated orphan products are also sampled.
*
* @return The weight that the particle is emitted into mu *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Reaction::angleBiasing( int a_index, double a_temperature, double a_energy_in, double a_mu_lab, double &a_energy_out,
RNG && a_rng, double *a_cumulative_weight, bool a_checkOrphanProducts ) const {
double cumulative_weight1 = 0.0;
if( a_cumulative_weight == nullptr ) a_cumulative_weight = &cumulative_weight1;
double weight1 = 0.0;
#ifdef MCGIDI_USE_OUTPUT_CHANNEL
m_outputChannel->angleBiasing( this, a_index, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng, *a_cumulative_weight );
#else
for( auto productIter = m_products.begin( ); productIter != m_products.end( ); ++productIter ) {
(*productIter)->angleBiasing( this, a_index, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng, *a_cumulative_weight );
}
if( ( m_totalDelayedNeutronMultiplicity != nullptr ) && ( a_index == m_neutronIndex ) ) {
for( std::size_t i1 = 0; i1 < (std::size_t) m_delayedNeutrons.size( ); ++i1 ) {
DelayedNeutron const *delayedNeutron1 = m_delayedNeutrons[i1];
Product const &product = delayedNeutron1->product( );
product.angleBiasing( this, a_index, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng, *a_cumulative_weight );
}
}
#endif
if( a_checkOrphanProducts ) {
for( auto productIter = m_associatedOrphanProducts.begin( ); productIter != m_associatedOrphanProducts.end( ); ++productIter ) {
(*productIter)->angleBiasing( this, a_index, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng,
*a_cumulative_weight );
}
}
return( weight1 );
}
/* *********************************************************************************************************//**
* Returns the weight for a project with energy *a_energy_in* to cause this reaction to emitted a particle of intid
* *a_intid* at angle *a_mu_lab* as seen in the lab frame. If a particle is emitted, *a_energy_out* is its sampled outgoing energy.
*
* @param a_intid [in] The intid of the particle to emit.
* @param a_temperature [in] Specifies the temperature of the material.
* @param a_energy_in [in] The energy of the incident particle.
* @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
* @param a_energy_out [in] The energy of the emitted outgoing particle.
* @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
* @param a_cumulative_weight [in] The cumulative multiplicity.
* @param a_checkOrphanProducts [in] If true, associated orphan products are also sampled.
*
* @return The weight that the particle is emitted into mu *a_mu_lab*.
***********************************************************************************************************/
template <typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Reaction::angleBiasingViaIntid( int a_intid, double a_temperature, double a_energy_in, double a_mu_lab, double &a_energy_out,
RNG && a_rng, double *a_cumulative_weight, bool a_checkOrphanProducts ) const {
double cumulative_weight1 = 0.0;
if( a_cumulative_weight == nullptr ) a_cumulative_weight = &cumulative_weight1;
double weight1 = 0.0;
#ifdef MCGIDI_USE_OUTPUT_CHANNEL
m_outputChannel->angleBiasingViaIntid( this, a_intid, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng, *a_cumulative_weight );
#else
for( auto productIter = m_products.begin( ); productIter != m_products.end( ); ++productIter ) {
(*productIter)->angleBiasingViaIntid( this, a_intid, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng, *a_cumulative_weight );
}
if( ( m_totalDelayedNeutronMultiplicity != nullptr ) && ( a_intid == PoPI::Intids::neutron ) ) {
for( std::size_t i1 = 0; i1 < (std::size_t) m_delayedNeutrons.size( ); ++i1 ) {
DelayedNeutron const *delayedNeutron1 = m_delayedNeutrons[i1];
Product const &product = delayedNeutron1->product( );
product.angleBiasingViaIntid( this, a_intid, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng, *a_cumulative_weight );
}
}
#endif
if( a_checkOrphanProducts ) {
for( auto productIter = m_associatedOrphanProducts.begin( ); productIter != m_associatedOrphanProducts.end( ); ++productIter ) {
(*productIter)->angleBiasingViaIntid( this, a_intid, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng,
*a_cumulative_weight );
}
}
return( weight1 );
}
// From file: MCGIDI_sampling.cpp
template <typename RNG, typename PUSHBACK>
LUPI_HOST_DEVICE void MCGIDI::Sampling::ProductHandler::add( double a_projectileEnergy, int a_productIntid, int a_productIndex, int a_userProductIndex,
double a_productMass, Input &a_input, RNG && a_rng, PUSHBACK && a_push_back, bool a_isPhoton ) {
Product product;
if( a_isPhoton && ( a_input.m_sampledType != SampledType::unspecified ) ) a_input.m_sampledType = SampledType::photon;
product.m_sampledType = a_input.m_sampledType;
product.m_isVelocity = a_input.wantVelocity( );
product.m_productIntid = a_productIntid;
product.m_productIndex = a_productIndex;
product.m_userProductIndex = a_userProductIndex;
product.m_numberOfDBRC_rejections = a_input.m_numberOfDBRC_rejections;
product.m_productMass = a_productMass;
product.m_delayedNeutronIndex = a_input.m_delayedNeutronIndex;
product.m_delayedNeutronDecayRate = a_input.m_delayedNeutronDecayRate;
product.m_birthTimeSec = 0.;
if( product.m_delayedNeutronDecayRate > 0. ) {
product.m_birthTimeSec = -log( a_rng( ) ) / product.m_delayedNeutronDecayRate;
}
if( a_input.m_sampledType == SampledType::unspecified ) {
product.m_kineticEnergy = 0.0;
product.m_px_vx = 0.0;
product.m_py_vy = 0.0;
product.m_pz_vz = 0.0; }
else if( a_input.m_sampledType == SampledType::uncorrelatedBody ) {
if( a_input.m_frame == GIDI::Frame::centerOfMass ) {
a_input.m_frame = GIDI::Frame::lab;
double massRatio = a_input.m_projectileMass + a_input.m_targetMass;
massRatio = a_input.m_projectileMass * a_productMass / ( massRatio * massRatio );
double modifiedProjectileEnergy = massRatio * a_projectileEnergy;
double sqrtModifiedProjectileEnergy = sqrt( modifiedProjectileEnergy );
double sqrtEnergyOut_com = a_input.m_mu * sqrt( a_input.m_energyOut1 );
a_input.m_energyOut1 += modifiedProjectileEnergy + 2. * sqrtModifiedProjectileEnergy * sqrtEnergyOut_com;
if( a_input.m_energyOut1 != 0 ) a_input.m_mu = ( sqrtModifiedProjectileEnergy + sqrtEnergyOut_com ) / sqrt( a_input.m_energyOut1 );
}
product.m_kineticEnergy = a_input.m_energyOut1;
double p_v = sqrt( a_input.m_energyOut1 * ( a_input.m_energyOut1 + 2. * a_productMass ) );
if( product.m_isVelocity ) p_v *= MCGIDI_speedOfLight_cm_sec / ( a_input.m_energyOut1 + a_productMass );
product.m_pz_vz = p_v * a_input.m_mu;
p_v *= sqrt( 1. - a_input.m_mu * a_input.m_mu );
product.m_px_vx = p_v * sin( a_input.m_phi );
product.m_py_vy = p_v * cos( a_input.m_phi ); }
else if( a_input.m_sampledType == SampledType::firstTwoBody ) {
product.m_kineticEnergy = a_input.m_energyOut1;
product.m_px_vx = a_input.m_px_vx1;
product.m_py_vy = a_input.m_py_vy1;
product.m_pz_vz = a_input.m_pz_vz1;
a_input.m_sampledType = SampledType::secondTwoBody; }
else if( a_input.m_sampledType == SampledType::secondTwoBody ) {
product.m_kineticEnergy = a_input.m_energyOut2;
product.m_px_vx = a_input.m_px_vx2;
product.m_py_vy = a_input.m_py_vy2;
product.m_pz_vz = a_input.m_pz_vz2; }
else if( a_input.m_sampledType == SampledType::photon ) {
product.m_kineticEnergy = a_input.m_energyOut1;
double pz_vz_factor = a_input.m_energyOut1;
if( product.m_isVelocity ) pz_vz_factor = MCGIDI_speedOfLight_cm_sec;
product.m_pz_vz = a_input.m_mu * pz_vz_factor;
double v_perp = sqrt( 1.0 - a_input.m_mu * a_input.m_mu ) * pz_vz_factor;
product.m_px_vx = cos( a_input.m_phi ) * v_perp;
product.m_py_vy = sin( a_input.m_phi ) * v_perp; }
else {
product.m_kineticEnergy = a_input.m_energyOut2;
product.m_px_vx = a_input.m_px_vx2;
product.m_py_vy = a_input.m_py_vy2;
product.m_pz_vz = a_input.m_pz_vz2;
}
if( a_input.m_dataInTargetFrame && ( a_input.m_sampledType != SampledType::photon ) ) upScatterModelABoostParticle( a_input, a_rng, product );
a_push_back( product );
}
#endif // End of MCGIDI_headerSource_hpp_included
Reference in New Issue View Git Blame Copy Permalink