Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering Class Reference

#include <MCGIDI_distributions.hpp>

Inheritance diagram for MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering:

Public Member Functions

LUPI_HOST_DEVICE IncoherentBoundToFreePhotoAtomicScattering ()
LUPI_HOST IncoherentBoundToFreePhotoAtomicScattering (GIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering const &a_incoherentPhotoAtomicScattering, SetupInfo &a_setupInfo)
LUPI_HOST_DEVICE ~IncoherentBoundToFreePhotoAtomicScattering ()
LUPI_HOST_DEVICE double energyRatio (double a_energyIn, double a_mu) const
LUPI_HOST_DEVICE double evaluateKleinNishina (double a_energyIn, double a_mu) const
LUPI_HOST_DEVICE double evaluateOccupationNumber (double a_X, double a_mu) const
template<typename RNG>
LUPI_HOST_DEVICE void sample (double a_X, Sampling::Input &a_input, RNG &&a_rng) const
template<typename RNG>
LUPI_HOST_DEVICE double angleBiasing (Reaction const *a_reaction, double a_temperature, double a_energy_in, double a_mu_lab, RNG &&a_rng, double &a_energy_out) const
LUPI_HOST_DEVICE void serialize (LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
template<typename RNG>
LUPI_HOST_DEVICE double 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
Public Member Functions inherited from MCGIDI::Distributions::Distribution
LUPI_HOST_DEVICE Distribution ()
LUPI_HOST Distribution (Type a_type, GIDI::Distributions::Distribution const &a_distribution, SetupInfo &a_setupInfo)
LUPI_HOST Distribution (Type a_type, GIDI::Frame a_productFrame, SetupInfo &a_setupInfo)
LUPI_HOST_DEVICE MCGIDI_VIRTUAL_FUNCTION ~Distribution () MCGIDI_TRUE_VIRTUAL
LUPI_HOST_DEVICE Type type () const
LUPI_HOST_DEVICE GIDI::Frame productFrame () const
LUPI_HOST_DEVICE double projectileMass () const
LUPI_HOST_DEVICE double targetMass () const
LUPI_HOST_DEVICE double productMass () const
LUPI_HOST void setModelDBRC_data (Sampling::Upscatter::ModelDBRC_data *a_modelDBRC_data)
template<typename RNG>
LUPI_HOST_DEVICE MCGIDI_VIRTUAL_FUNCTION void sample (double a_X, Sampling::Input &a_input, RNG &&a_rng) const MCGIDI_TRUE_VIRTUAL
template<typename RNG>
LUPI_HOST_DEVICE MCGIDI_VIRTUAL_FUNCTION double angleBiasing (Reaction const *a_reaction, double a_temperature, double a_energy_in, double a_mu_lab, RNG &&a_rng, double &a_energy_out) const MCGIDI_TRUE_VIRTUAL
LUPI_HOST_DEVICE void serialize (LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
template<typename RNG>
LUPI_HOST_DEVICE void sample (double a_X, MCGIDI::Sampling::Input &a_input, RNG &&a_rng) const
template<typename RNG>
LUPI_HOST_DEVICE double angleBiasing (Reaction const *a_reaction, double a_temperature, double a_energy_in, double a_mu_lab, RNG &&a_rng, double &a_energy_out) const

Detailed Description

This class represents the distribution for an outgoing photon via incoherent photo-atomic elastic scattering.

Definition at line 309 of file MCGIDI_distributions.hpp.

Constructor & Destructor Documentation

◆ IncoherentBoundToFreePhotoAtomicScattering() [1/2]

LUPI_HOST_DEVICE MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::IncoherentBoundToFreePhotoAtomicScattering ( )

Default constructor used when broadcasting a Protare as needed by MPI or GPUs.

Definition at line 979 of file MCGIDI_distributions.cc.

979 {
980
981}

◆ IncoherentBoundToFreePhotoAtomicScattering() [2/2]

LUPI_HOST MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::IncoherentBoundToFreePhotoAtomicScattering ( GIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering const & a_incoherentBoundToFreePhotoAtomicScattering,
SetupInfo & a_setupInfo )
Parameters
a_incoherentBoundToFreePhotoAtomicScattering[in] The GIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering instance whose data is to be used to construct this.
a_setupInfo[in] Used internally when constructing a Protare to pass information to other constructors.

Definition at line 988 of file MCGIDI_distributions.cc.

990 :
991 Distribution( Type::incoherentBoundToFreePhotoAtomicScattering, a_incoherentBoundToFreePhotoAtomicScattering, a_setupInfo ),
992 m_bindingEnergy( 0.0 ) {
993
994 GIDI::ProtareSingle const &GIDI_protare = a_setupInfo.m_GIDI_protare;
995 auto monikers = GIDI_protare.styles( ).findAllOfMoniker( GIDI_MonteCarlo_cdfStyleChars );
996 std::string MonteCarlo_cdf = "";
997 if( monikers.size( ) == 1 ) {
998 MonteCarlo_cdf = monikers[0][0]->label( );
999 }
1000
1001 std::string Compton_href = a_incoherentBoundToFreePhotoAtomicScattering.href( );
1002
1003 if( Compton_href.find( MonteCarlo_cdf ) != std::string::npos ) {
1004 const GUPI::Ancestry *link = a_incoherentBoundToFreePhotoAtomicScattering.findInAncestry( Compton_href );
1005 GIDI::DoubleDifferentialCrossSection::IncoherentBoundToFreePhotoAtomicScattering const &dd = *static_cast<GIDI::DoubleDifferentialCrossSection::IncoherentBoundToFreePhotoAtomicScattering const *>( link );
1006 GIDI::Functions::Xs_pdf_cdf1d const *xpcCompton;
1007 GIDI::Functions::Function1dForm const *ComptonProfile = dd.ComptonProfile( );
1008 xpcCompton = static_cast<GIDI::Functions::Xs_pdf_cdf1d const *>( ComptonProfile );
1009
1010 std::vector<double> occupationNumbers = xpcCompton->cdf( );
1011 std::vector<double> pz_grid = xpcCompton->Xs( );
1012 std::size_t dataSize = pz_grid.size( );
1013 m_occupationNumber.resize( dataSize );
1014 m_pz.resize( dataSize );
1015 for( std::size_t index = 0; index < occupationNumbers.size( ); ++index ) {
1016 m_occupationNumber[index] = occupationNumbers[index];
1017 m_pz[index] = pz_grid[index];
1018 }
1019
1020 m_bindingEnergy = a_setupInfo.m_Q;
1021 }
1022}
GIDI_MonteCarlo_cdfStyleChars
#define GIDI_MonteCarlo_cdfStyleChars
Definition GIDI.hpp:245
GIDI::Functions::Xs_pdf_cdf1d::cdf
std::vector< double > const & cdf() const
Definition GIDI.hpp:1313
GIDI::Functions::Xs_pdf_cdf1d::Xs
std::vector< double > const & Xs() const
Definition GIDI.hpp:1311
GIDI::ProtareSingle::styles
Styles::Suite & styles()
Definition GIDI.hpp:4802
GIDI::Suite::findAllOfMoniker
std::vector< iterator > findAllOfMoniker(std::string const &a_moniker)
Definition GIDI_suite.cc:285

◆ ~IncoherentBoundToFreePhotoAtomicScattering()

LUPI_HOST_DEVICE MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::~IncoherentBoundToFreePhotoAtomicScattering ( )

Definition at line 1027 of file MCGIDI_distributions.cc.

1027 {
1028
1029}

Referenced by ~IncoherentBoundToFreePhotoAtomicScattering().

Member Function Documentation

◆ angleBiasing() [1/2]

template<typename RNG>
LUPI_HOST_DEVICE double MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::angleBiasing ( Reaction const * a_reaction,
double a_temperature,
double a_energy_in,
double a_mu_lab,
RNG && a_rng,
double & a_energy_out ) const

◆ angleBiasing() [2/2]

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

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.

Parameters
a_reaction[in] The reaction containing the particle which this distribution describes.
a_temperature[in] The temperature of the material.
a_energy_in[in] The energy of the incident particle.
a_mu_lab[in] The desired mu in the lab frame for the emitted particle.
a_userrng[in] A random number generator that takes the state a_rngState and returns a double in the range [0.0, 1.0).
a_rngState[in] The current state for the random number generator.
a_energy_out[in] The energy of the emitted outgoing particle.
Returns
The probability of emitting outgoing particle into lab angle a_mu_lab.

Definition at line 1248 of file MCGIDI_headerSource.hpp.

1249 {
1250
1251 URR_protareInfos URR_protareInfos1;
1252 double sigma = a_reaction->protareSingle( )->reactionCrossSection( a_reaction->reactionIndex( ), URR_protareInfos1, 0.0, a_energy_in );
1253
1254 double norm = M_PI * MCGIDI_classicalElectronRadius * MCGIDI_classicalElectronRadius / sigma;
1255
1256 double one_minus_mu = 1.0 - a_mu_lab;
1257 double alpha_in = a_energy_in / PoPI_electronMass_MeV_c2;
1258
1259 double quad_a, quad_b, quad_c, alpha_ratio, pz, occupation_pz, occupationNumberMax;
1260
1261 bool energetically_possible = false;
1262 int ep_it = 0;
1263 while( energetically_possible == false && ep_it < 1000 ){
1264
1265 // Sample electron momentum projection, pz
1266 occupationNumberMax = evaluateOccupationNumber( a_energy_in, -1.0 );
1267 occupation_pz = occupationNumberMax*a_rng();
1268 int intLowerIndex = binarySearchVector( occupation_pz, m_occupationNumber );
1269 pz = 0;
1270 if( intLowerIndex == -1 ){
1271 pz = m_pz.back();
1272 }
1273 else{
1274 std::size_t lowerIndex = static_cast<std::size_t>( intLowerIndex );
1275 pz = m_pz[lowerIndex] + (occupation_pz-m_occupationNumber[lowerIndex])*(m_pz[lowerIndex+1]-m_pz[lowerIndex])/(m_occupationNumber[lowerIndex+1]-m_occupationNumber[lowerIndex]);
1276 }
1277
1278 // Convert pz to outgoing photon energy
1279 alpha_ratio = energyRatio(a_energy_in, a_mu_lab);
1280 quad_a = pz*pz - (1/alpha_ratio)*(1/alpha_ratio);
1281 quad_b = -2*alpha_in*( pz*pz * a_mu_lab - (1/alpha_ratio));
1282 quad_c = alpha_in*alpha_in*( pz*pz - 1 );
1283 if(quad_b*quad_b - 4*quad_a*quad_c > 0){
1284 energetically_possible = true;
1285 }
1286 ep_it = ep_it + 1;
1287 }
1288
1289 const double quad_1 = -quad_b/(2*quad_a) + sqrt( quad_b*quad_b - 4*quad_a*quad_c )/( 2*quad_a );
1290 const double quad_2 = -quad_b/(2*quad_a) - sqrt( quad_b*quad_b - 4*quad_a*quad_c )/( 2*quad_a );
1291
1292 // Select the correct outgoing energy based on the pz value
1293 double alpha_out = 0;
1294 if(pz >= 0){
1295 if(quad_1 >= quad_2){
1296 alpha_out = quad_1;
1297 }
1298 else{
1299 alpha_out = quad_2;
1300 }
1301 }
1302 else{
1303 if(quad_1 >= quad_2){
1304 alpha_out = quad_2;
1305 }
1306 else{
1307 alpha_out = quad_1;
1308 }
1309 }
1310 alpha_ratio = alpha_out / alpha_in;
1311 a_energy_out = alpha_out * PoPI_electronMass_MeV_c2;
1312 double probability = evaluateOccupationNumber( a_energy_in, a_mu_lab );
1313 probability *= alpha_ratio * alpha_ratio * ( 1.0 + a_mu_lab * a_mu_lab + alpha_in * alpha_out * one_minus_mu * one_minus_mu ) * norm;
1314
1315 return( probability );
1316}
MCGIDI_classicalElectronRadius
#define MCGIDI_classicalElectronRadius
Definition MCGIDI.hpp:69
PoPI_electronMass_MeV_c2
#define PoPI_electronMass_MeV_c2
Definition PoPI.hpp:31
M_PI
#define M_PI
Definition SbMath.h:33
MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::energyRatio
LUPI_HOST_DEVICE double energyRatio(double a_energyIn, double a_mu) const
Definition MCGIDI_distributions.cc:1038
MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::evaluateOccupationNumber
LUPI_HOST_DEVICE double evaluateOccupationNumber(double a_X, double a_mu) const
Definition MCGIDI_distributions.cc:1074
MCGIDI::binarySearchVector
LUPI_HOST_DEVICE int binarySearchVector(double a_x, Vector< double > const &a_Xs, bool a_boundIndex=false)
Definition MCGIDI.hpp:318

◆ energyRatio()

LUPI_HOST_DEVICE double MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::energyRatio ( double a_energyIn,
double a_mu ) const

FIX ME.

Parameters
a_energyIn[in]
a_mu[in]

Definition at line 1038 of file MCGIDI_distributions.cc.

1038 {
1039
1040 double relativeEnergy = a_energyIn / PoPI_electronMass_MeV_c2;
1041
1042 return( 1.0 / ( 1.0 + relativeEnergy * ( 1.0 - a_mu ) ) );
1043}

Referenced by angleBiasing(), energyRatio(), evaluateKleinNishina(), and sample().

◆ evaluateKleinNishina()

LUPI_HOST_DEVICE double MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::evaluateKleinNishina ( double a_energyIn,
double a_mu ) const

This method evaluates the Klein-Nishina. FIX ME. This should be a function as it does not use member data.

Parameters
a_energyIn[in]
a_mu[in]

Definition at line 1052 of file MCGIDI_distributions.cc.

1052 {
1053
1054 double relativeEnergy = a_energyIn / PoPI_electronMass_MeV_c2;
1055 double _energyRatio = energyRatio( a_energyIn, a_mu );
1056 double one_minus_mu = 1.0 - a_mu;
1057
1058 double norm = ( 1.0 + 2.0 * relativeEnergy );
1059 norm = 2.0 * relativeEnergy * ( 2.0 + relativeEnergy * ( 1.0 + relativeEnergy ) * ( 8.0 + relativeEnergy ) ) / ( norm * norm );
1060 norm += ( ( relativeEnergy - 2.0 ) * relativeEnergy - 2.0 ) * log( 1.0 + 2.0 * relativeEnergy );
1061 norm /= relativeEnergy * relativeEnergy * relativeEnergy;
1062
1063 return( _energyRatio * _energyRatio * ( _energyRatio + a_mu * a_mu + relativeEnergy * one_minus_mu * one_minus_mu ) / norm );
1064}

Referenced by evaluateKleinNishina().

◆ evaluateOccupationNumber()

LUPI_HOST_DEVICE double MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::evaluateOccupationNumber ( double a_energyIn,
double a_mu ) const

This method evaluates the occupation number cdf.

Parameters
a_energyIn[in]
a_mu[in]

Definition at line 1074 of file MCGIDI_distributions.cc.

1074 {
1075
1076 const double alpha_in = a_energyIn / PoPI_electronMass_MeV_c2;
1077 //const double mec = 2.7309245307378233e-22; // m_e * c in SI units
1078 const double alpha_binding = -m_bindingEnergy/PoPI_electronMass_MeV_c2; // BE [MeV] / 0.511 [MeV]
1079 const double pzmax = ( -alpha_binding + alpha_in*(alpha_in - alpha_binding)*(1-a_mu) )/( sqrt( 2*alpha_in*(alpha_in-alpha_binding)*(1-a_mu) + alpha_binding*alpha_binding ) ); // *mec
1080
1081 int intLowerIndex = binarySearchVector( pzmax, m_pz );
1082 std::size_t lowerIndex = static_cast<std::size_t>( intLowerIndex );
1083 int size1 = static_cast<int>( m_occupationNumber.size( ) );
1084
1085 if( intLowerIndex == -1 || intLowerIndex == ( size1 - 1 ) ) {
1086 return( m_occupationNumber.back( ) );
1087 }
1088 if( intLowerIndex == -2 ){
1089 return( m_occupationNumber[0] );
1090 }
1091
1092 return(m_occupationNumber[lowerIndex] + (pzmax-m_pz[lowerIndex])*(m_occupationNumber[lowerIndex+1]-m_occupationNumber[lowerIndex])/(m_pz[lowerIndex+1]-m_pz[lowerIndex]) );
1093
1094}

Referenced by angleBiasing(), evaluateOccupationNumber(), and sample().

◆ sample()

template<typename RNG>
LUPI_HOST_DEVICE void MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::sample ( double a_X,
Sampling::Input & a_input,
RNG && a_rng ) const

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.

Parameters
a_X[in] The energy of the projectile.
a_input[in] Sample options requested by user.
a_userrng[in] A random number generator that takes the state a_rngState and returns a double in the range [0.0, 1.0).
a_rngState[in] The current state for the random number generator.

Definition at line 1124 of file MCGIDI_headerSource.hpp.

1124 {
1125
1126 double energyOut, mu, occupationNumber;
1127 // Convert incident photon energy [MeV] to units of rest mass energy of the electron
1128 const double alpha_in = a_X / PoPI_electronMass_MeV_c2;
1129 double alpha_ratio, occupation_pz, occupationNumberMax;
1130 double quad_a = 0, quad_b = 0, quad_c = 0, pz = 0; // Initialize with dummy values to silence compiler warnings
1131
1132 bool energetically_possible = false;
1133 int ep_it = 0;
1134 while( energetically_possible == false && ep_it < 1000 ){
1135 // Sample outgoing angle
1136 occupationNumberMax = evaluateOccupationNumber( a_X, -1.0 );
1137 if( a_X >= 10.0 ) { // This condition is not yet correct
1138 MCGIDI_sampleKleinNishina( alpha_in, a_rng, &energyOut, &mu ); }
1139 else {
1140 do {
1141 MCGIDI_sampleKleinNishina( alpha_in, a_rng, &energyOut, &mu );
1142 occupationNumber = evaluateOccupationNumber( a_X, mu );
1143 } while( occupationNumber < occupationNumberMax * a_rng( ) );
1144 }
1145
1146 // Sample electron momentum projection, pz
1147 occupation_pz = occupationNumberMax*a_rng();
1148 int intLowerIndex = binarySearchVector( occupation_pz, m_occupationNumber );
1149 if( intLowerIndex == -1 ){
1150 pz = m_pz.back();
1151 }
1152 else{
1153 std::size_t lowerIndex = static_cast<std::size_t>( intLowerIndex );
1154 pz = m_pz[lowerIndex] + (occupation_pz-m_occupationNumber[lowerIndex])*(m_pz[lowerIndex+1]-m_pz[lowerIndex])/(m_occupationNumber[lowerIndex+1]-m_occupationNumber[lowerIndex]);
1155 }
1156
1157 // Convert pz to outgoing photon energy
1158 alpha_ratio = energyRatio(a_X, mu);
1159 quad_a = pz*pz - (1/alpha_ratio)*(1/alpha_ratio);
1160 quad_b = -2*alpha_in*( pz*pz * mu - (1/alpha_ratio));
1161 quad_c = alpha_in*alpha_in*( pz*pz - 1 );
1162
1163 if(quad_b*quad_b - 4*quad_a*quad_c > 0){
1164 energetically_possible = true;
1165 }
1166 ep_it = ep_it + 1;
1167 }
1168
1169 const double quad_1 = -quad_b/(2*quad_a) + sqrt( quad_b*quad_b - 4*quad_a*quad_c )/( 2*quad_a );
1170 const double quad_2 = -quad_b/(2*quad_a) - sqrt( quad_b*quad_b - 4*quad_a*quad_c )/( 2*quad_a );
1171
1172 // Select the correct outgoing energy based on the pz value
1173 if(pz >= 0){
1174 if(quad_1 >= quad_2){
1175 energyOut = quad_1;
1176 }
1177 else{
1178 energyOut = quad_2;
1179 }
1180 }
1181 else{
1182 if(quad_1 >= quad_2){
1183 energyOut = quad_2;
1184 }
1185 else{
1186 energyOut = quad_1;
1187 }
1188 }
1189
1190 a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
1191 a_input.m_energyOut1 = energyOut * PoPI_electronMass_MeV_c2;
1192 a_input.m_mu = mu;
1193 a_input.m_phi = 2.0 * M_PI * a_rng( );
1194 a_input.m_frame = productFrame( );
1195}
MCGIDI_sampleKleinNishina
LUPI_HOST_DEVICE void MCGIDI_sampleKleinNishina(double a_energyIn, RNG &&a_rng, double *a_energyOut, double *a_mu)
Definition MCGIDI_headerSource.hpp:58
MCGIDI::Distributions::Distribution::productFrame
LUPI_HOST_DEVICE GIDI::Frame productFrame() const
Definition MCGIDI_distributions.hpp:44

◆ serialize()

LUPI_HOST_DEVICE void MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::serialize ( LUPI::DataBuffer & a_buffer,
LUPI::DataBuffer::Mode a_mode )

This method serializes this for broadcasting as needed for MPI and GPUs. The method can count the number of required bytes, pack this or unpack this depending on a_mode.

Parameters
a_buffer[in] The buffer to read or write data to depending on a_mode.
a_mode[in] Specifies the action of this method.

Definition at line 1104 of file MCGIDI_distributions.cc.

1104 {
1105
1106 Distribution::serialize( a_buffer, a_mode );
1107
1108 DATA_MEMBER_VECTOR_DOUBLE( m_pz, a_buffer, a_mode );
1109 DATA_MEMBER_VECTOR_DOUBLE( m_occupationNumber, a_buffer, a_mode );
1110
1111}
DATA_MEMBER_VECTOR_DOUBLE
#define DATA_MEMBER_VECTOR_DOUBLE(member, buf, mode)
Definition LUPI_dataBuffer.hpp:371
MCGIDI::Distributions::Distribution::serialize
LUPI_HOST_DEVICE void serialize(LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode)
Definition MCGIDI_distributions.cc:95

Referenced by serialize().

The documentation for this class was generated from the following files: