MCGIDI_headerSource.hpp File Reference

#include <climits>

Go to the source code of this file.

Namespaces
namespace	MCGIDI
	Simple C++ string class, useful as replacement for std::string if this cannot be used, or just for fun.

Macros
#define	crossSectionSumError   1e-6

Functions
template<typename RNG>
LUPI_HOST_DEVICE void	MCGIDI_sampleKleinNishina (double a_energyIn, RNG &&a_rng, double *a_energyOut, double *a_mu)
LUPI_HOST_DEVICE void	kinetics_COMKineticEnergy2LabEnergyAndMomentum (double a_beta, double a_kinetic_com, double a_m3cc, double a_m4cc, MCGIDI::Sampling::Input &a_input)
template<typename RNG>
LUPI_HOST_DEVICE double	sampleBetaFromMaxwellian (RNG &&a_rng)
template<typename RNG>
LUPI_HOST_DEVICE void	MCGIDI::upScatterModelABoostParticle (Sampling::Input &a_input, RNG &&a_rng, Sampling::Product &a_product)

Macro Definition Documentation

crossSectionSumError

#define crossSectionSumError   1e-6

Definition at line 19 of file MCGIDI_headerSource.hpp.

Referenced by MCGIDI::HeatedCrossSectionsContinuousEnergy::sampleReaction().

Function Documentation

kinetics_COMKineticEnergy2LabEnergyAndMomentum()

LUPI_HOST_DEVICE void kinetics_COMKineticEnergy2LabEnergyAndMomentum ( double a_beta, double a_kinetic_com, double a_m3cc, double a_m4cc, MCGIDI::Sampling::Input & a_input )

inline

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.

Parameters

a_beta	[in] The velocity/speedOflight of the com frame relative to the lab frame.
a_kinetic_com	[in] Total kinetic energy (K1 + K2) in the COM frame.
a_m3cc	[in] The mass of the first product.
a_m4cc	[in] The mass of the second product.
a_input	[in] Sample options requested by user and where the products' outgoing data are returned.

Definition at line 283 of file MCGIDI_headerSource.hpp.

                                                                     {
/*
    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;
    }
}

Referenced by MCGIDI::Distributions::AngularTwoBody::sample(), and MCGIDI::GRIN_inelastic::sampleProducts().

MCGIDI_sampleKleinNishina()

template<typename RNG>

LUPI_HOST_DEVICE void MCGIDI_sampleKleinNishina	(	double	a_energyIn,
		RNG &&	a_rng,
		double *	a_energyOut,
		double *	a_mu )

This function samples an energy and cosine of the angle for a photon for Klein Nishina scattering (i.e, incoherent photo-atomic scattering).

Parameters

a_energyIn	[in] The energy of the incoming photon.
a_rng	[in] The random number generator function that returns a double in the range [0, 1.0).
a_energyOut	[in] The energy of the scattered photon.
a_mu	[in] The cosine of the angle of the scattered photon's z-axis and the incoming photon's z-axis.

Definition at line 58 of file MCGIDI_headerSource.hpp.

                                                                                                                      {
/*
  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;
}

Referenced by MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::sample(), and MCGIDI::Distributions::IncoherentPhotoAtomicScattering::sample().

sampleBetaFromMaxwellian()

template<typename RNG>

LUPI_HOST_DEVICE double sampleBetaFromMaxwellian ( RNG && a_rng )

inline

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".

Parameters

a_rng

[in] The random number generator function that returns a double in the range [0, 1.0).

Returns

The sampled normalized Maxwellian speed.

Definition at line 2333 of file MCGIDI_headerSource.hpp.

                                                                        {
 
    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 );
}

Referenced by MCGIDI::ProtareSingle::sampleTargetBetaForUpscatterModelA().