EIC code displayed by LXR master/​include/​Geant4/​MCGIDI_headerSource.hpp	Source navigation Diff markup Identifier search General search
	Version: master test Revert   Change

Warning, file /include/Geant4/MCGIDI_headerSource.hpp was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

 /*
 # <<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