Geant4 11.4.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
MCGIDI_headerSource.hpp
Go to the documentation of this file.
1/*
2# <<BEGIN-copyright>>
3# Copyright 2019, Lawrence Livermore National Security, LLC.
4# This file is part of the gidiplus package (https://github.com/LLNL/gidiplus).
5# gidiplus is licensed under the MIT license (see https://opensource.org/licenses/MIT).
6# SPDX-License-Identifier: MIT
7# <<END-copyright>>
8*/
9
10#ifndef MCGIDI_headerSource_hpp_included
11#define MCGIDI_headerSource_hpp_included 1
12
13#include <climits>
14
15
16#ifdef MCGIDI_USE_DOUBLES
17 #define crossSectionSumError 1e-8
18#else
19 #define crossSectionSumError 1e-6
20#endif
21
22// From file: MCGIDI_URR.cpp
23
24/* *********************************************************************************************************//**
25 * Updates *this* if *a_protare* has a non-negative *URR_index*.
26 *
27 * @param a_protare [in] The protare whose *URR_index* is used to see if *this* needs updating.
28 * @param a_energy [in] The energy of the projectile.
29 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
30 ***********************************************************************************************************/
31
32template <typename RNG>
33LUPI_HOST_DEVICE void MCGIDI::URR_protareInfos::updateProtare( MCGIDI::Protare const *a_protare, double a_energy, RNG && a_rng ) {
34
35 for( std::size_t i1 = 0; i1 < a_protare->numberOfProtares( ); ++i1 ) {
36 ProtareSingle *protareSingle = const_cast<ProtareSingle *>( a_protare->protare( i1 ) );
37
38 if( protareSingle->URR_index( ) >= 0 ) {
39 URR_protareInfo &URR_protare_info = m_URR_protareInfos[static_cast<std::size_t>(protareSingle->URR_index())];
40
41 URR_protare_info.m_inURR = protareSingle->inURR( a_energy );
42 if( URR_protare_info.inURR( ) ) URR_protare_info.m_rng_Value = a_rng( );
43 }
44 }
45}
46
47
48/* *********************************************************************************************************//**
49 * This function samples an energy and cosine of the angle for a photon for Klein Nishina scattering (i.e, incoherent photo-atomic scattering).
50 *
51 * @param a_energyIn [in] The energy of the incoming photon.
52 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
53 * @param a_energyOut [in] The energy of the scattered photon.
54 * @param a_mu [in] The cosine of the angle of the scattered photon's z-axis and the incoming photon's z-axis.
55 ***********************************************************************************************************/
56
57template <typename RNG>
58LUPI_HOST_DEVICE void MCGIDI_sampleKleinNishina( double a_energyIn, RNG && a_rng, double *a_energyOut, double *a_mu ) {
59/*
60 Description
61 Sample the Klein-Nishina distribution.
62 The unit of energy is the rest mass of the electron.
63 Reference: R. N. Blomquist and E. N. Gelbard, Nuclear Science
64 and Engineering, 83, 380-384 (1983)
65
66 This routine was taken from MCAPM which was from MCNP with only cosmetic changes.
67
68 Input
69 a_energyIn - incident photon energy ( in electron rest mass units )
70 *rng - user supplied random number generator
71 Output
72 *a_energyOut - exiting photon energy ( in electron rest mass units )
73 *a_mu - exiting photon cosine
74*/
75
76 double a1, b1, t1, s1, r1, mu, energyOut;
77
78 a1 = 1.0 / a_energyIn;
79 b1 = 1.0 / ( 1.0 + 2.0 * a_energyIn );
80
81 if( a_energyIn < 3.0 ) { // Kahn''s method ( e < 1.5 MeV ) AECU-3259.
82 bool reject = true;
83
84 t1 = 1.0 / ( 1.0 + 8.0 * b1 );
85 do {
86 if( a_rng( ) <= t1 ) {
87 r1 = 2.0 * a_rng( );
88 s1 = 1.0 / ( 1.0 + a_energyIn * r1 );
89 mu = 1.0 - r1;
90 reject = a_rng( ) > 4.0 * s1 * ( 1.0 - s1 ); }
91 else {
92 s1 = ( 1.0 + 2.0 * a_energyIn * a_rng( ) ) * b1;
93 mu = 1.0 + a1 * ( 1.0 - 1.0 / s1 );
94 reject = a_rng( ) > 0.5 * ( mu * mu + s1 );
95 }
96 } while( reject );
97 energyOut = a_energyIn / ( 1 + a_energyIn * ( 1 - mu ) ); }
98 else { // Koblinger''s method ( e > 1.5 MeV ) NSE 56, 218 ( 1975 ).
99 t1 = a_rng( ) * ( 4.0 * a1 + 0.5 * ( 1.0 - b1 * b1 ) - ( 1.0 - 2.0 * ( 1.0 + a_energyIn ) * ( a1 * a1 ) ) * log( b1 ) );
100 if( t1 > 2.0 * a1 ) {
101 if( t1 > 4.0 * a1 ) {
102 if( t1 > 4.0 * a1 + 0.5 * ( 1.0 - b1 * b1 ) ) {
103 energyOut = a_energyIn * pow( b1, a_rng( ) );
104 mu = 1.0 + a1 - 1.0 / energyOut; }
105 else {
106 energyOut = a_energyIn * sqrt( 1.0 - a_rng( ) * ( 1.0 - b1 * b1 ) );
107 mu = 1.0 + a1 - 1.0 / energyOut;
108 } }
109 else {
110 energyOut = a_energyIn * ( 1.0 + a_rng( ) * ( b1 - 1.0 ) );
111 mu = 1.0 + a1 - 1.0 / energyOut; } }
112 else {
113 r1 = 2.0 * a_rng( );
114 mu = 1.0 - r1;
115 energyOut = 1.0 / ( a1 + r1 );
116 }
117 }
118
119 *a_mu = mu;
120 *a_energyOut = energyOut;
121
122 return;
123}
124
125/* *********************************************************************************************************//**
126 * This method samples the outgoing product data for the two outgoing particles in a two-body outgoing channel.
127 * First, is samples *mu*, the cosine of the product's outgoing angle, since this is for two-body interactions, *mu*
128 * is in the center-of-mass frame. It then calls kinetics_COMKineticEnergy2LabEnergyAndMomentum.
129 *
130 * @param a_X [in] The energy of the projectile in the lab frame.
131 * @param a_input [in] Sample options requested by user.
132 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
133 ***********************************************************************************************************/
134template <typename RNG>
135LUPI_HOST_DEVICE void MCGIDI::Distributions::Distribution::sample( double a_X, MCGIDI::Sampling::Input &a_input, RNG && a_rng ) const {
136
137 switch( type( ) ) {
138 case Distributions::Type::none:
139 break;
140 case Distributions::Type::unspecified:
141 static_cast<Distributions::Unspecified const *>( this )->sample( a_X, a_input, a_rng );
142 break;
143 case Distributions::Type::angularTwoBody:
144 static_cast<Distributions::AngularTwoBody const *>( this )->sample( a_X, a_input, a_rng );
145 break;
146 case Distributions::Type::KalbachMann:
147 static_cast<Distributions::KalbachMann const *>( this )->sample( a_X, a_input, a_rng );
148 break;
149 case Distributions::Type::uncorrelated:
150 static_cast<Distributions::Uncorrelated const *>( this )->sample( a_X, a_input, a_rng );
151 break;
152 case Distributions::Type::branching3d:
153 static_cast<Distributions::Branching3d const *>( this )->sample( a_X, a_input, a_rng );
154 break;
155 case Distributions::Type::energyAngularMC:
156 static_cast<Distributions::EnergyAngularMC const *>( this )->sample( a_X, a_input, a_rng );
157 break;
158 case Distributions::Type::angularEnergyMC:
159 static_cast<Distributions::AngularEnergyMC const *>( this )->sample( a_X, a_input, a_rng );
160 break;
161 case Distributions::Type::coherentPhotoAtomicScattering:
162 static_cast<Distributions::CoherentPhotoAtomicScattering const *>( this )->sample( a_X, a_input, a_rng );
163 break;
164 case Distributions::Type::incoherentPhotoAtomicScattering:
165 static_cast<Distributions::IncoherentPhotoAtomicScattering const *>( this )->sample( a_X, a_input, a_rng );
166 break;
167 case Distributions::Type::incoherentBoundToFreePhotoAtomicScattering:
168 static_cast<Distributions::IncoherentBoundToFreePhotoAtomicScattering const *>( this )->sample( a_X, a_input, a_rng );
169 break;
170 case Distributions::Type::incoherentPhotoAtomicScatteringElectron:
171 static_cast<Distributions::IncoherentPhotoAtomicScatteringElectron const *>( this )->sample( a_X, a_input, a_rng );
172 break;
173 case Distributions::Type::pairProductionGamma:
174 static_cast<Distributions::PairProductionGamma const *>( this )->sample( a_X, a_input, a_rng );
175 break;
176 case Distributions::Type::coherentElasticTNSL:
177 static_cast<Distributions::CoherentElasticTNSL const *>( this )->sample( a_X, a_input, a_rng );
178 break;
179 case Distributions::Type::incoherentElasticTNSL:
180 static_cast<Distributions::IncoherentElasticTNSL const *>( this )->sample( a_X, a_input, a_rng );
181 break;
182 }
183}
184
185// From file: MCGIDI_distributions.cpp
186
187/* *********************************************************************************************************//**
188 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
189 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
190 *
191 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
192 * @param a_temperature [in] The temperature of the material.
193 * @param a_energy_in [in] The energy of the incident particle.
194 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
195 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
196 * @param a_energy_out [out] The energy of the emitted outgoing particle.
197 *
198 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
199 ***********************************************************************************************************/
200
201template <typename RNG>
202LUPI_HOST_DEVICE double MCGIDI::Distributions::Distribution::angleBiasing( Reaction const *a_reaction, double a_temperature, double a_energy_in, double a_mu_lab,
203 RNG && a_rng, double &a_energy_out ) const {
204
205 double probability = 0.0;
206 a_energy_out = 0.0;
207
208 switch( type( ) ) {
209 case Distributions::Type::none:
210 break;
211 case Distributions::Type::unspecified:
212 probability = static_cast<Distributions::Unspecified const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
213 a_rng, a_energy_out );
214 break;
215 case Distributions::Type::angularTwoBody:
216 probability = static_cast<Distributions::AngularTwoBody const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
217 a_rng, a_energy_out );
218 break;
219 case Distributions::Type::KalbachMann:
220 probability = static_cast<Distributions::KalbachMann const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
221 a_rng, a_energy_out );
222 break;
223 case Distributions::Type::uncorrelated:
224 probability = static_cast<Distributions::Uncorrelated const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
225 a_rng, a_energy_out );
226 break;
227 case Distributions::Type::branching3d:
228 probability = static_cast<Distributions::Branching3d const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
229 a_rng, a_energy_out );
230 break;
231 case Distributions::Type::energyAngularMC:
232 probability = static_cast<Distributions::EnergyAngularMC const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
233 a_rng, a_energy_out );
234 break;
235 case Distributions::Type::angularEnergyMC:
236 probability = static_cast<Distributions::AngularEnergyMC const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
237 a_rng, a_energy_out );
238 break;
239 case Distributions::Type::coherentPhotoAtomicScattering:
240 probability = static_cast<Distributions::CoherentPhotoAtomicScattering const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
241 a_rng, a_energy_out );
242 break;
243 case Distributions::Type::incoherentPhotoAtomicScattering:
244 probability = static_cast<Distributions::IncoherentPhotoAtomicScattering const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
245 a_rng, a_energy_out );
246 break;
247 case Distributions::Type::incoherentBoundToFreePhotoAtomicScattering:
248 probability = static_cast<Distributions::IncoherentBoundToFreePhotoAtomicScattering const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
249 a_rng, a_energy_out );
250 break;
251 case Distributions::Type::incoherentPhotoAtomicScatteringElectron:
252 probability = static_cast<Distributions::IncoherentPhotoAtomicScatteringElectron const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
253 a_rng, a_energy_out );
254 break;
255 case Distributions::Type::pairProductionGamma:
256 probability = static_cast<Distributions::PairProductionGamma const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
257 a_rng, a_energy_out );
258 break;
259 case Distributions::Type::coherentElasticTNSL:
260 probability = static_cast<Distributions::CoherentElasticTNSL const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
261 a_rng, a_energy_out );
262 break;
263 case Distributions::Type::incoherentElasticTNSL:
264 probability = static_cast<Distributions::IncoherentElasticTNSL const *>( this )->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab,
265 a_rng, a_energy_out );
266 break;
267 }
268
269 return( probability );
270}
271
272/* *********************************************************************************************************//**
273 * This function calculates the products outgoing data (i.e., energy, velocity/momentum) for the two products of
274 * a two-body interaction give the cosine of the first product's outgoing angle.
275 *
276 * @param a_beta [in] The velocity/speedOflight of the com frame relative to the lab frame.
277 * @param a_kinetic_com [in] Total kinetic energy (K1 + K2) in the COM frame.
278 * @param a_m3cc [in] The mass of the first product.
279 * @param a_m4cc [in] The mass of the second product.
280 * @param a_input [in] Sample options requested by user and where the products' outgoing data are returned.
281 ***********************************************************************************************************/
282
283inline LUPI_HOST_DEVICE void kinetics_COMKineticEnergy2LabEnergyAndMomentum( double a_beta, double a_kinetic_com,
284 double a_m3cc, double a_m4cc, MCGIDI::Sampling::Input &a_input ) {
285/*
286 Relativity:
287 E = K + m, E^2 = K^2 + 2 K m + m^2, E^2 - m^2 = p^2 = K^2 + 2 K m
288
289 pc p v
290 ---- = v, --- = --- = beta = b
291 E E c
292
293 K ( K + 2 m )
294 b^2 = ---------------
295 ( K + m )^2
296*/
297 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;
298
299 p = sqrt( a_kinetic_com * ( a_kinetic_com + 2. * a_m3cc ) * ( a_kinetic_com + 2. * a_m4cc ) *
300 ( a_kinetic_com + 2. * ( a_m3cc + a_m4cc ) ) ) / ( 2. * ( a_kinetic_com + a_m3cc + a_m4cc ) );
301
302 py3 = p * sqrt( 1 - a_input.m_mu * a_input.m_mu );
303 px3 = py3 * cos( a_input.m_phi );
304 py3 *= sin( a_input.m_phi );
305 pz = p * a_input.m_mu;
306 if( 1 ) { // FIXME Assuming the answer is wanted in the lab frame for now.
307 a_input.m_frame = GIDI::Frame::lab;
308 E3 = sqrt( p * p + m3cc2 );
309 E4 = sqrt( p * p + m4cc2 );
310 gamma = sqrt( 1. / ( 1. - a_beta * a_beta ) );
311 pz3 = gamma * ( pz + a_beta * E3 );
312 pz4 = gamma * ( -pz + a_beta * E4 ); }
313 else { // COM frame.
314 a_input.m_frame = GIDI::Frame::centerOfMass;
315 pz3 = pz;
316 pz4 = -pz;
317 }
318
319 p_perp2 = px3 * px3 + py3 * py3;
320
321 a_input.m_px_vx1 = px3;
322 a_input.m_py_vy1 = py3;
323 a_input.m_pz_vz1 = pz3;
324 pp3 = p_perp2 + pz3 * pz3;
325 x = ( a_m3cc > 0 ) ? pp3 / ( 2 * m3cc2 ) : 1.;
326 if( x < 1e-5 ) {
327 a_input.m_energyOut1 = a_m3cc * x * ( 1 - 0.5 * x * ( 1 - x ) ); }
328 else {
329 a_input.m_energyOut1 = sqrt( m3cc2 + pp3 ) - a_m3cc;
330 }
331 a_input.m_px_vx2 = -px3;
332 a_input.m_py_vy2 = -py3;
333 a_input.m_pz_vz2 = pz4;
334 pp4 = p_perp2 + pz4 * pz4;
335 x = ( a_m4cc > 0 ) ? pp4 / ( 2 * m4cc2 ) : 1.;
336 if( x < 1e-5 ) {
337 a_input.m_energyOut2 = a_m4cc * x * ( 1 - 0.5 * x * ( 1 - x ) ); }
338 else {
339 a_input.m_energyOut2 = sqrt( m4cc2 + pp4 ) - a_m4cc;
340 }
341
342 if( a_input.wantVelocity( ) ) {
343 v_p = MCGIDI_speedOfLight_cm_sec / sqrt( pp3 + m3cc2 );
344 a_input.m_px_vx1 *= v_p;
345 a_input.m_py_vy1 *= v_p;
346 a_input.m_pz_vz1 *= v_p;
347
348 v_p = MCGIDI_speedOfLight_cm_sec / sqrt( pp4 + m4cc2 );
349 a_input.m_px_vx2 *= v_p;
350 a_input.m_py_vy2 *= v_p;
351 a_input.m_pz_vz2 *= v_p;
352 }
353}
354
355
356/* *********************************************************************************************************//**
357 * This method samples the outgoing product data for the two outgoing particles in a two-body outgoing channel.
358 * First, is samples *mu*, the cosine of the product's outgoing angle, since this is for two-body interactions, *mu*
359 * is in the center-of-mass frame. It then calls kinetics_COMKineticEnergy2LabEnergyAndMomentum.
360 *
361 * @param a_X [in] The energy of the projectile in the lab frame.
362 * @param a_input [in] Sample options requested by user.
363 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
364 ***********************************************************************************************************/
365template <typename RNG>
366LUPI_HOST_DEVICE void MCGIDI::Distributions::AngularTwoBody::sample( double a_X, MCGIDI::Sampling::Input &a_input, RNG && a_rng ) const {
367
368 double initialMass = projectileMass( ) + targetMass( ), finalMass = productMass( ) + m_residualMass;
369 double beta = sqrt( a_X * ( a_X + 2. * projectileMass( ) ) ) / ( a_X + initialMass ); // beta = v/c.
370 double _x = targetMass( ) * ( a_X - m_twoBodyThreshold ) / ( finalMass * finalMass );
371 double Kp; // Kp is the total kinetic energy for m3 and m4 in the COM frame.
372
373 a_input.setSampledType( Sampling::SampledType::firstTwoBody );
374
375 if( m_Upscatter ) {
376 if( ( a_input.m_upscatterModel == Sampling::Upscatter::Model::B ) || ( a_input.m_upscatterModel == Sampling::Upscatter::Model::BSnLimits )
377 || ( a_input.m_upscatterModel == Sampling::Upscatter::Model::DBRC ) ) {
378 if( upscatterModelB( a_X, a_input, a_rng ) ) return;
379 }
380 }
381
382 if( _x < 2e-5 ) {
383 Kp = finalMass * _x * ( 1 - 0.5 * _x * ( 1 - _x ) ); }
384 else { // This is the relativistic formula derived from E^2 - (pc)^2 is frame independent.
385 Kp = sqrt( finalMass * finalMass + 2 * targetMass( ) * ( a_X - m_twoBodyThreshold ) ) - finalMass;
386 }
387 if( Kp < 0 ) Kp = 0.; // FIXME There needs to be a better test here.
388
389 a_input.m_mu = m_angular->sample( a_X, a_rng( ), a_rng );
390 a_input.m_phi = 2. * M_PI * a_rng( );
391 kinetics_COMKineticEnergy2LabEnergyAndMomentum( beta, Kp, productMass( ), m_residualMass, a_input );
392}
393
394
395/* *********************************************************************************************************//**
396 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
397 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
398 *
399 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
400 * @param a_temperature [in] The temperature of the material.
401 * @param a_energy_in [in] The energy of the incident particle.
402 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
403 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
404 * @param a_energy_out [out] The energy of the emitted outgoing particle.
405 *
406 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
407 ***********************************************************************************************************/
408
409template <typename RNG>
410LUPI_HOST_DEVICE double MCGIDI::Distributions::AngularTwoBody::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
411 double a_energy_in, double a_mu_lab, RNG && a_rng, double &a_energy_out ) const {
412
413 a_energy_out = 0.0;
414
415 double initialMass = projectileMass( ) + targetMass( ), finalMass = productMass( ) + m_residualMass;
416 double _x = targetMass( ) * ( a_energy_in - m_twoBodyThreshold ) / ( finalMass * finalMass );
417 double Kp; // Total kinetic energy of products in the center-of-mass.
418
419 if( _x < 2e-5 ) {
420 Kp = finalMass * _x * ( 1 - 0.5 * _x * ( 1 - _x ) ); }
421 else { // This is the relativistic formula derived from E^2 - (pc)^2 which is frame independent (i.e., an invariant).
422 Kp = sqrt( finalMass * finalMass + 2.0 * targetMass( ) * ( a_energy_in - m_twoBodyThreshold ) ) - finalMass;
423 }
424 if( Kp < 0 ) Kp = 0.; // FIXME There needs to be a better test here.
425
426 double energy_product_com = 0.5 * Kp * ( Kp + 2.0 * m_residualMass ) / ( Kp + productMass( ) + m_residualMass );
427
428 if( productMass( ) == 0.0 ) {
429 double boostBeta = sqrt( a_energy_in * ( a_energy_in + 2. * projectileMass( ) ) ) / ( a_energy_in + initialMass ); // Good, even for projectileMass = 0.
430 double one_mu_beta = 1.0 - a_mu_lab * boostBeta;
431 double mu_com = ( a_mu_lab - boostBeta ) / one_mu_beta;
432 double Jacobian = ( 1.0 - boostBeta * boostBeta ) / ( one_mu_beta * one_mu_beta );
433
434 a_energy_out = sqrt( 1.0 - boostBeta * boostBeta ) * energy_product_com * ( 1.0 + mu_com * boostBeta );
435
436 return( Jacobian * m_angular->evaluate( a_energy_in, mu_com ) );
437 }
438
439 double productBeta = MCGIDI_particleBeta( productMass( ), energy_product_com );
440 double boostBeta = sqrt( a_energy_in * ( a_energy_in + 2. * projectileMass( ) ) ) / ( a_energy_in + initialMass ); // beta = v/c.
441 double muPlus = 0.0, JacobianPlus = 0.0, muMinus = 0.0, JacobianMinus = 0.0;
442
443 int numberOfMus = muCOM_From_muLab( a_mu_lab, boostBeta, productBeta, muPlus, JacobianPlus, muMinus, JacobianMinus );
444
445 if( numberOfMus == 0 ) return( 0.0 );
446
447 double probability = JacobianPlus * m_angular->evaluate( a_energy_in, muPlus );
448
449 if( numberOfMus == 2 ) {
450 double probabilityMinus = JacobianMinus * m_angular->evaluate( a_energy_in, muMinus );
451 probability += probabilityMinus;
452 if( probabilityMinus > a_rng( ) * probability ) {
453 muPlus = muMinus;
454 }
455 }
456
457 double productBeta2 = productBeta * productBeta;
458 double productBetaLab2 = productBeta2 + boostBeta * boostBeta * ( 1.0 - productBeta2 * ( 1.0 - muPlus * muPlus ) ) + 2.0 * muPlus * productBeta * boostBeta;
459 productBetaLab2 /= 1.0 - muPlus * productBeta * boostBeta;
460 a_energy_out = MCGIDI::particleKineticEnergyFromBeta2( productMass( ), productBetaLab2 );
461
462 return( probability );
463}
464
465/* *********************************************************************************************************//**
466 * This method samples a targets velocity for elastic upscattering for upscatter model B and then calculates the outgoing
467 * product data for the projectile and target.
468 *
469 * @param a_kineticLab [in] The kinetic energy of the projectile in the lab frame.
470 * @param a_input [in] Sample options requested by user.
471 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
472 ***********************************************************************************************************/
473
474template <typename RNG>
475LUPI_HOST_DEVICE bool MCGIDI::Distributions::AngularTwoBody::upscatterModelB( double a_kineticLab, Sampling::Input &a_input, RNG && a_rng ) const {
476
477 const double Two_sqrtPi = 1.1283791670955125739;
478 const double C0 = 1.0410423479, C1 = 3.9626339162e-4, C2 =-1.8654539193e-3, C3 = 1.0264818153e-4;
479 double neutronMass = projectileMass( ); // Mass are in incident energy unit / c**2.
480 double _targetMass = targetMass( );
481 double temperature = a_input.temperature( );
482 double kineticLabMax = 1e4 * temperature;
483
484 if( a_input.m_upscatterModel == Sampling::Upscatter::Model::BSnLimits ) {
485 double kineticLabMax200 = 200.0 * temperature;
486
487 kineticLabMax = 1e3 * temperature * neutronMass / _targetMass;
488 if( kineticLabMax < kineticLabMax200 ) kineticLabMax = kineticLabMax200;
489 if( a_kineticLab >= 0.1 ) kineticLabMax = 0.9 * a_kineticLab; }
490 else {
491 if( kineticLabMax > 1e-2 ) {
492 kineticLabMax = 1e-2;
493 if( kineticLabMax < 100.0 * temperature ) {
494 kineticLabMax = 100.0 * temperature;
495 if( kineticLabMax > 10.0 ) kineticLabMax = 10.0; // Assumes energy is in MeV.
496 }
497 }
498 }
499
500 if( a_kineticLab > kineticLabMax ) return( false ); // Only for low neutron energy.
501
502 a_input.m_frame = GIDI::Frame::lab;
503
504 double muProjectileTarget, relativeBeta, targetBeta;
505 double targetThermalBeta = MCGIDI_particleBeta( _targetMass, temperature );
506 double neutronBeta = MCGIDI_particleBeta( neutronMass, a_kineticLab );
507
508 a_input.m_numberOfDBRC_rejections = 0;
509 bool continueLoop = false;
510 double crossSectionMax = 0.0;
511 if( a_input.m_upscatterModel == Sampling::Upscatter::Model::DBRC ) {
512 double targetThermalSpeed = m_modelDBRC_data->targetThermalSpeed( temperature ); // Non-relativistic calculations, unlike targetThermalBeta.
513 crossSectionMax = m_modelDBRC_data->crossSectionMax( a_kineticLab, targetThermalSpeed );
514 }
515 do {
516 continueLoop = false;
517 ++a_input.m_numberOfDBRC_rejections;
518 do {
519 int MethodP1orP2 = 0; /* Assume P2 */
520 if( a_rng( ) * ( neutronBeta + Two_sqrtPi * targetThermalBeta ) < neutronBeta ) MethodP1orP2 = 1;
521 muProjectileTarget = 1.0 - 2.0 * a_rng( );
522 if( MethodP1orP2 == 0 ) { // x Exp( -x ) term.
523 targetBeta = targetThermalBeta * sqrt( -log( ( 1.0 - a_rng( ) ) * ( 1.0 - a_rng( ) ) ) ); }
524 else { // x^2 Exp( -x^2 ) term.
525 double x1;
526 do {
527 x1 = a_rng( );
528 x1 = sqrt( -log( ( 1.0 - a_rng( ) ) * ( 1.0 - x1 * x1 ) ) );
529 x1 = x1 / ( ( ( C3 * x1 + C2 ) * x1 + C1 ) * x1 + C0 );
530 } while( x1 > 4.0 );
531 targetBeta = targetThermalBeta * x1;
532 }
533 relativeBeta = sqrt( targetBeta * targetBeta + neutronBeta * neutronBeta - 2 * muProjectileTarget * targetBeta * neutronBeta );
534 } while( relativeBeta < ( targetBeta + neutronBeta ) * a_rng( ) );
535 if( a_input.m_upscatterModel == Sampling::Upscatter::Model::DBRC ) {
536 double relativeNeutronEnergy = 0.5 * productMass( ) * relativeBeta * relativeBeta;
537 if( m_modelDBRC_data->evaluate( relativeNeutronEnergy ) < a_rng( ) * crossSectionMax ) continueLoop = true;
538 }
539 } while( continueLoop );
540
541 double m1_12 = neutronMass / ( neutronMass + _targetMass );
542 double m2_12 = _targetMass / ( neutronMass + _targetMass );
543
544 double cosRelative = 0.0; // Cosine of angle between projectile velocity and relative velocity.
545 if( relativeBeta != 0.0 ) cosRelative = ( neutronBeta - muProjectileTarget * targetBeta ) / relativeBeta;
546 if( cosRelative > 1.0 ) {
547 cosRelative = 1.0; }
548 else if( cosRelative < -1.0 ) {
549 cosRelative = -1.0;
550 }
551 double sinRelative = sqrt( 1.0 - cosRelative * cosRelative ); // Sine of angle between projectile velocity and relative velocity.
552
553 a_input.m_muLab = muProjectileTarget;
554 a_input.m_targetBeta = targetBeta;
555 a_input.m_relativeBeta = relativeBeta;
556
557 double betaNeutronOut = m2_12 * relativeBeta;
558 double kineticEnergyRelative = particleKineticEnergy( neutronMass, betaNeutronOut );
559 double muCOM = m_angular->sample( kineticEnergyRelative, a_rng( ), a_rng );
560 double phiCOM = 2.0 * M_PI * a_rng( );
561 double SCcom = sqrt( 1.0 - muCOM * muCOM );
562 double SScom = SCcom * sin( phiCOM );
563 SCcom *= cos( phiCOM );
564
565 a_input.m_pz_vz1 = betaNeutronOut * ( muCOM * cosRelative - SCcom * sinRelative );
566 a_input.m_px_vx1 = betaNeutronOut * ( muCOM * sinRelative + SCcom * cosRelative );
567 a_input.m_py_vy1 = betaNeutronOut * SScom;
568
569 double massRatio = -neutronMass / _targetMass;
570 a_input.m_pz_vz2 = massRatio * a_input.m_pz_vz1;
571 a_input.m_px_vx2 = massRatio * a_input.m_px_vx1;
572 a_input.m_py_vy2 = massRatio * a_input.m_py_vy1;
573
574 double vCOMz = m1_12 * neutronBeta + m2_12 * muProjectileTarget * targetBeta; // Boost from center-of-mass to lab frame.
575 double vCOMx = m2_12 * sqrt( 1.0 - muProjectileTarget * muProjectileTarget ) * targetBeta;
576 a_input.m_pz_vz1 += vCOMz;
577 a_input.m_px_vx1 += vCOMx;
578 a_input.m_pz_vz2 += vCOMz;
579 a_input.m_px_vx2 += vCOMx;
580
581 double vx2_vy2 = a_input.m_px_vx1 * a_input.m_px_vx1 + a_input.m_py_vy1 * a_input.m_py_vy1;
582 double v2 = a_input.m_pz_vz1 * a_input.m_pz_vz1 + vx2_vy2;
583 a_input.m_mu = 0.0;
584 if( v2 != 0.0 ) a_input.m_mu = a_input.m_pz_vz1 / sqrt( v2 );
585 a_input.m_phi = atan2( a_input.m_py_vy1, a_input.m_px_vx1 );
586
587 a_input.m_energyOut1 = MCGIDI::particleKineticEnergyFromBeta2( neutronMass, v2 );
588 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 );
589
590 a_input.m_px_vx1 *= MCGIDI_speedOfLight_cm_sec;
591 a_input.m_py_vy1 *= MCGIDI_speedOfLight_cm_sec;
592 a_input.m_pz_vz1 *= MCGIDI_speedOfLight_cm_sec;
593
594 a_input.m_px_vx2 *= MCGIDI_speedOfLight_cm_sec;
595 a_input.m_py_vy2 *= MCGIDI_speedOfLight_cm_sec;
596 a_input.m_pz_vz2 *= MCGIDI_speedOfLight_cm_sec;
597
598 if( !a_input.wantVelocity( ) ) { // Return momenta.
599 a_input.m_px_vx1 *= neutronMass; // Non-relativistic.
600 a_input.m_py_vy1 *= neutronMass;
601 a_input.m_pz_vz1 *= neutronMass;
602
603 a_input.m_px_vx2 *= _targetMass;
604 a_input.m_py_vy2 *= _targetMass;
605 a_input.m_pz_vz2 *= _targetMass;
606 }
607
608 double phi = 2.0 * M_PI * a_rng( );
609 double sine = sin( phi );
610 double cosine = cos( phi );
611
612 double saved = a_input.m_px_vx1;
613 a_input.m_px_vx1 = cosine * a_input.m_px_vx1 - sine * a_input.m_py_vy1;
614 a_input.m_py_vy1 = sine * saved + cosine * a_input.m_py_vy1;
615
616 return( true );
617}
618
619
620/* *********************************************************************************************************//**
621 * This method samples the outgoing product data by sampling the outgoing energy E' and mu from the uncorrelated
622 * E and mu probabilities. It also samples the outgoing phi uniformly between 0 and 2 pi.
623 *
624 * @param a_X [in] The energy of the projectile.
625 * @param a_input [in] Sample options requested by user.
626 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
627 ***********************************************************************************************************/
628
629template <typename RNG>
630LUPI_HOST_DEVICE void MCGIDI::Distributions::Uncorrelated::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
631
632 a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
633 a_input.m_mu = m_angular->sample( a_X, a_rng( ), a_rng );
634 a_input.m_energyOut1 = m_energy->sample( a_X, a_rng( ), a_rng );
635 a_input.m_phi = 2. * M_PI * a_rng( );
636 a_input.m_frame = productFrame( );
637}
638
639/* *********************************************************************************************************//**
640 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
641 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
642 *
643 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
644 * @param a_temperature [in] The temperature of the material.
645 * @param a_energy_in [in] The energy of the incident particle.
646 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
647 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
648 * @param a_energy_out [in] The energy of the emitted outgoing particle.
649 *
650 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
651 ***********************************************************************************************************/
652
653template <typename RNG>
654LUPI_HOST_DEVICE double MCGIDI::Distributions::Uncorrelated::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
655 double a_energy_in, double a_mu_lab, RNG && a_rng, double &a_energy_out ) const {
656
657 if( productFrame( ) != GIDI::Frame::lab ) {
658 a_energy_out = 0.0;
659
660 double initialMass = projectileMass( ) + targetMass( );
661 double boostBeta = sqrt( a_energy_in * ( a_energy_in + 2. * projectileMass( ) ) ) / ( a_energy_in + initialMass ); // Good, even for projectileMass = 0.
662 double energy_out_com = m_energy->sample( a_energy_in, a_rng( ), a_rng );
663
664 if( productMass( ) == 0.0 ) {
665 double one_mu_beta = 1.0 - a_mu_lab * boostBeta;
666 double mu_com = ( a_mu_lab - boostBeta ) / one_mu_beta;
667 double Jacobian = ( 1.0 - boostBeta * boostBeta ) / ( one_mu_beta * one_mu_beta );
668
669 a_energy_out = sqrt( 1.0 - boostBeta * boostBeta ) * energy_out_com * ( 1.0 + mu_com * boostBeta );
670
671 return( Jacobian * m_angular->evaluate( a_energy_in, mu_com ) );
672 }
673
674 double productBeta = MCGIDI_particleBeta( productMass( ), energy_out_com );
675 double muPlus = 0.0, JacobianPlus = 0.0, muMinus = 0.0, JacobianMinus = 0.0;
676 int numberOfMus = muCOM_From_muLab( a_mu_lab, boostBeta, productBeta, muPlus, JacobianPlus, muMinus, JacobianMinus );
677
678 if( numberOfMus == 0 ) return( 0.0 );
679
680 double probability = JacobianPlus * m_angular->evaluate( a_energy_in, muPlus );
681
682 if( numberOfMus == 2 ) {
683 double probabilityMinus = JacobianMinus * m_angular->evaluate( a_energy_in, muMinus );
684
685 probability += probabilityMinus;
686 if( probabilityMinus > a_rng( ) * probability ) muPlus = muMinus;
687 }
688
689 double productBeta2 = productBeta * productBeta;
690 double productBetaLab2 = productBeta2 + boostBeta * boostBeta * ( 1.0 - productBeta2 * ( 1.0 - muPlus * muPlus ) ) + 2.0 * muPlus * productBeta * boostBeta;
691 productBetaLab2 /= 1.0 - muPlus * productBeta * boostBeta;
692 a_energy_out = MCGIDI::particleKineticEnergyFromBeta2( productMass( ), productBetaLab2 );
693
694 return( probability );
695 }
696
697 a_energy_out = m_energy->sample( a_energy_in, a_rng( ), a_rng );
698 return( m_angular->evaluate( a_energy_in, a_mu_lab ) );
699}
700
701/* *********************************************************************************************************//**
702 * This method samples the outgoing branching photons.
703 *
704 * @param a_X [in] The energy of the projectile in the lab frame.
705 * @param a_input [in] Sample options requested by user.
706 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
707 ***********************************************************************************************************/
708
709template <typename RNG>
713
714/* *********************************************************************************************************//**
715 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
716 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
717 *
718 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
719 * @param a_temperature [in] The temperature of the material.
720 * @param a_energy_in [in] The energy of the incident particle.
721 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
722 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
723 * @param a_energy_out [in] The energy of the emitted outgoing particle.
724 *
725 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
726 ***********************************************************************************************************/
727
728template <typename RNG>
729LUPI_HOST_DEVICE double MCGIDI::Distributions::Branching3d::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
730 LUPI_maybeUnused double a_energy_in, LUPI_maybeUnused double a_mu_lab, LUPI_maybeUnused RNG && a_rng, LUPI_maybeUnused double &a_energy_out ) const {
731
732 double probability = 0.0;
733
734 return( probability );
735}
736
737/* *********************************************************************************************************//**
738 * This method samples the outgoing product data by sampling the outgoing energy E' from the probability P(E'|E) and then samples mu from
739 * the probability P(mu|E,E'). It also samples the outgoing phi uniformly between 0 and 2 pi.
740 *
741 * @param a_X [in] The energy of the projectile.
742 * @param a_input [in] Sample options requested by user.
743 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
744 ***********************************************************************************************************/
745
746template <typename RNG>
747LUPI_HOST_DEVICE void MCGIDI::Distributions::EnergyAngularMC::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
748
749 double energyOut_1, energyOut_2;
750
751 a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
752 a_input.m_energyOut1 = m_energy->sample2dOf3d( a_X, a_rng( ), a_rng, &energyOut_1, &energyOut_2 );
753 a_input.m_mu = m_angularGivenEnergy->sample( a_X, energyOut_1, energyOut_2, a_rng( ), a_rng );
754 a_input.m_phi = 2. * M_PI * a_rng( );
755 a_input.m_frame = productFrame( );
756}
757
758/* *********************************************************************************************************//**
759 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
760 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
761 *
762 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
763 * @param a_temperature [in] The temperature of the material.
764 * @param a_energy_in [in] The energy of the incident particle.
765 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
766 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
767 * @param a_energy_out [in] The energy of the emitted outgoing particle.
768 *
769 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
770 ***********************************************************************************************************/
771
772template <typename RNG>
773LUPI_HOST_DEVICE double MCGIDI::Distributions::EnergyAngularMC::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
774 double a_energy_in, double a_mu_lab, RNG && a_rng, double &a_energy_out ) const {
775
776 double probability = 0.0;
777
778 if( productFrame( ) == GIDI::Frame::centerOfMass ) {
779 a_energy_out = m_energy->sample( a_energy_in, a_rng( ), a_rng );
780
781 double initialMass = projectileMass( ) + targetMass( );
782 double boostBeta = sqrt( a_energy_in * ( a_energy_in + 2. * projectileMass( ) ) ) / ( a_energy_in + initialMass ); // Good, even for projectileMass = 0.
783 double energy_out_com = m_energy->sample( a_energy_in, a_rng( ), a_rng );
784
785 if( productMass( ) == 0.0 ) {
786 double one_mu_beta = 1.0 - a_mu_lab * boostBeta;
787 double mu_com = ( a_mu_lab - boostBeta ) / one_mu_beta;
788 double Jacobian = ( 1.0 - boostBeta * boostBeta ) / ( one_mu_beta * one_mu_beta );
789
790 a_energy_out = sqrt( 1.0 - boostBeta * boostBeta ) * energy_out_com * ( 1.0 + mu_com * boostBeta );
791
792 return( Jacobian * m_angularGivenEnergy->evaluate( a_energy_in, energy_out_com, mu_com ) );
793 }
794
795 double productBeta = MCGIDI_particleBeta( productMass( ), energy_out_com );
796 double muPlus = 0.0, JacobianPlus = 0.0, muMinus = 0.0, JacobianMinus = 0.0;
797 int numberOfMus = muCOM_From_muLab( a_mu_lab, boostBeta, productBeta, muPlus, JacobianPlus, muMinus, JacobianMinus );
798
799 if( numberOfMus == 0 ) return( 0.0 );
800
801 probability = JacobianPlus * m_angularGivenEnergy->evaluate( a_energy_in, energy_out_com, muPlus );
802
803 if( numberOfMus == 2 ) {
804 double probabilityMinus = JacobianMinus * m_angularGivenEnergy->evaluate( a_energy_in, energy_out_com, muMinus );
805
806 probability += probabilityMinus;
807 if( probabilityMinus > a_rng( ) * probability ) muPlus = muMinus;
808 }
809
810 double productBeta2 = productBeta * productBeta;
811 double productBetaLab2 = productBeta2 + boostBeta * boostBeta * ( 1.0 - productBeta2 * ( 1.0 - muPlus * muPlus ) ) + 2.0 * muPlus * productBeta * boostBeta;
812 productBetaLab2 /= 1.0 - muPlus * productBeta * boostBeta;
813 a_energy_out = MCGIDI::particleKineticEnergyFromBeta2( productMass( ), productBetaLab2 ); }
814 else {
815 a_energy_out = m_energy->sample( a_energy_in, a_rng( ), a_rng );
816 probability = m_angularGivenEnergy->evaluate( a_energy_in, a_energy_out, a_mu_lab );
817 }
818
819 return( probability );
820}
821
822/* *********************************************************************************************************//**
823 * This method samples the outgoing product data by sampling the outgoing mu from the probability P(mu|E) and then samples E' from
824 * the probability P(E'|E,mu). It also samples the outgoing phi uniformly between 0 and 2 pi.
825 *
826 * @param a_X [in] The energy of the projectile.
827 * @param a_input [in] Sample options requested by user.
828 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
829 ***********************************************************************************************************/
830
831template <typename RNG>
832LUPI_HOST_DEVICE void MCGIDI::Distributions::AngularEnergyMC::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
833
834 double mu_1, mu_2;
835
836 a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
837 a_input.m_mu = m_angular->sample2dOf3d( a_X, a_rng( ), a_rng, &mu_1, &mu_2 );
838 a_input.m_energyOut1 = m_energyGivenAngular->sample( a_X, mu_1, mu_2, a_rng( ), a_rng );
839 a_input.m_phi = 2. * M_PI * a_rng( );
840 a_input.m_frame = productFrame( );
841}
842
843/* *********************************************************************************************************//**
844 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
845 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
846 *
847 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
848 * @param a_temperature [in] The temperature of the material.
849 * @param a_energy_in [in] The energy of the incident particle.
850 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
851 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
852 * @param a_energy_out [in] The energy of the emitted outgoing particle.
853 *
854 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
855 ***********************************************************************************************************/
856
857template <typename RNG>
858LUPI_HOST_DEVICE double MCGIDI::Distributions::AngularEnergyMC::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
859 double a_energy_in, double a_mu_lab, RNG && a_rng, double &a_energy_out ) const {
860
861 if( productFrame( ) != GIDI::Frame::lab ) LUPI_THROW( "AngularEnergyMC::angleBiasing: center-of-mass not supported." );
862
863 a_energy_out = m_energyGivenAngular->sample( a_energy_in, a_mu_lab, a_mu_lab, a_rng( ), a_rng );
864 return( m_angular->evaluate( a_energy_in, a_mu_lab ) );
865}
866
867
868/* *********************************************************************************************************//**
869 * This method samples the outgoing product data using the Kalbach-Mann formalism.
870 *
871 * @param a_X [in] The energy of the projectile.
872 * @param a_input [in] Sample options requested by user.
873 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
874 ***********************************************************************************************************/
875
876template <typename RNG>
877LUPI_HOST_DEVICE void MCGIDI::Distributions::KalbachMann::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
878
879 a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
880 a_input.m_energyOut1 = m_f->sample( a_X, a_rng( ), a_rng );
881 double rValue = m_r->evaluate( a_X, a_input.m_energyOut1 );
882 double aValue = m_a->evaluate( a_X, a_input.m_energyOut1 );
883
884 // In the following: Cosh[ a mu ] + r Sinh[ a mu ] = ( 1 - r ) Cosh[ a mu ] + r ( Cosh[ a mu ] + Sinh[ a mu ] ).
885 if( a_rng( ) >= rValue ) { // Sample the '( 1 - r ) Cosh[ a mu ]' term.
886 double T = ( 2. * a_rng( ) - 1. ) * sinh( aValue );
887
888 a_input.m_mu = log( T + sqrt( T * T + 1. ) ) / aValue; }
889 else { // Sample the 'r ( Cosh[ a mu ] + Sinh[ a mu ] )' term.
890 double rng1 = a_rng( ), exp_a = exp( aValue );
891
892 a_input.m_mu = log( rng1 * exp_a + ( 1. - rng1 ) / exp_a ) / aValue;
893 }
894 if( a_input.m_mu < -1 ) a_input.m_mu = -1;
895 if( a_input.m_mu > 1 ) a_input.m_mu = 1;
896
897 a_input.m_phi = 2. * M_PI * a_rng( );
898 a_input.m_frame = productFrame( );
899}
900
901/* *********************************************************************************************************//**
902 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
903 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
904 *
905 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
906 * @param a_temperature [in] The temperature of the material.
907 * @param a_energy_in [in] The energy of the incident particle.
908 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
909 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
910 * @param a_energy_out [in] The energy of the emitted outgoing particle.
911 *
912 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
913 ***********************************************************************************************************/
914
915template <typename RNG>
916LUPI_HOST_DEVICE double MCGIDI::Distributions::KalbachMann::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
917 double a_energy_in, double a_mu_lab, RNG && a_rng, double &a_energy_out ) const {
918
919 a_energy_out = 0.0;
920
921 double initialMass = projectileMass( ) + targetMass( );
922 double energy_out_com = m_f->sample( a_energy_in, a_rng( ), a_rng );
923 double productBeta = MCGIDI_particleBeta( productMass( ), energy_out_com );
924 double boostBeta = sqrt( a_energy_in * ( a_energy_in + 2. * projectileMass( ) ) ) / ( a_energy_in + initialMass ); // beta = v/c.
925
926 double muPlus = 0.0, JacobianPlus = 0.0, muMinus = 0.0, JacobianMinus = 0.0;
927
928 int numberOfMus = muCOM_From_muLab( a_mu_lab, boostBeta, productBeta, muPlus, JacobianPlus, muMinus, JacobianMinus );
929
930 if( numberOfMus == 0 ) return( 0.0 );
931
932 double rAtEnergyEnergyPrime = m_r->evaluate( a_energy_in, energy_out_com );
933 double aAtEnergyEnergyPrime = m_a->evaluate( a_energy_in, energy_out_com );
934 double aMu = aAtEnergyEnergyPrime * muPlus;
935
936 double probability = 0.5 * JacobianPlus;
937 if( productMass( ) == 0.0 ) {
938 probability *= 1.0 - rAtEnergyEnergyPrime + rAtEnergyEnergyPrime * aAtEnergyEnergyPrime * exp( aMu ) / sinh( aAtEnergyEnergyPrime ); }
939 else {
940 probability *= aAtEnergyEnergyPrime * ( cosh( aMu ) + rAtEnergyEnergyPrime * cosh( aMu ) ) / sinh( aAtEnergyEnergyPrime );
941 }
942
943 if( numberOfMus == 2 ) {
944 aMu = aAtEnergyEnergyPrime * muMinus;
945
946 double probabilityMinus = 0.5 * JacobianMinus;
947 if( productMass( ) == 0.0 ) {
948 probabilityMinus *= 1.0 - rAtEnergyEnergyPrime + rAtEnergyEnergyPrime * aAtEnergyEnergyPrime * exp( aMu ) / sinh( aAtEnergyEnergyPrime ); }
949 else {
950 probabilityMinus *= aAtEnergyEnergyPrime * ( cosh( aMu ) + rAtEnergyEnergyPrime * cosh( aMu ) ) / sinh( aAtEnergyEnergyPrime );
951 }
952 probability += probabilityMinus;
953
954 if( probabilityMinus > a_rng( ) * probability ) muPlus = muMinus;
955 }
956
957 double productBeta2 = productBeta * productBeta;
958 double productBetaLab2 = productBeta2 + boostBeta * boostBeta * ( 1.0 - productBeta2 * ( 1.0 - muPlus * muPlus ) ) + 2.0 * muPlus * productBeta * boostBeta;
959 productBetaLab2 /= 1.0 - muPlus * productBeta * boostBeta;
960 a_energy_out = MCGIDI::particleKineticEnergyFromBeta2( productMass( ), productBetaLab2 );
961
962 return( probability );
963}
964
965/* *********************************************************************************************************//**
966 * This method samples the outgoing product data from the coherent photo-atomic scattering law.
967 * It also samples the outgoing phi uniformly between 0 and 2 pi.
968 *
969 * @param a_X [in] The energy of the projectile in the lab frame.
970 * @param a_input [in] Sample options requested by user.
971 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
972 ***********************************************************************************************************/
973
974template <typename RNG>
975LUPI_HOST_DEVICE void MCGIDI::Distributions::CoherentPhotoAtomicScattering::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
976
977 a_input.m_energyOut1 = a_X;
978
979 int intLowerIndex = binarySearchVector( a_X, m_energies );
980
981 if( intLowerIndex < 1 ) {
982 do {
983 a_input.m_mu = 1.0 - 2.0 * a_rng( );
984 } while( ( 1.0 + a_input.m_mu * a_input.m_mu ) < 2.0 * a_rng( ) ); }
985 else {
986 std::size_t lowerIndex = static_cast<std::size_t>( intLowerIndex );
987 double _a = m_a[lowerIndex];
988 double X_i = m_energies[lowerIndex];
989 double formFactor_i = m_formFactor[lowerIndex];
990 double formFactor_X_i = formFactor_i * X_i;
991 double Z = a_X / X_i;
992 double realAnomalousFactor = 0.0;
993 double imaginaryAnomalousFactor = 0.0;
994
995 if( m_anomalousDataPresent ) {
996 realAnomalousFactor = m_realAnomalousFactor->evaluate( a_X );
997 imaginaryAnomalousFactor = m_imaginaryAnomalousFactor->evaluate( a_X );
998 }
999
1000 double anomalousFactorSquared = realAnomalousFactor * realAnomalousFactor + imaginaryAnomalousFactor * imaginaryAnomalousFactor;
1001 double normalization = m_integratedFormFactorSquared[lowerIndex] + formFactor_X_i * formFactor_X_i * Z_a( Z, 2.0 * _a + 2.0 );
1002
1003 anomalousFactorSquared = 0.0;
1004 if( anomalousFactorSquared != 0.0 ) {
1005 double integratedFormFactor_i = m_integratedFormFactor[lowerIndex] + formFactor_X_i * Z_a( Z, _a + 2.0 );
1006
1007 normalization += 2.0 * integratedFormFactor_i * realAnomalousFactor + 0.5 * anomalousFactorSquared * a_X * a_X;
1008 }
1009
1010 do {
1011 double partialIntegral = a_rng( ) * normalization;
1012 double X;
1013 if( anomalousFactorSquared == 0.0 ) {
1014 intLowerIndex = binarySearchVector( partialIntegral, m_integratedFormFactorSquared );
1015 lowerIndex = static_cast<std::size_t>( intLowerIndex );
1016
1017 if( lowerIndex == 0 ) {
1018 X = sqrt( 2.0 * partialIntegral ) / m_formFactor[0]; }
1019 else {
1020 double remainer = partialIntegral - m_integratedFormFactorSquared[lowerIndex];
1021 double epsilon = 2.0 * m_a[lowerIndex] + 2.0;
1022
1023 X_i = m_energies[lowerIndex];
1024 formFactor_i = m_formFactor[lowerIndex];
1025 formFactor_X_i = formFactor_i * X_i;
1026
1027 remainer /= formFactor_X_i * formFactor_X_i;
1028 if( fabs( epsilon ) < 1e-6 ) {
1029 X = X_i * exp( remainer ); }
1030 else {
1031 X = X_i * pow( 1.0 + epsilon * remainer, 1.0 / epsilon );
1032 }
1033 } }
1034 else { // Currently not implemented.
1035 X = 0.5 * a_X;
1036 }
1037 double X_E = X / a_X;
1038 a_input.m_mu = 1.0 - 2.0 * X_E * X_E;
1039 } while( ( 1.0 + a_input.m_mu * a_input.m_mu ) < 2.0 * a_rng( ) );
1040 }
1041
1042 a_input.m_phi = 2.0 * M_PI * a_rng( );
1043 a_input.m_frame = productFrame( );
1044}
1045
1046/* *********************************************************************************************************//**
1047 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
1048 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
1049 *
1050 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
1051 * @param a_temperature [in] The temperature of the material.
1052 * @param a_energy_in [in] The energy of the incident particle.
1053 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
1054 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1055 * @param a_energy_out [in] The energy of the emitted outgoing particle.
1056 *
1057 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
1058 ***********************************************************************************************************/
1059
1060template <typename RNG>
1061LUPI_HOST_DEVICE double MCGIDI::Distributions::CoherentPhotoAtomicScattering::angleBiasing( Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
1062 double a_energy_in, double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
1063
1064 a_energy_out = a_energy_in;
1065
1066 URR_protareInfos URR_protareInfos1;
1067 double sigma = a_reaction->protareSingle( )->reactionCrossSection( a_reaction->reactionIndex( ), URR_protareInfos1, 0.0, a_energy_in );
1068 double formFactor = evaluateFormFactor( a_energy_in, a_mu_lab );
1069 double imaginaryAnomalousFactor = 0.0;
1070
1071 if( m_anomalousDataPresent ) {
1072 formFactor += m_realAnomalousFactor->evaluate( a_energy_in );
1073 imaginaryAnomalousFactor = m_imaginaryAnomalousFactor->evaluate( a_energy_in );
1074 }
1075
1076 double probability = M_PI * MCGIDI_classicalElectronRadius * MCGIDI_classicalElectronRadius * ( 1.0 + a_mu_lab * a_mu_lab )
1077 * ( formFactor * formFactor + imaginaryAnomalousFactor * imaginaryAnomalousFactor ) / sigma;
1078
1079 return( probability );
1080}
1081
1082/* *********************************************************************************************************//**
1083 * This method samples the outgoing product data by sampling the outgoing energy E' from the probability P(E'|E) and then samples mu from
1084 * the probability P(mu|E,E'). It also samples the outgoing phi uniformly between 0 and 2 pi.
1085 *
1086 * @param a_X [in] The energy of the projectile.
1087 * @param a_input [in] Sample options requested by user.
1088 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1089 ***********************************************************************************************************/
1090
1091template <typename RNG>
1092LUPI_HOST_DEVICE void MCGIDI::Distributions::IncoherentPhotoAtomicScattering::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
1093
1094 double k1 = a_X / PoPI_electronMass_MeV_c2;
1095 double energyOut, mu, scatteringFactor;
1096
1097 if( a_X >= m_energies.back( ) ) {
1098 MCGIDI_sampleKleinNishina( k1, a_rng, &energyOut, &mu ); }
1099 else {
1100 double scatteringFactorMax = evaluateScatteringFactor( a_X );
1101 do {
1102 MCGIDI_sampleKleinNishina( k1, a_rng, &energyOut, &mu );
1103 scatteringFactor = evaluateScatteringFactor( a_X * sqrt( 0.5 * ( 1.0 - mu ) ) );
1104 } while( scatteringFactor < a_rng( ) * scatteringFactorMax );
1105 }
1106
1107 a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
1108 a_input.m_energyOut1 = energyOut * PoPI_electronMass_MeV_c2;
1109 a_input.m_mu = mu;
1110 a_input.m_phi = 2.0 * M_PI * a_rng( );
1111 a_input.m_frame = productFrame( );
1112}
1113
1114/* *********************************************************************************************************//**
1115 * This method samples the outgoing product data by sampling the outgoing energy E' from the probability P(E'|E) and then samples mu from
1116 * the probability P(mu|E,E'). It also samples the outgoing phi uniformly between 0 and 2 pi.
1117 *
1118 * @param a_X [in] The energy of the projectile.
1119 * @param a_input [in] Sample options requested by user.
1120 * @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).
1121 * @param a_rngState [in] The current state for the random number generator.
1122 ***********************************************************************************************************/
1123template <typename RNG>
1124LUPI_HOST_DEVICE void MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::sample( double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
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}
1196
1197/* *********************************************************************************************************//**
1198 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
1199 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
1200 *
1201 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
1202 * @param a_temperature [in] The temperature of the material.
1203 * @param a_energy_in [in] The energy of the incident particle.
1204 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
1205 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1206 * @param a_energy_out [in] The energy of the emitted outgoing particle.
1207 *
1208 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
1209 ***********************************************************************************************************/
1210
1211template <typename RNG>
1212LUPI_HOST_DEVICE double MCGIDI::Distributions::IncoherentPhotoAtomicScattering::angleBiasing( Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
1213 double a_energy_in, double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
1214
1215 URR_protareInfos URR_protareInfos1;
1216 double sigma = a_reaction->protareSingle( )->reactionCrossSection( a_reaction->reactionIndex( ), URR_protareInfos1, 0.0, a_energy_in );
1217
1218 double norm = M_PI * MCGIDI_classicalElectronRadius * MCGIDI_classicalElectronRadius / sigma;
1219
1220 double one_minus_mu = 1.0 - a_mu_lab;
1221 double k_in = a_energy_in / PoPI_electronMass_MeV_c2;
1222 a_energy_out = a_energy_in / ( 1.0 + k_in * one_minus_mu );
1223 double k_out = a_energy_out / PoPI_electronMass_MeV_c2;
1224
1225 double k_ratio = k_out / k_in;
1226 double probability = evaluateScatteringFactor( a_energy_in * sqrt( 0.5 * one_minus_mu ) );
1227 probability *= k_ratio * k_ratio * ( 1.0 + a_mu_lab * a_mu_lab + k_in * k_out * one_minus_mu * one_minus_mu ) * norm;
1228
1229 return( probability );
1230}
1231
1232/* *********************************************************************************************************//**
1233 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
1234 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
1235 *
1236 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
1237 * @param a_temperature [in] The temperature of the material.
1238 * @param a_energy_in [in] The energy of the incident particle.
1239 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
1240 * @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).
1241 * @param a_rngState [in] The current state for the random number generator.
1242 * @param a_energy_out [in] The energy of the emitted outgoing particle.
1243 *
1244 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
1245 ***********************************************************************************************************/
1246
1247template <typename RNG>
1248LUPI_HOST_DEVICE double MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::angleBiasing( Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
1249 double a_energy_in, double a_mu_lab, RNG && a_rng, double &a_energy_out ) const {
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}
1317
1318/* *********************************************************************************************************//**
1319 * This method returns the outgoing electron energy and angle given that the photon when out at an angle of *a_input.m_mu*.
1320 * Ergo, this method must be called directly after the photon has been sampled.
1321 *
1322 * @param a_energy [in] The energy of the projectile.
1323 * @param a_input [in] Sample options requested by user.
1324 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1325 ***********************************************************************************************************/
1326
1327template <typename RNG>
1328LUPI_HOST_DEVICE void MCGIDI::Distributions::IncoherentPhotoAtomicScatteringElectron::sample( double a_energy, Sampling::Input &a_input, RNG && a_rng ) const {
1329
1330 double halfTheta = 0.5 * acos( a_input.m_mu );
1331 double cot_psi = ( 1.0 + a_energy / PoPI_electronMass_MeV_c2 ) * tan( halfTheta );
1332 double psi = atan( 1.0 / cot_psi );
1333
1334 double deltaE_photon = a_energy - a_input.m_energyOut1;
1335 double electronMomentum2 = deltaE_photon * ( deltaE_photon + 2.0 * PoPI_electronMass_MeV_c2 ); // Square of the electron outlgoing momentum.
1336
1337 a_input.m_energyOut1 = electronMomentum2 / ( sqrt( electronMomentum2 + PoPI_electronMass_MeV_c2 * PoPI_electronMass_MeV_c2 ) + PoPI_electronMass_MeV_c2 );
1338 a_input.m_mu = cos( psi );
1339 a_input.m_phi = 2.0 * M_PI * a_rng( );
1340 a_input.m_frame = productFrame( );
1341}
1342
1343/* *********************************************************************************************************//**
1344 * Returns the probability for a projectile with energy *a_energy_in* causing a particle to be emitted
1345 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
1346 * Currently, this method only returns 0.0 for the probability and outgoing energy.
1347 *
1348 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
1349 * @param a_temperature [in] The temperature of the material.
1350 * @param a_energy_in [in] The energy of the incident particle.
1351 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
1352 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1353 * @param a_energy_out [in] The energy of the emitted outgoing particle.
1354 *
1355 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
1356 ***********************************************************************************************************/
1357
1358template <typename RNG>
1359LUPI_HOST_DEVICE double MCGIDI::Distributions::IncoherentPhotoAtomicScatteringElectron::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
1360 LUPI_maybeUnused double a_energy_in, LUPI_maybeUnused double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
1361
1362 a_energy_out = 0;
1363 return( 0.0 );
1364}
1365
1366/* *********************************************************************************************************//**
1367 * This method samples the outgoing photon by assigning the electron rest mass energy as the photon's energy and,
1368 * if m_firstSampled is true, randomly picking mu and phi. If m_firstSampled is false, the previous sampled particle
1369 * that filled in a_input must be the other sampled photon, then, the mu and phi for the second-sampled photon is such that
1370 * it is back-to-back with the other photon.
1371 *
1372 * @param a_X [in] The energy of the projectile.
1373 * @param a_input [in] Sample options requested by user.
1374 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1375 ***********************************************************************************************************/
1376
1377template <typename RNG>
1378LUPI_HOST_DEVICE void MCGIDI::Distributions::PairProductionGamma::sample( LUPI_maybeUnused double a_X, Sampling::Input &a_input, RNG && a_rng ) const {
1379
1380 if( m_firstSampled ) {
1381 a_input.m_mu = 1.0 - 2.0 * a_rng( );
1382 a_input.m_phi = M_PI * a_rng( ); }
1383 else {
1384 a_input.m_mu *= -1.0;
1385 a_input.m_phi += M_PI;
1386 }
1387 a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
1388 a_input.m_energyOut1 = PoPI_electronMass_MeV_c2;
1389 a_input.m_frame = productFrame( );
1390}
1391
1392/* *********************************************************************************************************//**
1393 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
1394 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
1395 *
1396 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
1397 * @param a_temperature [in] The temperature of the material.
1398 * @param a_energy_in [in] The energy of the incident particle.
1399 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
1400 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1401 * @param a_energy_out [in] The energy of the emitted outgoing particle.
1402 *
1403 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
1404 ***********************************************************************************************************/
1405
1406template <typename RNG>
1407LUPI_HOST_DEVICE double MCGIDI::Distributions::PairProductionGamma::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
1408 LUPI_maybeUnused double a_energy_in, LUPI_maybeUnused double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
1409
1410 a_energy_out = PoPI_electronMass_MeV_c2;
1411 return( 1.0 ); // 1.0 as there are two photons, each with 1/2 probability.
1412}
1413
1414/* *********************************************************************************************************//**
1415 * This method samples the outgoing neutron data for coherent elastic TSNL from the Debye/Waller function.
1416 *
1417 * @param a_energy [in] The energy of the projectile.
1418 * @param a_input [in] Sample options requested by user.
1419 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1420 ***********************************************************************************************************/
1421
1422template <typename RNG>
1423LUPI_HOST_DEVICE void MCGIDI::Distributions::CoherentElasticTNSL::sample( double a_energy, Sampling::Input &a_input,
1424 RNG && a_rng ) const {
1425
1426 if( a_energy <= m_energies[0] ) {
1427 a_input.m_mu = 1.0; }
1428 else {
1429 double temperature = a_input.temperature( );
1430 if( temperature < m_temperatures[0] ) temperature = m_temperatures[0];
1431 if( temperature > m_temperatures.back( ) ) temperature = m_temperatures.back( );
1432 std::size_t temperatureIndex = (std::size_t) MCGIDI::binarySearchVector( temperature, m_temperatures, true );
1433 double const *pointer1 = &m_S_table[temperatureIndex * m_energies.size( )];
1434
1435 double const *pointer2 = pointer1;
1436 double fractionFirstTemperature = 1.0;
1437 if( temperatureIndex != ( m_temperatures.size( ) - 1 ) ) {
1438 fractionFirstTemperature = ( m_temperatures[temperatureIndex+1] - temperature ) / ( m_temperatures[temperatureIndex+1] - m_temperatures[temperatureIndex] );
1439 pointer2 += m_energies.size( );
1440 }
1441 double fractionSecondTemperature = 1.0 - fractionFirstTemperature;
1442
1443 int intEnergyIndexMax = MCGIDI::binarySearchVector( a_energy, m_energies, true );
1444 std::size_t energyIndexMax = static_cast<std::size_t>( intEnergyIndexMax );
1445 if( a_energy == m_energies[energyIndexMax] ) --energyIndexMax;
1446
1447 double randomTotal = a_rng( ) * ( fractionFirstTemperature * pointer1[energyIndexMax] + fractionSecondTemperature * pointer2[energyIndexMax] );
1448 std::size_t energyIndex = 0;
1449 for( ; energyIndex < energyIndexMax; ++energyIndex ) {
1450 if( randomTotal <= fractionFirstTemperature * pointer1[energyIndex] + fractionSecondTemperature * pointer2[energyIndex] ) break;
1451 }
1452 a_input.m_mu = 1.0 - 2.0 * m_energies[energyIndex] / a_energy;
1453 }
1454
1455 a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
1456 a_input.m_energyOut1 = a_energy;
1457 a_input.m_phi = 2.0 * M_PI * a_rng( );
1458}
1459
1460/* *********************************************************************************************************//**
1461 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
1462 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
1463 *
1464 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
1465 * @param a_temperature [in] The temperature of the material.
1466 * @param a_energy_in [in] The energy of the incident particle.
1467 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
1468 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1469 * @param a_energy_out [in] The energy of the emitted outgoing particle.
1470 *
1471 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
1472 ***********************************************************************************************************/
1473
1474template <typename RNG>
1475LUPI_HOST_DEVICE double MCGIDI::Distributions::CoherentElasticTNSL::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
1476 double a_energy_in, LUPI_maybeUnused double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
1477
1478 double probability = 0.0;
1479
1480 a_energy_out = a_energy_in;
1481
1482 return( probability );
1483}
1484
1485/* *********************************************************************************************************//**
1486 * This method samples the outgoing neutron data for incoherent elastic TSNL from the Debye/Waller function.
1487 *
1488 * @param a_energy [in] The energy of the projectile.
1489 * @param a_input [in] Sample options requested by user.
1490 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1491 ***********************************************************************************************************/
1492
1493template <typename RNG>
1494LUPI_HOST_DEVICE void MCGIDI::Distributions::IncoherentElasticTNSL::sample( double a_energy, Sampling::Input &a_input,
1495 RNG && a_rng ) const {
1496
1497 double temperature = a_input.temperature( ) / m_temperatureToMeV_K;
1498 double W_prime = m_DebyeWallerIntegral->evaluate( temperature );
1499 double twoEW = 2 * a_energy * W_prime;
1500 double expOfTwice_twoEW = exp( -2 * twoEW );
1501 double sampled_cdf = a_rng( );
1502
1503 if( sampled_cdf > ( 1 - 1e-5 ) ) {
1504 double Variable = ( 1.0 - sampled_cdf ) * ( 1.0 - expOfTwice_twoEW );
1505 a_input.m_mu = 1.0 - Variable * ( 1.0 + 0.5 * Variable ) / twoEW; }
1506 else if( sampled_cdf < expOfTwice_twoEW ) {
1507 a_input.m_mu = -1.0 + log( sampled_cdf / expOfTwice_twoEW * ( 1.0 - expOfTwice_twoEW ) + 1.0 ) / twoEW; }
1508 else {
1509 a_input.m_mu = 1.0 + log( expOfTwice_twoEW + sampled_cdf * ( 1.0 - expOfTwice_twoEW ) ) / twoEW;
1510 }
1511
1512 a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
1513 a_input.m_energyOut1 = a_energy;
1514 a_input.m_phi = 2.0 * M_PI * a_rng( );
1515}
1516
1517/* *********************************************************************************************************//**
1518 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
1519 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy.
1520 *
1521 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
1522 * @param a_temperature [in] The temperature of the material.
1523 * @param a_energy_in [in] The energy of the incident particle.
1524 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
1525 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1526 * @param a_energy_out [in] The energy of the emitted outgoing particle.
1527 *
1528 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
1529 ***********************************************************************************************************/
1530
1531template <typename RNG>
1532LUPI_HOST_DEVICE double MCGIDI::Distributions::IncoherentElasticTNSL::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, double a_temperature,
1533 double a_energy_in, double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
1534
1535 double temperature = a_temperature / m_temperatureToMeV_K;
1536 double W_prime = m_DebyeWallerIntegral->evaluate( temperature );
1537 double twoEW = 2 * a_energy_in * W_prime;
1538 double probability = exp( -twoEW * ( 1.0 - a_mu_lab ) ) * twoEW / ( 1.0 - exp( -2 * twoEW ) );
1539
1540 a_energy_out = a_energy_in;
1541
1542 return( probability );
1543}
1544
1545/* *********************************************************************************************************//**
1546 * The method sets all outgoing product data to 0.0 and set the sampledType to Sampling::unspecified.
1547 *
1548 * @param a_X [in] The energy of the projectile.
1549 * @param a_input [in] Sample options requested by user.
1550 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1551 ***********************************************************************************************************/
1552
1553template <typename RNG>
1554LUPI_HOST_DEVICE void MCGIDI::Distributions::Unspecified::sample( LUPI_maybeUnused double a_X, Sampling::Input &a_input, LUPI_maybeUnused RNG && a_rng ) const {
1555
1556 a_input.setSampledType( Sampling::SampledType::unspecified );
1557 a_input.m_energyOut1 = 0.;
1558 a_input.m_mu = 0.;
1559 a_input.m_phi = 0.;
1560 a_input.m_frame = productFrame( );
1561}
1562
1563/* *********************************************************************************************************//**
1564 * Returns the probability for a projectile with energy *a_energy_in* to cause a particle to be emitted
1565 * at angle *a_mu_lab* as seen in the lab frame. *a_energy_out* is the sampled outgoing energy. This one should never
1566 * be called. If called, returns 0.0 for a probability.
1567 *
1568 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
1569 * @param a_temperature [in] Specifies the temperature of the material.
1570 * @param a_energy_in [in] The energy of the incident particle.
1571 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
1572 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1573 * @param a_energy_out [in] The energy of the emitted outgoing particle.
1574 *
1575 * @return The probability of emitting outgoing particle into lab angle *a_mu_lab*.
1576 ***********************************************************************************************************/
1577
1578template <typename RNG>
1579LUPI_HOST_DEVICE double MCGIDI::Distributions::Unspecified::angleBiasing( LUPI_maybeUnused Reaction const *a_reaction, LUPI_maybeUnused double a_temperature,
1580 LUPI_maybeUnused double a_energy_in, LUPI_maybeUnused double a_mu_lab, LUPI_maybeUnused RNG && a_rng, double &a_energy_out ) const {
1581
1582 a_energy_out = 0.0;
1583 return( 0.0 );
1584}
1585
1586
1587// From file: MCGIDI_functions.cpp
1588
1589/*
1590============================================================
1591*/
1592template <typename RNG >
1593LUPI_HOST_DEVICE int MCGIDI::Functions::Function1d::sampleBoundingInteger( double a_x1, RNG && a_rng ) const {
1594
1595 if( type( ) == Function1dType::TerrellFissionNeutronMultiplicityModel )
1596 return( static_cast<TerrellFissionNeutronMultiplicityModel const *>( this )->sampleBoundingInteger( a_x1, a_rng ) );
1597
1598 double d_value = evaluate( a_x1 );
1599 int iValue = (int) d_value;
1600 if( iValue == d_value ) return( iValue );
1601 if( d_value - iValue > a_rng( ) ) ++iValue;
1602
1603 return( iValue );
1604}
1605
1606/* *********************************************************************************************************//**
1607 * Sample the number of fission prompt neutrons using Terrell's modified Gaussian distribution.
1608 * Method uses Red Cullen's algoritm (see UCRL-TR-222526).
1609 *
1610 * @param a_energy [in] The energy of the projectile.
1611 * @param a_rng [in] The random number generator function the uses *a_rngState* to generator a double in the range [0, 1.0).
1612 * @param a_rngState [in/out] The random number generator state.
1613 *
1614 * @return The sampled number of emitted, prompt neutrons for fission.
1615 ***********************************************************************************************************/
1616template <typename RNG>
1617LUPI_HOST_DEVICE int MCGIDI::Functions::TerrellFissionNeutronMultiplicityModel::sampleBoundingInteger( double a_energy, RNG && a_rng ) const {
1618
1619 const double Terrell_BSHIFT = -0.43287;
1620 double width = M_SQRT2 * m_width;
1621 double temp1 = m_multiplicity->evaluate( a_energy ) + 0.5;
1622 double temp2 = temp1 / width;
1623 double expo = exp( -temp2 * temp2 );
1624 double cshift = temp1 + Terrell_BSHIFT * m_width * expo / ( 1.0 - expo );
1625
1626 double multiplicity = 1.0;
1627 do {
1628 double rw = sqrt( -log( a_rng( ) ) );
1629 double theta = ( 2.0 * M_PI ) * a_rng();
1630
1631 multiplicity = width * rw * cos( theta ) + cshift;
1632 } while ( multiplicity < 0.0 );
1633
1634 return( static_cast<int>( floor( multiplicity ) ) );
1635}
1636
1637/* *********************************************************************************************************//**
1638 * Returns the x-value corresponding cumulative probability *a_rngValue*.
1639 *
1640 * @param a_rngValue [in] The x-value to evaluate the function at.
1641 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1642 *
1643 * @return The value of the function at *a_x1*.
1644 ***********************************************************************************************************/
1645
1646template <typename RNG>
1647LUPI_HOST_DEVICE double MCGIDI::Probabilities::ProbabilityBase1d::sample( double a_rngValue, RNG && a_rng ) const {
1648
1649 return( static_cast<Xs_pdf_cdf1d const *>( this )->sample( a_rngValue, a_rng ) );
1650}
1651
1652template <typename RNG>
1653LUPI_HOST_DEVICE double MCGIDI::Probabilities::Xs_pdf_cdf1d::sample( double a_rngValue, LUPI_maybeUnused RNG && a_rng ) const {
1654
1655 int intLower = binarySearchVector( a_rngValue, m_cdf );
1656 double domainValue = 0;
1657
1658 if( intLower < 0 ) { // This should never happen.
1659 LUPI_THROW( "Xs_pdf_cdf1d::sample: intLower < 0." );
1660 }
1661
1662 std::size_t lower = static_cast<std::size_t>( intLower );
1663
1664 if( interpolation( ) == Interpolation::FLAT ) {
1665 double fraction = ( m_cdf[lower+1] - a_rngValue ) / ( m_cdf[lower+1] - m_cdf[lower] );
1666 domainValue = fraction * m_Xs[lower] + ( 1 - fraction ) * m_Xs[lower+1]; }
1667 else { // Assumes lin-lin interpolation.
1668 double slope = m_pdf[lower+1] - m_pdf[lower];
1669
1670 if( slope == 0.0 ) {
1671 if( m_pdf[lower] == 0.0 ) {
1672 domainValue = m_Xs[lower];
1673 if( lower == 0 ) domainValue = m_Xs[1]; }
1674 else {
1675 double fraction = ( m_cdf[lower+1] - a_rngValue ) / ( m_cdf[lower+1] - m_cdf[lower] );
1676 domainValue = fraction * m_Xs[lower] + ( 1 - fraction ) * m_Xs[lower+1];
1677 } }
1678 else {
1679 double d1, d2;
1680
1681 slope = slope / ( m_Xs[lower+1] - m_Xs[lower] );
1682 d1 = a_rngValue - m_cdf[lower];
1683 d2 = m_cdf[lower+1] - a_rngValue;
1684 if( d2 > d1 ) { // Closer to lower.
1685 domainValue = m_Xs[lower] + ( sqrt( m_pdf[lower] * m_pdf[lower] + 2. * slope * d1 ) - m_pdf[lower] ) / slope; }
1686 else { // Closer to lower + 1.
1687 domainValue = m_Xs[lower+1] - ( m_pdf[lower+1] - sqrt( m_pdf[lower+1] * m_pdf[lower+1] - 2. * slope * d2 ) ) / slope;
1688 }
1689 }
1690 }
1691 return( domainValue );
1692}
1693
1694/* *********************************************************************************************************//**
1695 * This method samples an x1 from a pdf(x1|x2) given x2 and the cumulative value of the pdf as *a_rngValue*.
1696 *
1697 * @param a_x2 [in] The value of x2.
1698 * @param a_rngValue [in] The value of the cumulative used to determine the x1 value.
1699 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1700 ***********************************************************************************************************/
1701
1702template <typename RNG>
1703LUPI_HOST_DEVICE double MCGIDI::Probabilities::ProbabilityBase2d::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
1704
1705 double value = 0.0;
1706
1707 switch( type( ) ) {
1708 case ProbabilityBase2dType::none:
1709 break;
1710 case ProbabilityBase2dType::weightedFunctionals:
1711 value = static_cast<WeightedFunctionals2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
1712 break;
1713 default:
1714 value = static_cast<ProbabilityBase2d_d1 const *>( this )->sample( a_x2, a_rngValue, a_rng );
1715 break;
1716 }
1717
1718 return( value );
1719}
1720
1721/* *********************************************************************************************************//**
1722 * Returns the value of x1, given x2 and the cumulative probability *a_rngValue*.
1723 *
1724 * @param a_x2 [in] Value of the outer most independent variable (i.e., *x2*).
1725 * @param a_rngValue [in] The value of the cumulative probability used to determine the x1 value.
1726 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1727 *
1728 * @return The value *x1* where the cumulative probability is *a_rngValue* for x2 = *a_x2*.
1729 ***********************************************************************************************************/
1730
1731template <typename RNG>
1732LUPI_HOST_DEVICE double MCGIDI::Probabilities::ProbabilityBase2d_d1::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
1733
1734 double value = 0.0;
1735
1736 switch( type( ) ) {
1737 case ProbabilityBase2dType::XYs:
1738 case ProbabilityBase2dType::isotropic:
1739 case ProbabilityBase2dType::discreteGamma:
1740 case ProbabilityBase2dType::primaryGamma:
1741 case ProbabilityBase2dType::recoil:
1742 case ProbabilityBase2dType::NBodyPhaseSpace:
1743 case ProbabilityBase2dType::evaporation:
1744 case ProbabilityBase2dType::generalEvaporation:
1745 case ProbabilityBase2dType::simpleMaxwellianFission:
1746 case ProbabilityBase2dType::Watt:
1747 value = static_cast<ProbabilityBase2d_d2 const *>( this )->sample( a_x2, a_rngValue, a_rng );
1748 break;
1749 case ProbabilityBase2dType::regions:
1750 value = static_cast<Regions2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
1751 break;
1752 case ProbabilityBase2dType::none:
1753 case ProbabilityBase2dType::weightedFunctionals:
1754 LUPI_THROW( "ProbabilityBase2d_d1::sample: This should never happen." );
1755 }
1756
1757 return( value );
1758}
1759
1760/* *********************************************************************************************************//**
1761 * This method returns two x1 values for use with ProbabilityBase3d functions.
1762 *
1763 * @param a_x2 [in] The value of x2.
1764 * @param a_rngValue [in] The value of the cumulative value used to determine the x1 value.
1765 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1766 * @param a_x1_1 [in] The lower value of the x1 value.
1767 * @param a_x1_2 [in] The upper value of the x1 value.
1768 ***********************************************************************************************************/
1769
1770template <typename RNG>
1771LUPI_HOST_DEVICE double MCGIDI::Probabilities::ProbabilityBase2d_d1::sample2dOf3d( double a_x2, double a_rngValue, RNG && a_rng,
1772 double *a_x1_1, double *a_x1_2 ) const {
1773
1774 double value = 0.0;
1775
1776 switch( type( ) ) {
1777 case ProbabilityBase2dType::XYs:
1778 value = static_cast<XYs2d const *>( this )->sample2dOf3d( a_x2, a_rngValue, a_rng, a_x1_1, a_x1_2 );
1779 break;
1780 default:
1781 LUPI_THROW( "ProbabilityBase2d_d1::sample2dOf3d: not implemented." );
1782 }
1783
1784 return( value );
1785}
1786
1787/* *********************************************************************************************************//**
1788 * Returns the value of x1, given x2 and the cumulative probability *a_rngValue*.
1789 *
1790 * @param a_x2 [in] Value of the outer most independent variable (i.e., *x2*).
1791 * @param a_rngValue [in] The value of the cumulative probability used to determine the x1 value.
1792 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1793 *
1794 * @return The value *x1* where the cumulative probability is *a_rngValue* for x2 = *a_x2*.
1795 ***********************************************************************************************************/
1796
1797template <typename RNG>
1798LUPI_HOST_DEVICE double MCGIDI::Probabilities::ProbabilityBase2d_d2::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
1799
1800 double value = 0.0;
1801
1802 switch( type( ) ) {
1803 case ProbabilityBase2dType::XYs:
1804 value = static_cast<XYs2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
1805 break;
1806 case ProbabilityBase2dType::isotropic:
1807 value = static_cast<Isotropic2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
1808 break;
1809 case ProbabilityBase2dType::discreteGamma:
1810 value = static_cast<DiscreteGamma2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
1811 break;
1812 case ProbabilityBase2dType::primaryGamma:
1813 value = static_cast<PrimaryGamma2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
1814 break;
1815 case ProbabilityBase2dType::recoil:
1816 value = static_cast<Recoil2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
1817 break;
1818 case ProbabilityBase2dType::NBodyPhaseSpace:
1819 value = static_cast<NBodyPhaseSpace2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
1820 break;
1821 case ProbabilityBase2dType::evaporation:
1822 value = static_cast<Evaporation2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
1823 break;
1824 case ProbabilityBase2dType::generalEvaporation:
1825 value = static_cast<GeneralEvaporation2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
1826 break;
1827 case ProbabilityBase2dType::simpleMaxwellianFission:
1828 value = static_cast<SimpleMaxwellianFission2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
1829 break;
1830 case ProbabilityBase2dType::Watt:
1831 value = static_cast<Watt2d const *>( this )->sample( a_x2, a_rngValue, a_rng );
1832 break;
1833 case ProbabilityBase2dType::none:
1834 case ProbabilityBase2dType::weightedFunctionals:
1835 case ProbabilityBase2dType::regions:
1836 LUPI_THROW( "ProbabilityBase2d_d2::sample: This should never happen." );
1837 }
1838
1839 return( value );
1840}
1841
1842/* *********************************************************************************************************//**
1843 * This method returns two x1 values for use with ProbabilityBase3d functions.
1844 *
1845 * @param a_x2 [in] The value of x2.
1846 * @param a_rngValue [in] The value of the cumulative value used to determine the x1 value.
1847 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1848 * @param a_x1_1 [in] The lower value of the x1 value.
1849 * @param a_x1_2 [in] The upper value of the x1 value.
1850 ***********************************************************************************************************/
1851
1852template <typename RNG>
1853LUPI_HOST_DEVICE double MCGIDI::Probabilities::ProbabilityBase2d_d2::sample2dOf3d( double a_x2, double a_rngValue, RNG && a_rng,
1854 double *a_x1_1, double *a_x1_2 ) const {
1855
1856 double value = 0.0;
1857
1858 switch( type( ) ) {
1859 case ProbabilityBase2dType::XYs:
1860 value = static_cast<XYs2d const *>( this )->sample2dOf3d( a_x2, a_rngValue, a_rng, a_x1_1, a_x1_2 );
1861 break;
1862 default:
1863 LUPI_THROW( "ProbabilityBase2d_d2::sample2dOf3d: not implemented." );
1864 }
1865
1866 return( value );
1867}
1868
1869template <typename RNG>
1870LUPI_HOST_DEVICE double MCGIDI::Probabilities::XYs2d::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
1871/*
1872C Samples from a pdf(x1|x2). First determine which pdf(s) to sample from given x2.
1873C Then use rngValue to sample from pdf1(x1) and maybe pdf2(x1) and interpolate to
1874C determine x1.
1875*/
1876 double sampledValue = 0;
1877 int intLower = binarySearchVector( a_x2, m_Xs );
1878
1879 if( intLower == -2 ) {
1880 sampledValue = m_probabilities[0]->sample( a_rngValue, a_rng ); }
1881 else if( intLower == -1 ) {
1882 sampledValue = m_probabilities.back( )->sample( a_rngValue, a_rng ); }
1883 else {
1884 auto lower = static_cast<std::size_t>( intLower );
1885 double sampled1 = m_probabilities[lower]->sample( a_rngValue, a_rng );
1886
1887 if( interpolation( ) == Interpolation::FLAT ) {
1888 sampledValue = sampled1; }
1889 else {
1890 double sampled2 = m_probabilities[lower+1]->sample( a_rngValue, a_rng );
1891
1892 if( interpolation( ) == Interpolation::LINLIN ) {
1893 double fraction = ( m_Xs[lower+1] - a_x2 ) / ( m_Xs[lower+1] - m_Xs[lower] );
1894 sampledValue = fraction * sampled1 + ( 1 - fraction ) * sampled2; }
1895 else if( interpolation( ) == Interpolation::LOGLIN ) {
1896 double fraction = ( m_Xs[lower+1] - a_x2 ) / ( m_Xs[lower+1] - m_Xs[lower] );
1897 sampledValue = sampled2 * pow( sampled2 / sampled1, fraction ); }
1898 else if( interpolation( ) == Interpolation::LINLOG ) {
1899 double fraction = log( m_Xs[lower+1] / a_x2 ) / log( m_Xs[lower+1] / m_Xs[lower] );
1900 sampledValue = fraction * sampled1 + ( 1 - fraction ) * sampled2; }
1901 else if( interpolation( ) == Interpolation::LOGLOG ) {
1902 double fraction = log( m_Xs[lower+1] / a_x2 ) / log( m_Xs[lower+1] / m_Xs[lower] );
1903 sampledValue = sampled2 * pow( sampled2 / sampled1, fraction ); }
1904 else { // This should never happen.
1905 LUPI_THROW( "XYs2d::sample: unsupported interpolation." );
1906 }
1907 }
1908 }
1909 return( sampledValue );
1910}
1911
1912/* *********************************************************************************************************//**
1913 * This method returns two x1 values for use with ProbabilityBase3d functions.
1914 *
1915 * @param a_x2 [in] The value of x2.
1916 * @param a_rngValue [in] The value of the cumulative value used to determine the x1 value.
1917 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
1918 * @param a_x1_1 [in] The lower value of the x1 value.
1919 * @param a_x1_2 [in] The upper value of the x1 value.
1920 ***********************************************************************************************************/
1921
1922template <typename RNG>
1923LUPI_HOST_DEVICE double MCGIDI::Probabilities::XYs2d::sample2dOf3d( double a_x2, double a_rngValue, RNG && a_rng,
1924 double *a_x1_1, double *a_x1_2 ) const {
1925/*
1926C Samples from a pdf(x1|x2). First determine which pdf(s) to sample from given x2. Then use rngValue to sample from pdf1(x1)
1927C and maybe pdf2(x1) and interpolate to determine x1.
1928*/
1929 double sampledValue = 0;
1930 int intLower = binarySearchVector( a_x2, m_Xs );
1931
1932 if( intLower == -2 ) {
1933 sampledValue = m_probabilities[0]->sample( a_rngValue, a_rng );
1934 *a_x1_2 = *a_x1_1 = sampledValue; }
1935 else if( intLower == -1 ) {
1936 sampledValue = m_probabilities.back( )->sample( a_rngValue, a_rng );
1937 *a_x1_2 = *a_x1_1 = sampledValue; }
1938 else {
1939 std::size_t lower = static_cast<std::size_t>( intLower );
1940 *a_x1_1 = m_probabilities[lower]->sample( a_rngValue, a_rng );
1941
1942 if( interpolation( ) == Interpolation::FLAT ) {
1943 sampledValue = *a_x1_2 = *a_x1_1; }
1944 else {
1945 *a_x1_2 = m_probabilities[lower+1]->sample( a_rngValue, a_rng );
1946
1947 if( interpolation( ) == Interpolation::LINLIN ) {
1948 double fraction = ( m_Xs[lower+1] - a_x2 ) / ( m_Xs[lower+1] - m_Xs[lower] );
1949 sampledValue = fraction * *a_x1_1 + ( 1 - fraction ) * *a_x1_2; }
1950 else if( interpolation( ) == Interpolation::LOGLIN ) {
1951 double fraction = ( m_Xs[lower+1] - a_x2 ) / ( m_Xs[lower+1] - m_Xs[lower] );
1952 sampledValue = *a_x1_2 * pow( *a_x1_2 / *a_x1_1, fraction ); }
1953 else if( interpolation( ) == Interpolation::LINLOG ) {
1954 double fraction = log( m_Xs[lower+1] / a_x2 ) / log( m_Xs[lower+1] / m_Xs[lower] );
1955 sampledValue = fraction * *a_x1_1 + ( 1 - fraction ) * *a_x1_2; }
1956 else if( interpolation( ) == Interpolation::LOGLOG ) {
1957 double fraction = log( m_Xs[lower+1] / a_x2 ) / log( m_Xs[lower+1] / m_Xs[lower] );
1958 sampledValue = *a_x1_2 * pow( *a_x1_2 / *a_x1_1 , fraction ); }
1959 else { // This should never happen.
1960 LUPI_THROW( "XYs2d::sample: unsupported interpolation." );
1961 }
1962 }
1963 }
1964 return( sampledValue );
1965}
1966
1967template <typename RNG>
1968LUPI_HOST_DEVICE double MCGIDI::Probabilities::Regions2d::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
1969
1970 int intLower = binarySearchVector( a_x2, m_Xs );
1971
1972 if( intLower < 0 ) {
1973 if( intLower == -1 ) { // a_x2 > last value of m_Xs.
1974 return( m_probabilities.back( )->sample( a_x2, a_rngValue, a_rng ) );
1975 }
1976 intLower = 0; // a_x2 < first value of m_Xs.
1977 }
1978
1979 std::size_t lower = static_cast<std::size_t>( intLower );
1980
1981 return( m_probabilities[lower]->sample( a_x2, a_rngValue, a_rng ) );
1982}
1983
1984template <typename RNG>
1985LUPI_HOST_DEVICE double MCGIDI::Probabilities::Recoil2d::sample( LUPI_maybeUnused double a_x2, LUPI_maybeUnused double a_rngValue, LUPI_maybeUnused RNG && a_rng ) const {
1986
1987#if !defined(__NVCC__) && !defined(__HIP__)
1988 LUPI_THROW( "Recoil2d::sample: not implemented." );
1989#endif
1990
1991 return( 0.0 );
1992}
1993
1994/* *********************************************************************************************************
1995 * Sampling for NBody phase space.
1996 *
1997 * @param a_x2 [in] Incident energy of the projectile.
1998 * @param a_rngValue [in] The GIDI::Protare whose data is to be used to construct *this*.
1999 * @param a_rng [in] This argument is not used by this method but needed to match ProbabilityBase1d::sample's definition.
2000 ***********************************************************************************************************/
2001
2002template <typename RNG>
2003LUPI_HOST_DEVICE double MCGIDI::Probabilities::NBodyPhaseSpace2d::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
2004
2005 double energyMax = m_energy_in_COMFactor * a_x2 + m_Q;
2006
2007 if( energyMax < 0.0 ) return( 0.0 ); // Kludge for now until for upscatter model A until sampling below threshold is fixed.
2008
2009 return( energyMax * m_massFactor * m_dist->sample( a_rngValue, a_rng ) );
2010}
2011
2012inline LUPI_HOST_DEVICE static double MCGIDI_sampleEvaporation( double a_xMax, double a_rngValue ) {
2013
2014 double b1, c1, xMid, norm, xMin = 0.;
2015
2016 norm = 1 - ( 1 + a_xMax ) * exp( -a_xMax );
2017 b1 = 1. - norm * a_rngValue;
2018 for( int i1 = 0; i1 < 16; i1++ ) {
2019 xMid = 0.5 * ( xMin + a_xMax );
2020 c1 = ( 1 + xMid ) * exp( -xMid );
2021 if( b1 > c1 ) {
2022 a_xMax = xMid; }
2023 else {
2024 xMin = xMid;
2025 }
2026 }
2027 return( xMid );
2028}
2029
2030
2031template <typename RNG>
2032LUPI_HOST_DEVICE double MCGIDI::Probabilities::Evaporation2d::sample( double a_x2, double a_rngValue, LUPI_maybeUnused RNG && a_rng ) const {
2033
2034 double theta = m_theta->evaluate( a_x2 );
2035
2036 return( theta * MCGIDI_sampleEvaporation( ( a_x2 - m_U ) / theta, a_rngValue ) );
2037}
2038
2039template <typename RNG>
2040LUPI_HOST_DEVICE double MCGIDI::Probabilities::GeneralEvaporation2d::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
2041
2042 return( m_theta->evaluate( a_x2 ) * m_g->sample( a_rngValue, a_rng ) );
2043}
2044
2045inline LUPI_HOST_DEVICE static double MCGIDI_sampleSimpleMaxwellianFission( double a_xMax, double a_rngValue ) {
2046
2047 double b1, c1, xMid, norm, xMin = 0., sqrt_xMid, sqrt_pi_2 = 0.5 * sqrt( M_PI );
2048
2049 sqrt_xMid = sqrt( a_xMax );
2050 norm = sqrt_pi_2 * erf( sqrt_xMid ) - sqrt_xMid * exp( -a_xMax );
2051 b1 = norm * a_rngValue;
2052 for( int i1 = 0; i1 < 16; i1++ ) {
2053 xMid = 0.5 * ( xMin + a_xMax );
2054 sqrt_xMid = sqrt( xMid );
2055 c1 = sqrt_pi_2 * erf( sqrt_xMid ) - sqrt_xMid * exp( -xMid );
2056 if( b1 < c1 ) {
2057 a_xMax = xMid; }
2058 else {
2059 xMin = xMid;
2060 }
2061 }
2062 return( xMid );
2063}
2064
2065template <typename RNG>
2066LUPI_HOST_DEVICE double MCGIDI::Probabilities::SimpleMaxwellianFission2d::sample( double a_x2, double a_rngValue, LUPI_maybeUnused RNG && a_rng ) const {
2067
2068 double theta = m_theta->evaluate( a_x2 );
2069
2070 return( theta * MCGIDI_sampleSimpleMaxwellianFission( ( a_x2 - m_U ) / theta, a_rngValue ) );
2071}
2072
2073template <typename RNG>
2074LUPI_HOST_DEVICE double MCGIDI::Probabilities::Watt2d::sample( double a_x2, LUPI_maybeUnused double a_rngValue, RNG && a_rng ) const {
2075/*
2076* From MCAPM via Sample Watt Spectrum as in TART ( Kalos algorithm ).
2077*/
2078 double WattMin = 0., WattMax = a_x2 - m_U, x, y, z, energyOut, rand1, rand2;
2079 double Watt_a = 1./m_a->evaluate( a_x2 ); // Kalos algorithm uses the inverse of the 'a' parameter stored in GNDS
2080 double Watt_b = m_b->evaluate( a_x2 );
2081
2082 x = 1. + ( Watt_b / ( 8. * Watt_a ) );
2083 y = ( x + sqrt( x * x - 1. ) ) / Watt_a;
2084 z = Watt_a * y - 1.;
2085 do {
2086 rand1 = -log( a_rng( ) );
2087 rand2 = -log( a_rng( ) );
2088 energyOut = y * rand1;
2089 } while( ( ( rand2 - z * ( rand1 + 1. ) ) * ( rand2 - z * ( rand1 + 1. ) ) > Watt_b * y * rand1 ) ||
2090 ( energyOut < WattMin ) || ( energyOut > WattMax ) );
2091 return( energyOut );
2092}
2093
2094
2095template <typename RNG>
2096LUPI_HOST_DEVICE double MCGIDI::Probabilities::WeightedFunctionals2d::sample( double a_x2, double a_rngValue, RNG && a_rng ) const {
2097/*
2098c This routine assumes that the weights sum to 1.
2099*/
2100 std::size_t i1;
2101 std::size_t n1 = m_weight.size( ) - 1; // Take last point if others do not add to randomWeight.
2102 double randomWeight = a_rng( ), cumulativeWeight = 0.;
2103
2104 for( i1 = 0; i1 < n1; ++i1 ) {
2105 cumulativeWeight += m_weight[i1]->evaluate( a_x2 );
2106 if( cumulativeWeight >= randomWeight ) break;
2107 }
2108 return( m_energy[i1]->sample( a_x2, a_rngValue, a_rng) );
2109}
2110
2111/* *********************************************************************************************************//**
2112 * This method samples an x1 from a pdf(x1|x2) given x2 and the cumulative value of the pdf as *a_rngValue*.
2113 *
2114 * @param a_x3 [in] The value of x3.
2115 * @param a_x2_1 [in] The value of ?.
2116 * @param a_x2_2 [in] The value of ?.
2117 * @param a_rngValue [in] The value of the cumulative used to determine the x1 value.
2118 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2119 ***********************************************************************************************************/
2120
2121template <typename RNG>
2122LUPI_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 {
2123
2124 return( static_cast<XYs3d const *>( this )->sample( a_x3, a_x2_1, a_x2_2, a_rngValue, a_rng ) );
2125}
2126
2127template <typename RNG>
2128LUPI_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 {
2129/*
2130C Samples from a pdf(x1|x3,x2). First determine which pdf(s) to sample from given x3
2131C Then use rngValue to sample from pdf2_1(x2) and maybe pdf2_2(x2) and interpolate to
2132C determine x1.
2133*/
2134 double sampledValue = 0;
2135 int intLower = binarySearchVector( a_x3, m_Xs );
2136
2137 if( intLower == -2 ) { // x3 < first value of Xs.
2138 sampledValue = m_probabilities[0]->sample( a_x2_1, a_rngValue, a_rng ); }
2139 else if( intLower == -1 ) { // x3 > last value of Xs.
2140 sampledValue = m_probabilities.back( )->sample( a_x2_1, a_rngValue, a_rng ); }
2141 else {
2142 std::size_t lower = static_cast<std::size_t>( intLower );
2143 double sampled1 = m_probabilities[lower]->sample( a_x2_1, a_rngValue, a_rng );
2144
2145 if( interpolation( ) == Interpolation::FLAT ) {
2146 sampledValue = sampled1; }
2147 else {
2148 double sampled2 = m_probabilities[lower+1]->sample( a_x2_2, a_rngValue, a_rng );
2149
2150 if( interpolation( ) == Interpolation::LINLIN ) {
2151 double fraction = ( m_Xs[lower+1] - a_x3 ) / ( m_Xs[lower+1] - m_Xs[lower] );
2152 sampledValue = fraction * sampled1 + ( 1 - fraction ) * sampled2; }
2153 else if( interpolation( ) == Interpolation::LOGLIN ) {
2154 double fraction = ( m_Xs[lower+1] - a_x3 ) / ( m_Xs[lower+1] - m_Xs[lower] );
2155 sampledValue = sampled2 * pow( sampled2 / sampled1, fraction ); }
2156 else if( interpolation( ) == Interpolation::LINLOG ) {
2157 double fraction = log( m_Xs[lower+1] / a_x3 ) / log( m_Xs[lower+1] / m_Xs[lower] );
2158 sampledValue = fraction * sampled1 + ( 1 - fraction ) * sampled2; }
2159 else if( interpolation( ) == Interpolation::LOGLOG ) {
2160 double fraction = log( m_Xs[lower+1] / a_x3 ) / log( m_Xs[lower+1] / m_Xs[lower] );
2161 sampledValue = sampled2 * pow( sampled2 / sampled1, fraction ); }
2162 else { // This should never happen.
2163 LUPI_THROW( "XYs3d::sample: unsupported interpolation." );
2164 }
2165 }
2166 }
2167
2168 return( sampledValue );
2169}
2170
2171
2172// From file: MCGIDI_heatedCrossSections.cpp
2173
2174/* *********************************************************************************************************//**
2175 * Returns the requested reaction's multi-group cross section for target temperature *a_temperature* and projectile multi-group *a_hashIndex*.
2176 *
2177 * @param a_URR_protareInfos [in] URR information.
2178 * @param a_URR_index [in] If not negative, specifies the index in *a_URR_protareInfos*.
2179 * @param a_hashIndex [in] Specifies projectile energy hash index.
2180 * @param a_temperature [in] The temperature of the target.
2181 * @param a_energy [in] The energy of the projectile.
2182 * @param a_crossSection [in] The total cross section for the protare at *a_temperature* and *a_energy*.
2183 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2184 ***********************************************************************************************************/
2185
2186template <typename RNG>
2187LUPI_HOST_DEVICE std::size_t MCGIDI::HeatedCrossSectionsContinuousEnergy::sampleReaction( URR_protareInfos const &a_URR_protareInfos,
2188 int a_URR_index, std::size_t a_hashIndex, double a_temperature, double a_energy, double a_crossSection, RNG && a_rng ) const {
2189
2190 std::size_t sampled_reaction_index, temperatureIndex1, temperatureIndex2, number_of_temperatures = m_temperatures.size( );
2191 double sampleCrossSection = a_crossSection * a_rng( );
2192
2193 if( a_temperature <= m_temperatures[0] ) {
2194 temperatureIndex1 = 0;
2195 temperatureIndex2 = temperatureIndex1; }
2196 else if( a_temperature >= m_temperatures.back( ) ) {
2197 temperatureIndex1 = m_temperatures.size( ) - 1;
2198 temperatureIndex2 = temperatureIndex1; }
2199 else {
2200 std::size_t i1 = 0;
2201 for( ; i1 < number_of_temperatures; ++i1 ) if( a_temperature < m_temperatures[i1] ) break;
2202 temperatureIndex1 = i1 - 1;
2203 temperatureIndex2 = i1;
2204 }
2205
2206 std::size_t numberOfReactions = m_heatedCrossSections[0]->numberOfReactions( );
2207 double energyFraction1, energyFraction2, crossSectionSum = 0.0;
2208
2209 HeatedCrossSectionContinuousEnergy &heatedCrossSection1 = *m_heatedCrossSections[temperatureIndex1];
2210 std::size_t energyIndex1 = heatedCrossSection1.evaluationInfo( a_hashIndex, a_energy, &energyFraction1 );
2211
2212 if( temperatureIndex1 == temperatureIndex2 ) {
2213 for( sampled_reaction_index = 0; sampled_reaction_index < numberOfReactions; ++sampled_reaction_index ) {
2214 crossSectionSum += heatedCrossSection1.reactionCrossSection2( sampled_reaction_index, a_URR_protareInfos, a_URR_index, a_energy,
2215 energyIndex1, energyFraction1 );
2216 if( crossSectionSum >= sampleCrossSection ) break;
2217 } }
2218 else {
2219 double temperatureFraction2 = ( a_temperature - m_temperatures[temperatureIndex1] )
2220 / ( m_temperatures[temperatureIndex2] - m_temperatures[temperatureIndex1] );
2221 double temperatureFraction1 = 1.0 - temperatureFraction2;
2222 HeatedCrossSectionContinuousEnergy &heatedCrossSection2 = *m_heatedCrossSections[temperatureIndex2];
2223 std::size_t energyIndex2 = heatedCrossSection2.evaluationInfo( a_hashIndex, a_energy, &energyFraction2 );
2224
2225 for( sampled_reaction_index = 0; sampled_reaction_index < numberOfReactions; ++sampled_reaction_index ) {
2226 if( m_thresholds[sampled_reaction_index] >= a_energy ) continue;
2227 crossSectionSum += temperatureFraction1 * heatedCrossSection1.reactionCrossSection2( sampled_reaction_index, a_URR_protareInfos,
2228 a_URR_index, a_energy, energyIndex1, energyFraction1 );
2229 crossSectionSum += temperatureFraction2 * heatedCrossSection2.reactionCrossSection2( sampled_reaction_index, a_URR_protareInfos,
2230 a_URR_index, a_energy, energyIndex2, energyFraction2 );
2231 if( crossSectionSum >= sampleCrossSection ) break;
2232 }
2233 }
2234
2235 if( sampled_reaction_index == numberOfReactions ) {
2236 if( crossSectionSum < ( 1.0 - crossSectionSumError ) * a_crossSection ) {
2237#if LUPI_ON_GPU
2238 MCGIDI_PRINTF( "HeatedCrossSectionsContinuousEnergy::sampleReaction: crossSectionSum %.17e less than a_crossSection = %.17e.",
2239 crossSectionSum, a_crossSection );
2240#else
2241 std::string errorString = "HeatedCrossSectionsContinuousEnergy::sampleReaction: crossSectionSum "
2242 + LUPI::Misc::doubleToString3( "%.17e", crossSectionSum ) + " less than a_crossSection = "
2243 + LUPI::Misc::doubleToString3( "%.17e", a_crossSection ) + ".";
2244 LUPI_THROW( errorString.c_str( ) );
2245#endif
2246 }
2247 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.
2248 if( heatedCrossSection1.reactionCrossSection2( sampled_reaction_index, a_URR_protareInfos, a_URR_index, a_energy, energyIndex1,
2249 energyFraction1, true ) > 0 ) break;
2250 }
2251 }
2252
2253 return( sampled_reaction_index );
2254}
2255
2256/* *********************************************************************************************************//**
2257 * Returns the requested reaction's multi-group cross section for target temperature *a_temperature* and projectile multi-group *a_hashIndex*.
2258 *
2259 * @param a_hashIndex [in] The multi-group index.
2260 * @param a_temperature [in] The temperature of the target.
2261 * @param a_energy [in] The energy of the projectile.
2262 * @param a_crossSection [in] The total cross section for the protare at *a_temperature* and *a_energy*.
2263 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2264 ***********************************************************************************************************/
2265
2266template <typename RNG>
2267LUPI_HOST_DEVICE std::size_t MCGIDI::HeatedCrossSectionsMultiGroup::sampleReaction( std::size_t a_hashIndex, double a_temperature, double a_energy, double a_crossSection,
2268 RNG && a_rng ) const {
2269
2270 std::size_t i1, sampled_reaction_index, temperatureIndex1, temperatureIndex2, numberOfTemperatures = m_temperatures.size( );
2271 double sampleCrossSection = a_crossSection * a_rng( );
2272
2273 if( a_temperature <= m_temperatures[0] ) {
2274 temperatureIndex1 = 0;
2275 temperatureIndex2 = temperatureIndex1; }
2276 else if( a_temperature >= m_temperatures.back( ) ) {
2277 temperatureIndex1 = m_temperatures.size( ) - 1;
2278 temperatureIndex2 = temperatureIndex1; }
2279 else {
2280 for( i1 = 0; i1 < numberOfTemperatures; ++i1 ) if( a_temperature < m_temperatures[i1] ) break;
2281 temperatureIndex1 = i1 - 1;
2282 temperatureIndex2 = i1;
2283 }
2284
2285 std::size_t numberOfReactions = m_heatedCrossSections[0]->numberOfReactions( );
2286 double crossSectionSum = 0;
2287 HeatedCrossSectionMultiGroup &heatedCrossSection1 = *m_heatedCrossSections[temperatureIndex1];
2288
2289 if( temperatureIndex1 == temperatureIndex2 ) {
2290 for( sampled_reaction_index = 0; sampled_reaction_index < numberOfReactions; ++sampled_reaction_index ) {
2291 crossSectionSum += heatedCrossSection1.reactionCrossSection( sampled_reaction_index, a_hashIndex, true );
2292 if( crossSectionSum >= sampleCrossSection ) break;
2293 } }
2294 else {
2295 double temperatureFraction2 = ( a_temperature - m_temperatures[temperatureIndex1] ) / ( m_temperatures[temperatureIndex2] - m_temperatures[temperatureIndex1] );
2296 double temperatureFraction1 = 1.0 - temperatureFraction2;
2297 HeatedCrossSectionMultiGroup &heatedCrossSection2 = *m_heatedCrossSections[temperatureIndex2];
2298
2299 for( sampled_reaction_index = 0; sampled_reaction_index < numberOfReactions; ++sampled_reaction_index ) {
2300 if( m_thresholds[sampled_reaction_index] >= a_energy ) continue;
2301 crossSectionSum += temperatureFraction1 * heatedCrossSection1.reactionCrossSection( sampled_reaction_index, a_hashIndex, true );
2302 crossSectionSum += temperatureFraction2 * heatedCrossSection2.reactionCrossSection( sampled_reaction_index, a_hashIndex, true );
2303 if( crossSectionSum >= sampleCrossSection ) break;
2304 }
2305 }
2306
2307 if( sampled_reaction_index == numberOfReactions ) return( MCGIDI_nullReaction );
2308
2309 if( m_multiGroupThresholdIndex[sampled_reaction_index] == static_cast<int>( a_hashIndex ) ) {
2310 double energyAboveThreshold = a_energy - m_thresholds[sampled_reaction_index];
2311
2312 if( energyAboveThreshold <= ( a_rng( ) * ( m_projectileMultiGroupBoundariesCollapsed[a_hashIndex+1] - m_thresholds[sampled_reaction_index] ) ) )
2313 return( MCGIDI_nullReaction );
2314 }
2315
2316 return( sampled_reaction_index );
2317}
2318
2319// From file: MCGIDI_misc.cpp
2320
2321/* *********************************************************************************************************//**
2322 * This function returns a normalized Maxwellian speed (i.e., v = |velocity|) in 3d (i.e., x^2 Exp( -x^2 ))
2323 * where v = sqrt(2 * T / m) * x.
2324 * Using formula in https://link.springer.com/content/pdf/10.1007%2Fs10955-011-0364-y.pdf.
2325 * Author Nader M.A. Mohamed, title "Efficient Algorithm for Generating Maxwell Random Variables".
2326 *
2327 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2328 *
2329 * @return The sampled normalized Maxwellian speed.
2330 ***********************************************************************************************************/
2331
2332template <typename RNG>
2333inline LUPI_HOST_DEVICE double sampleBetaFromMaxwellian( RNG && a_rng ) {
2334
2335 double _g = 2.0 / ( 1.37 * 0.5 * 1.772453850905516 ); // 1.772453850905516 = sqrt( pi ).
2336 double beta, r1;
2337
2338 do {
2339 r1 = a_rng( );
2340 beta = sqrt( -2.0 * log( r1 ) );
2341 } while( _g * r1 * beta < a_rng( ) );
2342
2343 return( beta );
2344}
2345
2346namespace MCGIDI {
2347
2348/* *********************************************************************************************************//**
2349 * This function boost a particle from one frame to another frame. The frames have a relative speed *a_boostSpeed*
2350 * and cosine of angle *a_boostMu* between their z-axes. BRB FIXME, currently it is the x-axis.
2351 *
2352 * @param a_input [in] Instance containing a random number generator that returns a double in the range [0, 1).
2353 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2354 * @param a_product [in] The particle to boost.
2355 ***********************************************************************************************************/
2356
2357template <typename RNG>
2358inline LUPI_HOST_DEVICE void upScatterModelABoostParticle( Sampling::Input &a_input, RNG && a_rng, Sampling::Product &a_product ) {
2359
2360 double C_rel = 1.0;
2361 if( a_input.m_relativeBeta != 0.0 ) {
2362 C_rel = ( a_input.m_projectileBeta - a_input.m_muLab * a_input.m_targetBeta ) / a_input.m_relativeBeta;
2363 if( C_rel > 1.0 ) C_rel = 1.0; // Handle round-off issue. Probably should check how big the issue is.
2364 if( C_rel < -1.0 ) C_rel = -1.0; // Handle round-off issue. Probably should check how big the issue is.
2365 }
2366 double S_rel = sqrt( 1.0 - C_rel * C_rel );
2367
2368 double pz_vz = a_product.m_pz_vz;
2369 a_product.m_pz_vz = C_rel * a_product.m_pz_vz + S_rel * a_product.m_px_vx;
2370 a_product.m_px_vx = -S_rel * pz_vz + C_rel * a_product.m_px_vx;
2371
2372 double targetSpeed = MCGIDI_speedOfLight_cm_sec * a_input.m_targetBeta;
2373 a_product.m_pz_vz += a_input.m_muLab * targetSpeed;
2374 a_product.m_px_vx += sqrt( 1.0 - a_input.m_muLab * a_input.m_muLab ) * targetSpeed;
2375
2376 double phi = 2.0 * M_PI * a_rng( );
2377 double sine = sin( phi );
2378 double cosine = cos( phi );
2379 double px_vx = a_product.m_px_vx;
2380 a_product.m_px_vx = cosine * a_product.m_px_vx - sine * a_product.m_py_vy;
2381 a_product.m_py_vy = sine * px_vx + cosine * a_product.m_py_vy;
2382
2383 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;
2384 speed2 /= MCGIDI_speedOfLight_cm_sec * MCGIDI_speedOfLight_cm_sec;
2385
2386 a_product.m_kineticEnergy = particleKineticEnergyFromBeta2( a_product.m_productMass, speed2 );
2387}
2388
2389}
2390
2391
2392// From file: MCGIDI_outputChannel.cpp
2393
2394
2395/* *********************************************************************************************************//**
2396 * This method adds sampled products to *a_products*.
2397 *
2398 * @param a_protare [in] The Protare this Reaction belongs to.
2399 * @param a_projectileEnergy [in] The energy of the projectile.
2400 * @param a_input [in] Sample options requested by user.
2401 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2402 * @param a_products [in] The object to add all sampled products to.
2403 ***********************************************************************************************************/
2404
2405template <typename RNG, typename PUSHBACK>
2406LUPI_HOST_DEVICE void MCGIDI::OutputChannel::sampleProducts( ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input,
2407 RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) const {
2408
2409 if( m_hasFinalStatePhotons ) {
2410 double random = a_rng( );
2411 double cumulative = 0.0;
2412 bool sampled = false;
2413 for( auto productIter = m_products.begin( ); productIter != m_products.end( ); ++productIter ) {
2414 cumulative += (*productIter)->multiplicity( )->evaluate( a_projectileEnergy );
2415 if( cumulative >= random ) {
2416 (*productIter)->sampleFinalState( a_protare, a_projectileEnergy, a_input, a_rng, a_push_back, a_products );
2417 sampled = true;
2418 break;
2419 }
2420 }
2421 if( !sampled ) { // BRB: FIXME: still need to code for continuum photon.
2422 } }
2423 else {
2424 for( Vector<Product *>::const_iterator iter = m_products.begin( ); iter != m_products.end( ); ++iter )
2425 (*iter)->sampleProducts( a_protare, a_projectileEnergy, a_input, a_rng, a_push_back, a_products );
2426 }
2427
2428 if( m_totalDelayedNeutronMultiplicity != nullptr ) {
2429 double totalDelayedNeutronMultiplicity = m_totalDelayedNeutronMultiplicity->evaluate( a_projectileEnergy );
2430
2431 if( a_rng( ) < totalDelayedNeutronMultiplicity ) { // Assumes that totalDelayedNeutronMultiplicity < 1.0, which it is.
2432 double sum = 0.0;
2433
2434 totalDelayedNeutronMultiplicity *= a_rng( );
2435 for( std::size_t i1 = 0; i1 < (std::size_t) m_delayedNeutrons.size( ); ++i1 ) {
2436 DelayedNeutron const *delayedNeutron1( delayedNeutron( i1 ) );
2437 Product const &product = delayedNeutron1->product( );
2438
2439 sum += product.multiplicity( )->evaluate( a_projectileEnergy );
2440 if( sum >= totalDelayedNeutronMultiplicity ) {
2441 product.distribution( )->sample( a_projectileEnergy, a_input, a_rng );
2442 a_input.m_delayedNeutronIndex = delayedNeutron1->delayedNeutronIndex( );
2443 a_input.m_delayedNeutronDecayRate = delayedNeutron1->rate( );
2444 a_products.add( a_projectileEnergy, product.intid( ), product.index( ), product.userParticleIndex( ), product.mass( ),
2445 a_input, a_rng, a_push_back, false );
2446 break;
2447 }
2448 }
2449 }
2450 }
2451}
2452
2453/* *********************************************************************************************************//**
2454 * Returns the probability for a project with energy *a_energy_in* to cause this channel to emitted a particle of index
2455 * *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.
2456 *
2457 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
2458 * @param a_index [in] The index of the particle to emit.
2459 * @param a_temperature [in] Specifies the temperature of the material.
2460 * @param a_energy_in [in] The energy of the incident particle.
2461 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
2462 * @param a_weight [in] The probability of emitting outgoing particle into lab angle *a_mu_lab*.
2463 * @param a_energy_out [in] The energy of the emitted outgoing particle.
2464 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2465 * @param a_cumulative_weight [in] The sum of the multiplicity for other outgoing particles with index *a_index*.
2466 ***********************************************************************************************************/
2467
2468template <typename RNG>
2469LUPI_HOST_DEVICE void MCGIDI::OutputChannel::angleBiasing( Reaction const *a_reaction, int a_index, double a_temperature, double a_energy_in, double a_mu_lab,
2470 double &a_weight, double &a_energy_out, RNG && a_rng, double &a_cumulative_weight ) const {
2471
2472 for( Vector<Product *>::const_iterator iter = m_products.begin( ); iter != m_products.end( ); ++iter )
2473 (*iter)->angleBiasing( a_reaction, a_index, a_temperature, a_energy_in, a_mu_lab, a_weight, a_energy_out, a_rng, a_cumulative_weight );
2474
2475 if( ( m_totalDelayedNeutronMultiplicity != nullptr ) && ( a_index == m_neutronIndex ) ) {
2476 for( std::size_t i1 = 0; i1 < (std::size_t) m_delayedNeutrons.size( ); ++i1 ) {
2477 DelayedNeutron const *delayedNeutron1( delayedNeutron( i1 ) );
2478 Product const &product = delayedNeutron1->product( );
2479
2480 product.angleBiasing( a_reaction, a_index, a_temperature, a_energy_in, a_mu_lab, a_weight, a_energy_out, a_rng, a_cumulative_weight );
2481 }
2482 }
2483}
2484
2485/* *********************************************************************************************************//**
2486 * Returns the probability for a project with energy *a_energy_in* to cause this channel to emitted a particle of intid
2487 * *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.
2488 *
2489 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
2490 * @param a_intid [in] The intid of the particle to emit.
2491 * @param a_temperature [in] Specifies the temperature of the material.
2492 * @param a_energy_in [in] The energy of the incident particle.
2493 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
2494 * @param a_weight [in] The probability of emitting outgoing particle into lab angle *a_mu_lab*.
2495 * @param a_energy_out [in] The energy of the emitted outgoing particle.
2496 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2497 * @param a_cumulative_weight [in] The sum of the multiplicity for other outgoing particles with intid *a_intid*.
2498 ***********************************************************************************************************/
2499
2500template <typename RNG>
2501LUPI_HOST_DEVICE void MCGIDI::OutputChannel::angleBiasingViaIntid( Reaction const *a_reaction, int a_intid, double a_temperature, double a_energy_in, double a_mu_lab,
2502 double &a_weight, double &a_energy_out, RNG && a_rng, double &a_cumulative_weight ) const {
2503
2504 for( Vector<Product *>::const_iterator iter = m_products.begin( ); iter != m_products.end( ); ++iter )
2505 (*iter)->angleBiasingViaIntid( a_reaction, a_intid, a_temperature, a_energy_in, a_mu_lab, a_weight, a_energy_out, a_rng, a_cumulative_weight );
2506
2507 if( ( m_totalDelayedNeutronMultiplicity != nullptr ) && ( a_intid == PoPI::Intids::neutron ) ) {
2508 for( std::size_t i1 = 0; i1 < (std::size_t) m_delayedNeutrons.size( ); ++i1 ) {
2509 DelayedNeutron const *delayedNeutron1( delayedNeutron( i1 ) );
2510 Product const &product = delayedNeutron1->product( );
2511
2512 product.angleBiasingViaIntid( a_reaction, a_intid, a_temperature, a_energy_in, a_mu_lab, a_weight, a_energy_out, a_rng, a_cumulative_weight );
2513 }
2514 }
2515}
2516
2517
2518// From file: MCGIDI_product.cpp
2519
2520/* *********************************************************************************************************//**
2521 * This method adds sampled products to *a_products*.
2522 *
2523 * @param a_protare [in] The Protare this Reaction belongs to.
2524 * @param a_projectileEnergy [in] The energy of the projectile.
2525 * @param a_input [in] Sample options requested by user.
2526 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2527 * @param a_products [in] The object to add all sampled products to.
2528 ***********************************************************************************************************/
2529
2530template <typename RNG, typename PUSHBACK>
2531LUPI_HOST_DEVICE void MCGIDI::Product::sampleProducts( ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input,
2532 RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) const {
2533
2534#ifdef MCGIDI_USE_OUTPUT_CHANNEL
2535 if( m_outputChannel != nullptr ) {
2536 m_outputChannel->sampleProducts( a_protare, a_projectileEnergy, a_input, a_rng, a_push_back, a_products ); }
2537 else {
2538#endif
2539 if( m_twoBodyOrder == TwoBodyOrder::secondParticle ) {
2540 a_products.add( a_projectileEnergy, intid( ), index( ), userParticleIndex( ), mass( ), a_input, a_rng, a_push_back, m_intid == PoPI::Intids::photon ); }
2541 else {
2542 int _multiplicity = m_multiplicity->sampleBoundingInteger( a_projectileEnergy, a_rng );
2543 int __multiplicity = _multiplicity;
2544
2545 for( ; _multiplicity > 0; --_multiplicity ) {
2546 m_distribution->sample( a_projectileEnergy, a_input, a_rng );
2547 a_input.m_delayedNeutronIndex = -1;
2548 a_input.m_delayedNeutronDecayRate = 0.0;
2549 a_products.add( a_projectileEnergy, intid( ), index( ), userParticleIndex( ), mass( ), a_input, a_rng, a_push_back, m_intid == PoPI::Intids::photon );
2550 }
2551 if( m_initialStateIndex >= 0 ) {
2552 if( __multiplicity == 0 ) {
2553 a_protare->sampleBranchingGammas( a_input, a_projectileEnergy, m_initialStateIndex, a_rng, a_push_back, a_products );
2554 }
2555 }
2556 }
2557#ifdef MCGIDI_USE_OUTPUT_CHANNEL
2558 }
2559#endif
2560}
2561
2562/* *********************************************************************************************************//**
2563 * This method adds sampled products to *a_products*. In particular, the product is a capture reaction
2564 * primary gamma what has a finalState attribute. This gamma is added as well as the gammas from the
2565 * gamma cascade.
2566 *
2567 * @param a_protare [in] The Protare this Reaction belongs to.
2568 * @param a_projectileEnergy [in] The energy of the projectile.
2569 * @param a_input [in] Sample options requested by user.
2570 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2571 * @param a_products [in] The object to add all sampled products to.
2572 ***********************************************************************************************************/
2573template <typename RNG, typename PUSHBACK>
2574LUPI_HOST_DEVICE void MCGIDI::Product::sampleFinalState( ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input,
2575 RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) const {
2576
2577 m_distribution->sample( a_projectileEnergy, a_input, a_rng );
2578 a_input.m_delayedNeutronIndex = -1;
2579 a_input.m_delayedNeutronDecayRate = 0.0;
2580 a_products.add( a_projectileEnergy, m_intid, m_index, m_userParticleIndex, mass( ), a_input, a_rng, a_push_back, m_intid == PoPI::Intids::photon );
2581
2582 if( m_initialStateIndex >= 0 ) {
2583 a_protare->sampleBranchingGammas( a_input, a_projectileEnergy, m_initialStateIndex, a_rng, a_push_back, a_products );
2584 }
2585}
2586
2587/* *********************************************************************************************************//**
2588 * Returns the weight for a projectile with energy *a_energy_in* to cause this channel to emitted a particle of index
2589 * *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.
2590 *
2591 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
2592 * @param a_index [in] The index of the particle to emit.
2593 * @param a_temperature [in] Specifies the temperature of the material.
2594 * @param a_energy_in [in] The energy of the incident particle.
2595 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
2596 * @param a_weight [in] The weight of emitting outgoing particle into lab angle *a_mu_lab*.
2597 * @param a_energy_out [in] The energy of the emitted outgoing particle.
2598 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2599 * @param a_cumulative_weight [in] The sum of the multiplicity for other outgoing particles with index *a_index*.
2600 ***********************************************************************************************************/
2601template <typename RNG>
2602LUPI_HOST_DEVICE void MCGIDI::Product::angleBiasing( Reaction const *a_reaction, int a_index, double a_temperature, double a_energy_in, double a_mu_lab,
2603 double &a_weight, double &a_energy_out, RNG && a_rng, double &a_cumulative_weight ) const {
2604
2605#ifdef MCGIDI_USE_OUTPUT_CHANNEL
2606 if( m_outputChannel != nullptr ) {
2607 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 ); }
2608 else {
2609#endif
2610 if( m_index != a_index ) return;
2611
2612 double probability = 0.0;
2613 double energy_out = 0.0;
2614
2615 if( a_cumulative_weight == 0.0 ) a_energy_out = 0.0;
2616
2617 if( m_multiplicity->type( ) == Function1dType::branching ) { // Needs to handle F1_Branching.
2618 }
2619 else {
2620 probability = m_distribution->angleBiasing( a_reaction, a_temperature, a_energy_in, a_mu_lab, a_rng, energy_out );
2621 }
2622
2623 double weight = m_multiplicity->evaluate( a_energy_in ) * probability;
2624 a_cumulative_weight += weight;
2625 if( weight > a_rng( ) * a_cumulative_weight ) {
2626 a_weight = weight;
2627 a_energy_out = energy_out;
2628 }
2629#ifdef MCGIDI_USE_OUTPUT_CHANNEL
2630 }
2631#endif
2632}
2633
2634/* *********************************************************************************************************//**
2635 * Returns the weight for a projectile with energy *a_energy_in* to cause this channel to emitted a particle of intid
2636 * *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.
2637 *
2638 * @param a_reaction [in] The reaction containing the particle which this distribution describes.
2639 * @param a_intid [in] The intid of the particle to emit.
2640 * @param a_temperature [in] Specifies the temperature of the material.
2641 * @param a_energy_in [in] The energy of the incident particle.
2642 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
2643 * @param a_weight [in] The weight of emitting outgoing particle into lab angle *a_mu_lab*.
2644 * @param a_energy_out [in] The energy of the emitted outgoing particle.
2645 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2646 * @param a_cumulative_weight [in] The sum of the multiplicity for other outgoing particles with intid *a_intid*.
2647 ***********************************************************************************************************/
2648
2649template <typename RNG>
2650LUPI_HOST_DEVICE void MCGIDI::Product::angleBiasingViaIntid( Reaction const *a_reaction, int a_intid, double a_temperature, double a_energy_in, double a_mu_lab,
2651 double &a_weight, double &a_energy_out, RNG && a_rng , double &a_cumulative_weight ) const {
2652
2653#ifdef MCGIDI_USE_OUTPUT_CHANNEL
2654 if( m_outputChannel != nullptr ) {
2655 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 ); }
2656 else {
2657#endif
2658 if( m_intid != a_intid ) return;
2659
2660 angleBiasing( a_reaction, m_index, a_temperature, a_energy_in, a_mu_lab, a_weight, a_energy_out, a_rng, a_cumulative_weight );
2661
2662#ifdef MCGIDI_USE_OUTPUT_CHANNEL
2663 }
2664#endif
2665}
2666
2667
2668
2669// From file: MCGIDI_protare.cpp
2670
2671
2672/* *********************************************************************************************************//**
2673 * Samples a reaction of *this* and returns its index.
2674 *
2675 * @param a_input [in/out] Sample options requested by user.
2676 * @param a_URR_protareInfos [in] URR information.
2677 * @param a_hashIndex [in] Specifies the continuous energy hash index or multi-group index.
2678 * @param a_crossSection [in] The total cross section for the protare at *a_input.temperature()* and *a_input.energy()*.
2679 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2680 *
2681 * @return The index of the sampled reaction.
2682 ***********************************************************************************************************/
2683
2684template <typename RNG>
2685LUPI_HOST_DEVICE std::size_t MCGIDI::Protare::sampleReaction( Sampling::Input &a_input, URR_protareInfos const &a_URR_protareInfos,
2686 std::size_t a_hashIndex, double a_crossSection, RNG && a_rng ) const {
2687
2688 std::size_t reactionIndex = MCGIDI_nullReaction;
2689
2690 switch( protareType( ) ) {
2691 case ProtareType::single:
2692 reactionIndex = static_cast<ProtareSingle const *>( this )->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex,
2693 a_crossSection, a_rng );
2694 break;
2695 case ProtareType::composite:
2696 reactionIndex = static_cast<ProtareComposite const *>( this )->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex,
2697 a_crossSection, a_rng );
2698 break;
2699 case ProtareType::TNSL:
2700 reactionIndex = static_cast<ProtareTNSL const *>( this )->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex,
2701 a_crossSection, a_rng );
2702 break;
2703 }
2704
2705 return( reactionIndex );
2706}
2707
2708/* *********************************************************************************************************//**
2709 * Samples gammas from a nuclide electro-magnetic decay.
2710 *
2711 * @param a_input [in] Sample options requested by user.
2712 * @param a_projectileEnergy [in] The energy of the projectile.
2713 * @param a_initialStateIndex [in] The index in *m_nuclideGammaBranchStateInfos* whose nuclide data are used for sampling.
2714 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2715 * @param a_products [in] The object to add all sampled gammas to.
2716 ***********************************************************************************************************/
2717template <typename RNG, typename PUSHBACK>
2718LUPI_HOST_DEVICE void MCGIDI::ProtareSingle::sampleBranchingGammas( Sampling::Input &a_input, double a_projectileEnergy, int a_initialStateIndex,
2719 RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) const {
2720
2721 int initialStateIndex = a_initialStateIndex;
2722 double energyLevelSampleWidthUpper = 0.0; // Used for GRIN continuum levels to add variaction to outgoing photons.
2723
2724 NuclideGammaBranchStateInfo *nuclideGammaBranchStateInfo = nullptr;
2725 if( initialStateIndex >= 0 ) nuclideGammaBranchStateInfo = m_nuclideGammaBranchStateInfos[static_cast<std::size_t>(initialStateIndex)];
2726 while( initialStateIndex >= 0 ) {
2727 auto const &branchIndices = nuclideGammaBranchStateInfo->branchIndices( );
2728
2729 double random = a_rng( );
2730 double sum = 0.0;
2731 initialStateIndex = -1; // Just in case the for loop never has "sum >= random".
2732 for( std::size_t i1 = 0; i1 < branchIndices.size( ); ++i1 ) {
2733 NuclideGammaBranchInfo *nuclideGammaBranchInfo = m_branches[branchIndices[i1]];
2734
2735 sum += nuclideGammaBranchInfo->probability( );
2736 if( sum >= random ) {
2737 double energyLevelSampleWidthLower = 0.0;
2738 initialStateIndex = nuclideGammaBranchInfo->residualStateIndex( );
2739 if( initialStateIndex >= 0 ) {
2740 nuclideGammaBranchStateInfo = m_nuclideGammaBranchStateInfos[static_cast<std::size_t>(initialStateIndex)];
2741 energyLevelSampleWidthLower = a_rng( ) * nuclideGammaBranchStateInfo->nuclearLevelEnergyWidth( );
2742 }
2743 if( nuclideGammaBranchInfo->photonEmissionProbability( ) > a_rng( ) ) {
2744 a_input.setSampledType( Sampling::SampledType::photon );
2745 a_input.m_dataInTargetFrame = false;
2746 a_input.m_frame = GIDI::Frame::lab;
2747
2748 a_input.m_energyOut1 = nuclideGammaBranchInfo->gammaEnergy( ) + energyLevelSampleWidthUpper - energyLevelSampleWidthLower;
2749 a_input.m_mu = 1.0 - 2.0 * a_rng( );
2750 a_input.m_phi = 2.0 * M_PI * a_rng( );
2751
2752 a_products.add( a_projectileEnergy, PoPI::Intids::photon, m_photonIndex, userPhotonIndex( ), 0.0, a_input, a_rng, a_push_back, true );
2753 }
2754 energyLevelSampleWidthUpper = energyLevelSampleWidthLower;
2755 break;
2756 }
2757 }
2758 }
2759}
2760
2761/* *********************************************************************************************************//**
2762 * This function is used internally to sample a target's velocity (speed and cosine of angle relative to projectile)
2763 * for a heated target using zero temperature, multi-grouped cross sections.
2764 *
2765 * @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*.
2766 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2767 *
2768 * @return Returns *true* if target velocity is sampled and false otherwise.
2769 ***********************************************************************************************************/
2770
2771template <typename RNG>
2772inline LUPI_HOST_DEVICE bool MCGIDI::ProtareSingle::sampleTargetBetaForUpscatterModelA( Sampling::Input &a_input, RNG && a_rng ) const {
2773
2774 double projectileBeta = MCGIDI_particleBeta( m_projectileMass, a_input.energy( ) );
2775 double targetThermalBeta = MCGIDI_particleBeta( m_targetMass, a_input.temperature( ) );
2776
2777 a_input.m_projectileBeta = projectileBeta;
2778 a_input.m_relativeBeta = projectileBeta;
2779 a_input.m_muLab = 0.0;
2780 a_input.m_targetBeta = 0.0;
2781
2782 if( targetThermalBeta < 1e-4 * projectileBeta ) return( false );
2783
2784 a_input.m_modelTemperature = 0.0;
2785
2786 double relativeBetaMin = projectileBeta - 2.0 * targetThermalBeta;
2787 double relativeBetaMax = projectileBeta + 2.0 * targetThermalBeta;
2788
2789 std::size_t maxIndex = m_upscatterModelAGroupVelocities.size( ) - 2;
2790 int intRelativeBetaMinIndex = binarySearchVector( relativeBetaMin, m_upscatterModelAGroupVelocities, true );
2791 std::size_t relativeBetaMinIndex = static_cast<std::size_t>( intRelativeBetaMinIndex );
2792 int intRelativeBetaMaxIndex = binarySearchVector( relativeBetaMax, m_upscatterModelAGroupVelocities, true );
2793 std::size_t relativeBetaMaxIndex = static_cast<std::size_t>( intRelativeBetaMaxIndex );
2794 double targetBeta, relativeBeta, mu;
2795
2796 if( relativeBetaMinIndex >= maxIndex ) relativeBetaMinIndex = maxIndex;
2797 if( relativeBetaMaxIndex >= maxIndex ) relativeBetaMaxIndex = maxIndex;
2798
2799 if( relativeBetaMinIndex == relativeBetaMaxIndex ) {
2800 targetBeta = targetThermalBeta * sampleBetaFromMaxwellian( a_rng );
2801 mu = 1.0 - 2.0 * a_rng( );
2802 relativeBeta = sqrt( targetBeta * targetBeta + projectileBeta * projectileBeta - 2.0 * mu * targetBeta * projectileBeta ); }
2803 else {
2804
2805 double reactionRate;
2806 double reactionRateMax = 0;
2807 for( std::size_t i1 = relativeBetaMinIndex; i1 <= relativeBetaMaxIndex; ++i1 ) {
2808 reactionRate = m_upscatterModelACrossSection[i1] * m_upscatterModelAGroupVelocities[i1+1];
2809 if( reactionRate > reactionRateMax ) reactionRateMax = reactionRate;
2810 }
2811
2812 do {
2813 targetBeta = targetThermalBeta * sampleBetaFromMaxwellian( a_rng );
2814 mu = 1.0 - 2.0 * a_rng( );
2815 relativeBeta = sqrt( targetBeta * targetBeta + projectileBeta * projectileBeta - 2.0 * mu * targetBeta * projectileBeta );
2816
2817 std::size_t index = static_cast<std::size_t>( binarySearchVector( relativeBeta, m_upscatterModelAGroupVelocities, true ) );
2818 if( index > maxIndex ) index = maxIndex;
2819 reactionRate = m_upscatterModelACrossSection[index] * relativeBeta;
2820 } while( reactionRate < a_rng( ) * reactionRateMax );
2821 }
2822
2823 a_input.m_modelEnergy = particleKineticEnergy( m_projectileMass, relativeBeta );
2824 a_input.m_relativeBeta = relativeBeta;
2825 a_input.m_muLab = mu;
2826 a_input.m_targetBeta = targetBeta;
2827
2828 return( true );
2829}
2830
2831/* *********************************************************************************************************//**
2832 * Returns the index of a sampled reaction for target temperature, projectile energy and total cross section
2833 * as specified via argument *a_input*. Random numbers are obtained via *a_rng*.
2834 *
2835 * @param a_input [in/out] Sample options requested by user. The values m_modelTemperature and m_modelEnergy are set by this method.
2836 * @param a_URR_protareInfos [in] URR information.
2837 * @param a_hashIndex [in] Specifies the continuous energy hash index or multi-group index.
2838 * @param a_crossSection [in] The total cross section for the protare at *a_input.temperature()* and *a_input.energy()*.
2839 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2840 ***********************************************************************************************************/
2841
2842template <typename RNG>
2843LUPI_HOST_DEVICE std::size_t MCGIDI::ProtareSingle::sampleReaction( Sampling::Input &a_input, URR_protareInfos const &a_URR_protareInfos,
2844 std::size_t a_hashIndex, double a_crossSection, RNG && a_rng ) const {
2845
2846 std::size_t hashIndex = a_hashIndex;
2847 double crossSection1 = a_crossSection;
2848
2849 a_input.m_dataInTargetFrame = false;
2850 a_input.m_modelTemperature = a_input.m_temperature;
2851 a_input.m_modelEnergy = a_input.m_energy;
2852
2853 if( upscatterModelASupported( ) && ( a_input.m_upscatterModel == Sampling::Upscatter::Model::A ) ) {
2854 a_input.m_dataInTargetFrame = sampleTargetBetaForUpscatterModelA( a_input, a_rng );
2855 if( a_input.m_dataInTargetFrame ) {
2856 if( m_continuousEnergy ) {
2857 hashIndex = m_domainHash.index( a_input.m_modelEnergy ); }
2858 else {
2859 hashIndex = m_multiGroupHash.index( a_input.m_modelEnergy );
2860 }
2861 crossSection1 = crossSection( a_URR_protareInfos, hashIndex, a_input.m_modelTemperature, a_input.m_modelEnergy, true );
2862 }
2863 }
2864
2865 if( m_continuousEnergy ) return( m_heatedCrossSections.sampleReaction( a_URR_protareInfos, m_URR_index, hashIndex,
2866 a_input.m_modelTemperature, a_input.m_modelEnergy, crossSection1, a_rng ) );
2867
2868 return( m_heatedMultigroupCrossSections.sampleReaction( hashIndex, a_input.m_modelTemperature, a_input.m_modelEnergy,
2869 crossSection1, a_rng ) );
2870}
2871
2872
2873// From file: MCGIDI_protareComposite.cpp
2874
2875/* *********************************************************************************************************//**
2876 * Samples a reaction of *this* and returns its index.
2877 *
2878 * @param a_input [in] Sample options requested by user.
2879 * @param a_URR_protareInfos [in] URR information.
2880 * @param a_hashIndex [in] Specifies the continuous energy hash index or multi-group index.
2881 * @param a_crossSection [in] The total cross section for the protare at *a_input.temperature()* and *a_input.energy()*.
2882 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2883 *
2884 * @return The index of the sampled reaction.
2885 ***********************************************************************************************************/
2886
2887template <typename RNG>
2888LUPI_HOST_DEVICE std::size_t MCGIDI::ProtareComposite::sampleReaction( Sampling::Input &a_input, URR_protareInfos const &a_URR_protareInfos,
2889 std::size_t a_hashIndex, double a_crossSection, RNG && a_rng ) const {
2890
2891 std::size_t length = static_cast<std::size_t>( m_protares.size( ) );
2892 std::size_t reaction_index = 0;
2893 double cross_section_sum = 0.0;
2894 double cross_section_rng = a_rng( ) * a_crossSection;
2895
2896 for( std::size_t i1 = 0; i1 < length; ++i1 ) {
2897 double cross_section = m_protares[i1]->crossSection( a_URR_protareInfos, a_hashIndex, a_input.temperature( ), a_input.energy( ), true );
2898
2899 cross_section_sum += cross_section;
2900 if( cross_section_sum > cross_section_rng ) {
2901 std::size_t reaction_index2 = m_protares[i1]->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex, cross_section, a_rng );
2902
2903 reaction_index += reaction_index2;
2904 if( reaction_index2 == MCGIDI_nullReaction ) reaction_index = MCGIDI_nullReaction;
2905 break;
2906 }
2907 reaction_index += m_protares[i1]->numberOfReactions( );
2908 }
2909
2910 return( reaction_index );
2911}
2912
2913
2914// From file: MCGIDI_protareTNSL.cpp
2915
2916/* *********************************************************************************************************//**
2917 * Returns the total cross section.
2918 *
2919 * @param a_input [in] Sample options requested by user.
2920 * @param a_URR_protareInfos [in] URR information.
2921 * @param a_hashIndex [in] Specifies the continuous energy hash index or multi-group index.
2922 * @param a_crossSection [in] The total cross section for the protare at *a_input.temperature()* and *a_input.energy()*.
2923 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2924 *
2925 * @return The index of the sampled reaction.
2926 ***********************************************************************************************************/
2927
2928template <typename RNG>
2929LUPI_HOST_DEVICE std::size_t MCGIDI::ProtareTNSL::sampleReaction( Sampling::Input &a_input, URR_protareInfos const &a_URR_protareInfos,
2930 std::size_t a_hashIndex, double a_crossSection, RNG && a_rng ) const {
2931
2932 std::size_t reactionIndex = 0;
2933
2934 if( ( a_input.energy( ) < m_TNSL_maximumEnergy ) && ( a_input.temperature( ) <= m_TNSL_maximumTemperature ) ) {
2935 double TNSL_crossSection = m_TNSL->crossSection( a_URR_protareInfos, a_hashIndex, a_input.temperature( ), a_input.energy( ), true );
2936
2937 if( TNSL_crossSection > a_rng( ) * a_crossSection ) {
2938 reactionIndex = m_TNSL->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex, TNSL_crossSection, a_rng ); }
2939 else {
2940 reactionIndex = m_protareWithoutElastic->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex,
2941 a_crossSection - TNSL_crossSection, a_rng );
2942 if( reactionIndex != MCGIDI_nullReaction ) reactionIndex += m_numberOfTNSLReactions + 1;
2943 } }
2944 else {
2945 reactionIndex = m_protareWithElastic->sampleReaction( a_input, a_URR_protareInfos, a_hashIndex, a_crossSection, a_rng );
2946 if( reactionIndex != MCGIDI_nullReaction ) reactionIndex += m_numberOfTNSLReactions;
2947 }
2948
2949 return( reactionIndex );
2950}
2951
2952
2953// From file: MCGIDI_reaction.cpp
2954
2955/* *********************************************************************************************************//**
2956 * This method adds sampled products to *a_products*.
2957 *
2958 * @param a_protare [in] The Protare this Reaction belongs to.
2959 * @param a_input [in] Sample options requested by user.
2960 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
2961 * @param a_products [in] The object to add all sampled products to.
2962 * @param a_checkOrphanProducts [in] If true, associated orphan products are also sampled.
2963 ***********************************************************************************************************/
2964
2965template <typename RNG, typename PUSHBACK>
2966LUPI_HOST_DEVICE void MCGIDI::Reaction::sampleProducts( Protare const *a_protare, Sampling::Input &a_input,
2967 RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products, bool a_checkOrphanProducts ) const {
2968
2969 double projectileEnergy = a_input.modelEnergy( );
2970
2971 a_input.m_GRIN_intermediateResidual = -1;
2972 a_input.m_reaction = this;
2973 a_input.m_projectileMass = m_projectileMass;
2974 a_input.m_targetMass = m_targetMass;
2975 a_input.m_relativeBeta = MCGIDI_particleBeta( m_projectileMass, projectileEnergy );
2976
2977 if( m_GRIN_specialSampleProducts ) {
2978 if( m_GRIN_capture != nullptr ) {
2979 if( projectileEnergy < m_GRIN_maximumCaptureIncidentEnergy ) {
2980 if( m_GRIN_capture->sampleProducts( (ProtareSingle const *) a_protare, projectileEnergy, a_input, a_rng, a_push_back, a_products ) ) {
2981 return;
2982 }
2983 } }
2984 else if( m_GRIN_inelastic != nullptr ) {
2985 if( m_GRIN_inelastic->sampleProducts( (ProtareSingle const *) a_protare, projectileEnergy, a_input, a_rng, a_push_back, a_products ) ) {
2986 return;
2987 }
2988 }
2989 }
2990
2991#ifdef MCGIDI_USE_OUTPUT_CHANNEL
2992 m_outputChannel->sampleProducts( m_protareSingle, projectileEnergy, a_input, a_rng, a_push_back, a_products );
2993#else
2994 if( m_hasFinalStatePhotons ) {
2995 double random = a_rng( );
2996 double cumulative = 0.0;
2997 bool sampled = false;
2998 for( auto productIter = m_products.begin( ); productIter != m_products.end( ); ++productIter ) {
2999 cumulative += (*productIter)->multiplicity( )->evaluate( projectileEnergy );
3000 if( cumulative >= random ) {
3001 (*productIter)->sampleFinalState( m_protareSingle, projectileEnergy, a_input, a_rng, a_push_back, a_products );
3002 sampled = true;
3003 break;
3004 }
3005 }
3006 if( !sampled ) { // BRB: FIXME: still need to code for continuum photon.
3007 } }
3008 else {
3009 for( auto productIter = m_products.begin( ); productIter != m_products.end( ); ++productIter ) {
3010 (*productIter)->sampleProducts( m_protareSingle, projectileEnergy, a_input, a_rng, a_push_back, a_products );
3011 }
3012 }
3013
3014 if( m_totalDelayedNeutronMultiplicity != nullptr ) {
3015 double totalDelayedNeutronMultiplicity = m_totalDelayedNeutronMultiplicity->evaluate( projectileEnergy );
3016
3017 if( a_rng( ) < totalDelayedNeutronMultiplicity ) { // Assumes that totalDelayedNeutronMultiplicity < 1.0, which it is.
3018 double sum = 0.0;
3019
3020 totalDelayedNeutronMultiplicity *= a_rng( );
3021 for( std::size_t i1 = 0; i1 < (std::size_t) m_delayedNeutrons.size( ); ++i1 ) {
3022 DelayedNeutron const *delayedNeutron1 = m_delayedNeutrons[i1];
3023 Product const &product = delayedNeutron1->product( );
3024
3025 sum += product.multiplicity( )->evaluate( projectileEnergy );
3026 if( sum >= totalDelayedNeutronMultiplicity ) {
3027 product.distribution( )->sample( projectileEnergy, a_input, a_rng );
3028 a_input.m_delayedNeutronIndex = delayedNeutron1->delayedNeutronIndex( );
3029 a_input.m_delayedNeutronDecayRate = delayedNeutron1->rate( );
3030 a_products.add( a_input.energy( ), product.intid( ), product.index( ), product.userParticleIndex( ), product.mass( ), a_input, a_rng, a_push_back, false );
3031 break;
3032 }
3033 }
3034 }
3035 }
3036
3037 if( m_fissionResiduaIntid != -1 ) { // Special treatment to add 2 ENDL 99120 or 99125 products.
3038 a_input.setSampledType( MCGIDI::Sampling::SampledType::unspecified );
3039 a_input.m_frame = GIDI::Frame::lab;
3040 a_input.m_energyOut1 = 0.0;
3041 a_input.m_mu = 0.0;
3042 a_input.m_phi = 0.0;
3043 a_input.m_delayedNeutronIndex = -1;
3044 a_input.m_delayedNeutronDecayRate = 0.0;
3045 a_products.add( 0.0, m_fissionResiduaIntid, m_fissionResiduaIndex, m_fissionResiduaUserIndex, m_fissionResidualMass, a_input, a_rng, a_push_back, false );
3046 a_products.add( 0.0, m_fissionResiduaIntid, m_fissionResiduaIndex, m_fissionResiduaUserIndex, m_fissionResidualMass, a_input, a_rng, a_push_back, false );
3047 }
3048
3049#endif
3050
3051 if( a_checkOrphanProducts ) {
3052 for( auto productIter = m_associatedOrphanProducts.begin( ); productIter != m_associatedOrphanProducts.end( ); ++productIter ) {
3053 (*productIter)->sampleProducts( m_protareSingle, projectileEnergy, a_input, a_rng, a_push_back, a_products );
3054 }
3055 }
3056}
3057
3058/* *********************************************************************************************************//**
3059 * This method adds sampled products to *a_products*.
3060 *
3061 * @param a_protare [in] The ProtareSingle this Reaction belongs to.
3062 * @param a_projectileEnergy [in] The energy of the projectile.
3063 * @param a_input [in] Sample options requested by user.
3064 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
3065 * @param a_products [in] The object to add all sampled products to.
3066 ***********************************************************************************************************/
3067
3068template <typename RNG, typename PUSHBACK>
3069LUPI_HOST_DEVICE bool MCGIDI::GRIN_capture::sampleProducts( ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input,
3070 RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) const {
3071
3072 std::size_t index = 0;
3073 double random = a_rng( );
3074 for( ; index < m_summedProbabilities.size( ) - 1; ++index ) {
3075 if( random < m_summedProbabilities[index] ) break;
3076 }
3077 GRIN_captureLevelProbability *GRIN_captureLevelProbability1 = m_captureLevelProbabilities[index];
3078
3079 double availableEnergy = m_captureNeutronSeparationEnergy + a_projectileEnergy;
3080 int primaryCaptureLevelIndex = GRIN_captureLevelProbability1->sampleCaptureLevel( a_protare, availableEnergy, a_rng );
3081 NuclideGammaBranchStateInfo const *nuclideGammaBranchStateInfo = a_protare->nuclideGammaBranchStateInfos( )[static_cast<std::size_t>(primaryCaptureLevelIndex)];
3082
3083 a_input.m_GRIN_intermediateResidual = nuclideGammaBranchStateInfo->intid( );
3084
3085 a_input.setSampledType( Sampling::SampledType::photon );
3086 a_input.m_dataInTargetFrame = false;
3087 a_input.m_frame = GIDI::Frame::lab;
3088
3089 a_input.m_energyOut1 = availableEnergy - nuclideGammaBranchStateInfo->nuclearLevelEnergy( );
3090 a_input.m_mu = 2 * a_rng( ) - 1.0;
3091 a_input.m_phi = 2.0 * M_PI * a_rng( );
3092
3093 a_products.add( a_projectileEnergy, PoPI::Intids::photon, a_protare->photonIndex( ), a_protare->userPhotonIndex( ),
3094 0.0, a_input, a_rng, a_push_back, true );
3095
3096 a_protare->sampleBranchingGammas( a_input, a_projectileEnergy, primaryCaptureLevelIndex, a_rng, a_push_back, a_products );
3097
3098 if( m_residualIntid != -1 ) {
3099 a_input.m_energyOut1 = 0.0;
3100 a_input.m_mu = 0.0;
3101 a_input.m_phi = 0.0;
3102 a_products.add( a_projectileEnergy, m_residualIntid, m_residualIndex, m_residualUserIndex, m_residualMass, a_input, a_rng, a_push_back, false );
3103 }
3104
3105 return( true );
3106}
3107
3108/* *********************************************************************************************************//**
3109 * This method adds sampled products to *a_products*.
3110 *
3111 * @param a_protare [in] The ProtareSingle this Reaction belongs to.
3112 * @param a_projectileEnergy [in] The energy of the projectile.
3113 * @param a_input [in] Sample options requested by user.
3114 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
3115 * @param a_products [in] The object to add all sampled products to.
3116 ***********************************************************************************************************/
3117
3118template <typename RNG, typename PUSHBACK>
3119LUPI_HOST_DEVICE bool MCGIDI::GRIN_inelastic::sampleProducts( ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input,
3120 RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) const {
3121
3122 std::size_t index = 1;
3123 for( ; index < m_energies.size( ); ++index ) {
3124 if( m_energies[index] > a_projectileEnergy ) break;
3125 }
3126 --index;
3127
3128 GRIN_inelasticForEnergy *inelasticForEnergy = m_inelasticForEnergy[index];
3129 int levelIndex = inelasticForEnergy->sampleLevelIndex( a_projectileEnergy, a_rng( ) );
3130 if( levelIndex < 0 ) return( false );
3131
3132 NuclideGammaBranchStateInfo const *nuclideGammaBranchStateInfo = a_protare->nuclideGammaBranchStateInfos( )[static_cast<std::size_t>(levelIndex)];
3133
3134 a_input.m_GRIN_intermediateResidual = nuclideGammaBranchStateInfo->intid( );
3135
3136 double residualMass = m_targetMass + nuclideGammaBranchStateInfo->nuclearLevelEnergy( );
3137 double initialMass = m_neutronMass + m_targetMass;
3138 double finalMass = m_neutronMass + residualMass;
3139 double twoBodyThreshold = 0.5 * ( finalMass * finalMass - initialMass * initialMass ) / m_targetMass;
3140 double betaBoast = sqrt( a_projectileEnergy * ( a_projectileEnergy + 2. * m_neutronMass ) )
3141 / ( a_projectileEnergy + m_neutronMass + m_targetMass ); // betaBoast = v/c.
3142 double _x = m_targetMass * ( a_projectileEnergy - twoBodyThreshold ) / ( finalMass * finalMass );
3143 if( _x < 0 ) _x = 0.; // FIXME There needs to be a better test here.
3144 double Kp;
3145 if( _x < 2e-5 ) {
3146 Kp = finalMass * _x * ( 1 - 0.5 * _x * ( 1 - _x ) ); }
3147 else { // This is the relativistic formula derived from E^2 - (pc)^2 is frame independent.
3148 Kp = sqrt( finalMass * finalMass + 2 * m_targetMass * ( a_projectileEnergy - twoBodyThreshold ) ) - finalMass;
3149 }
3150 if( Kp < 0 ) Kp = 0.; // FIXME There needs to be a better test here.
3151
3152 a_input.setSampledType( Sampling::SampledType::firstTwoBody );
3153 a_input.m_mu = 1.0 - 2.0 * a_rng( );
3154 a_input.m_phi = 2. * M_PI * a_rng( );
3155 kinetics_COMKineticEnergy2LabEnergyAndMomentum( betaBoast, Kp, m_neutronMass, residualMass, a_input );
3156
3157 a_input.m_delayedNeutronIndex = -1;
3158 a_input.m_delayedNeutronDecayRate = 0.0;
3159 a_products.add( a_projectileEnergy, PoPI::Intids::neutron, m_neutronIndex, m_neutronUserParticleIndex, m_neutronMass, a_input, a_rng,
3160 a_push_back, false );
3161 a_products.add( a_projectileEnergy, m_targetIntid, m_targetIndex, m_targetUserParticleIndex,
3162 m_targetMass, a_input, a_rng, a_push_back, false );
3163
3164 a_protare->sampleBranchingGammas( a_input, a_projectileEnergy, levelIndex, a_rng, a_push_back, a_products );
3165
3166 return( true );
3167}
3168
3169/* *********************************************************************************************************//**
3170 * This method samples a capture state level and returns an index into the a_protare->m_nuclideGammaBranchStateInfos vector
3171 * of the sampled state level.
3172 *
3173 * @param a_energy [in] The neutron separation energy plus the projectile energy.
3174 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
3175 *
3176 * @return An integer of the sampled state in a_protare->m_nuclideGammaBranchStateInfos.
3177 ***********************************************************************************************************/
3178
3179template <typename RNG>
3180LUPI_HOST_DEVICE int MCGIDI::GRIN_captureToCompound::sampleCaptureLevel( ProtareSingle const *a_protare, double a_energy,
3181 RNG && a_rng, bool a_checkEnergy ) const {
3182
3183 if( a_checkEnergy ) {
3184 NuclideGammaBranchStateInfo *nuclideGammaBranchStateInfo = a_protare->nuclideGammaBranchStateInfos( )[m_index];
3185 if( nuclideGammaBranchStateInfo->nuclearLevelEnergy( ) < a_energy ) return( -1 );
3186 }
3187
3188 double random = a_rng( );
3189 std::size_t index = 0;
3190 for( ; index < m_continuumIndices.m_levels.size( ) - 1; ++index ) {
3191 if( m_continuumIndices.m_summedProbabilities[index] >= random ) break;
3192 }
3193 return( m_continuumIndices.m_levels[index] );
3194}
3195
3196/* *********************************************************************************************************//**
3197 * This method samples a capture state level and returns an index into the a_protare->m_nuclideGammaBranchStateInfos vector
3198 * of the sampled state level.
3199 *
3200 * @param a_protare [in] The ProtareSingle this Reaction belongs to.
3201 * @param a_energy [in] The neutron separation energy plus the projectile energy.
3202 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
3203 *
3204 * @return An integer of the sampled state in a_protare->m_nuclideGammaBranchStateInfos.
3205 ***********************************************************************************************************/
3206
3207template <typename RNG>
3208LUPI_HOST_DEVICE int MCGIDI::GRIN_captureLevelProbability::sampleCaptureLevel( ProtareSingle const *a_protare, double a_energy, RNG && a_rng ) {
3209
3210 double random = a_rng( );
3211 for( std::size_t index = 0; index < m_knownLevelsAndProbabilities.m_levels.size( ); ++index ) {
3212 if( m_knownLevelsAndProbabilities.m_summedProbabilities[index] >= random ) {
3213 return( m_knownLevelsAndProbabilities.m_levels[index] );
3214 }
3215 }
3216
3217 for( std::size_t i1 = 0; i1 < m_captureToCompounds.size( ) - 1; ++i1 ) {
3218 GRIN_captureToCompound const *GRIN_captureToCompound1 = m_captureToCompounds[i1];
3219
3220 int index = GRIN_captureToCompound1->sampleCaptureLevel( a_protare, a_energy, a_rng, true );
3221 if( index > -1 ) return( index );
3222 }
3223
3224 return( m_captureToCompounds.back( )->sampleCaptureLevel( a_protare, a_energy, a_rng, false ) );
3225}
3226
3227/* *********************************************************************************************************//**
3228 * This method adds a null product to *a_products*. When running in multi-group mode, a sampled reaction may be rejected if the threshold
3229 * 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
3230 * in MCAPM but is probably not needed for MCGIDI.
3231 *
3232 * @param a_protare [in] The Protare this Reaction belongs to.
3233 * @param a_projectileEnergy [in] The energy of the projectile.
3234 * @param a_input [in] Sample options requested by user.
3235 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
3236 * @param a_products [in] The object to add all sampled products to.
3237 ***********************************************************************************************************/
3238
3239template <typename RNG, typename PUSHBACK>
3240LUPI_HOST_DEVICE void MCGIDI::Reaction::sampleNullProducts( Protare const &a_protare, double a_projectileEnergy, Sampling::Input &a_input,
3241 RNG && a_rng, PUSHBACK && a_push_back, Sampling::ProductHandler &a_products ) {
3242
3243 a_input.m_GRIN_intermediateResidual = -1;
3244 a_input.setSampledType( Sampling::SampledType::uncorrelatedBody );
3245 a_input.m_dataInTargetFrame = false;
3246 a_input.m_frame = GIDI::Frame::lab;
3247 a_input.m_delayedNeutronIndex = -1;
3248 a_input.m_delayedNeutronDecayRate = 0.0;
3249
3250 a_input.m_energyOut1 = a_projectileEnergy;
3251 a_input.m_mu = 1.0;
3252 a_input.m_phi = 0.0;
3253
3254 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 );
3255}
3256
3257/* *********************************************************************************************************//**
3258 * Returns the weight for a project with energy *a_energy_in* to cause this reaction to emitted a particle of index
3259 * *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.
3260 *
3261 * @param a_index [in] The index of the particle to emit.
3262 * @param a_temperature [in] Specifies the temperature of the material.
3263 * @param a_energy_in [in] The energy of the incident particle.
3264 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
3265 * @param a_energy_out [in] The energy of the emitted outgoing particle.
3266 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
3267 * @param a_cumulative_weight [in] The cumulative multiplicity.
3268 * @param a_checkOrphanProducts [in] If true, associated orphan products are also sampled.
3269 *
3270 * @return The weight that the particle is emitted into mu *a_mu_lab*.
3271 ***********************************************************************************************************/
3272
3273template <typename RNG>
3274LUPI_HOST_DEVICE double MCGIDI::Reaction::angleBiasing( int a_index, double a_temperature, double a_energy_in, double a_mu_lab, double &a_energy_out,
3275 RNG && a_rng, double *a_cumulative_weight, bool a_checkOrphanProducts ) const {
3276
3277 double cumulative_weight1 = 0.0;
3278 if( a_cumulative_weight == nullptr ) a_cumulative_weight = &cumulative_weight1;
3279 double weight1 = 0.0;
3280
3281#ifdef MCGIDI_USE_OUTPUT_CHANNEL
3282 m_outputChannel->angleBiasing( this, a_index, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng, *a_cumulative_weight );
3283#else
3284
3285 for( auto productIter = m_products.begin( ); productIter != m_products.end( ); ++productIter ) {
3286 (*productIter)->angleBiasing( this, a_index, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng, *a_cumulative_weight );
3287 }
3288
3289 if( ( m_totalDelayedNeutronMultiplicity != nullptr ) && ( a_index == m_neutronIndex ) ) {
3290 for( std::size_t i1 = 0; i1 < (std::size_t) m_delayedNeutrons.size( ); ++i1 ) {
3291 DelayedNeutron const *delayedNeutron1 = m_delayedNeutrons[i1];
3292 Product const &product = delayedNeutron1->product( );
3293
3294 product.angleBiasing( this, a_index, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng, *a_cumulative_weight );
3295 }
3296 }
3297#endif
3298
3299 if( a_checkOrphanProducts ) {
3300 for( auto productIter = m_associatedOrphanProducts.begin( ); productIter != m_associatedOrphanProducts.end( ); ++productIter ) {
3301 (*productIter)->angleBiasing( this, a_index, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng,
3302 *a_cumulative_weight );
3303 }
3304 }
3305
3306 return( weight1 );
3307}
3308
3309/* *********************************************************************************************************//**
3310 * Returns the weight for a project with energy *a_energy_in* to cause this reaction to emitted a particle of intid
3311 * *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.
3312 *
3313 * @param a_intid [in] The intid of the particle to emit.
3314 * @param a_temperature [in] Specifies the temperature of the material.
3315 * @param a_energy_in [in] The energy of the incident particle.
3316 * @param a_mu_lab [in] The desired mu in the lab frame for the emitted particle.
3317 * @param a_energy_out [in] The energy of the emitted outgoing particle.
3318 * @param a_rng [in] The random number generator function that returns a double in the range [0, 1.0).
3319 * @param a_cumulative_weight [in] The cumulative multiplicity.
3320 * @param a_checkOrphanProducts [in] If true, associated orphan products are also sampled.
3321 *
3322 * @return The weight that the particle is emitted into mu *a_mu_lab*.
3323 ***********************************************************************************************************/
3324
3325template <typename RNG>
3326LUPI_HOST_DEVICE double MCGIDI::Reaction::angleBiasingViaIntid( int a_intid, double a_temperature, double a_energy_in, double a_mu_lab, double &a_energy_out,
3327 RNG && a_rng, double *a_cumulative_weight, bool a_checkOrphanProducts ) const {
3328
3329 double cumulative_weight1 = 0.0;
3330 if( a_cumulative_weight == nullptr ) a_cumulative_weight = &cumulative_weight1;
3331 double weight1 = 0.0;
3332
3333#ifdef MCGIDI_USE_OUTPUT_CHANNEL
3334 m_outputChannel->angleBiasingViaIntid( this, a_intid, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng, *a_cumulative_weight );
3335#else
3336
3337 for( auto productIter = m_products.begin( ); productIter != m_products.end( ); ++productIter ) {
3338 (*productIter)->angleBiasingViaIntid( this, a_intid, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng, *a_cumulative_weight );
3339 }
3340
3341 if( ( m_totalDelayedNeutronMultiplicity != nullptr ) && ( a_intid == PoPI::Intids::neutron ) ) {
3342 for( std::size_t i1 = 0; i1 < (std::size_t) m_delayedNeutrons.size( ); ++i1 ) {
3343 DelayedNeutron const *delayedNeutron1 = m_delayedNeutrons[i1];
3344 Product const &product = delayedNeutron1->product( );
3345
3346 product.angleBiasingViaIntid( this, a_intid, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng, *a_cumulative_weight );
3347 }
3348 }
3349#endif
3350
3351 if( a_checkOrphanProducts ) {
3352 for( auto productIter = m_associatedOrphanProducts.begin( ); productIter != m_associatedOrphanProducts.end( ); ++productIter ) {
3353 (*productIter)->angleBiasingViaIntid( this, a_intid, a_temperature, a_energy_in, a_mu_lab, weight1, a_energy_out, a_rng,
3354 *a_cumulative_weight );
3355 }
3356 }
3357
3358 return( weight1 );
3359}
3360
3361// From file: MCGIDI_sampling.cpp
3362
3363template <typename RNG, typename PUSHBACK>
3364LUPI_HOST_DEVICE void MCGIDI::Sampling::ProductHandler::add( double a_projectileEnergy, int a_productIntid, int a_productIndex, int a_userProductIndex,
3365 double a_productMass, Input &a_input, RNG && a_rng, PUSHBACK && a_push_back, bool a_isPhoton ) {
3366
3367 Product product;
3368
3369 if( a_isPhoton && ( a_input.m_sampledType != SampledType::unspecified ) ) a_input.m_sampledType = SampledType::photon;
3370
3371 product.m_sampledType = a_input.m_sampledType;
3372 product.m_isVelocity = a_input.wantVelocity( );
3373 product.m_productIntid = a_productIntid;
3374 product.m_productIndex = a_productIndex;
3375 product.m_userProductIndex = a_userProductIndex;
3376 product.m_numberOfDBRC_rejections = a_input.m_numberOfDBRC_rejections;
3377 product.m_productMass = a_productMass;
3378
3379 product.m_delayedNeutronIndex = a_input.m_delayedNeutronIndex;
3380 product.m_delayedNeutronDecayRate = a_input.m_delayedNeutronDecayRate;
3381 product.m_birthTimeSec = 0.;
3382 if( product.m_delayedNeutronDecayRate > 0. ) {
3383 product.m_birthTimeSec = -log( a_rng( ) ) / product.m_delayedNeutronDecayRate;
3384 }
3385
3386 if( a_input.m_sampledType == SampledType::unspecified ) {
3387 product.m_kineticEnergy = 0.0;
3388 product.m_px_vx = 0.0;
3389 product.m_py_vy = 0.0;
3390 product.m_pz_vz = 0.0; }
3391 else if( a_input.m_sampledType == SampledType::uncorrelatedBody ) {
3392 if( a_input.m_frame == GIDI::Frame::centerOfMass ) {
3393 a_input.m_frame = GIDI::Frame::lab;
3394
3395 double massRatio = a_input.m_projectileMass + a_input.m_targetMass;
3396 massRatio = a_input.m_projectileMass * a_productMass / ( massRatio * massRatio );
3397 double modifiedProjectileEnergy = massRatio * a_projectileEnergy;
3398
3399 double sqrtModifiedProjectileEnergy = sqrt( modifiedProjectileEnergy );
3400 double sqrtEnergyOut_com = a_input.m_mu * sqrt( a_input.m_energyOut1 );
3401
3402 a_input.m_energyOut1 += modifiedProjectileEnergy + 2. * sqrtModifiedProjectileEnergy * sqrtEnergyOut_com;
3403 if( a_input.m_energyOut1 != 0 ) a_input.m_mu = ( sqrtModifiedProjectileEnergy + sqrtEnergyOut_com ) / sqrt( a_input.m_energyOut1 );
3404 }
3405
3406 product.m_kineticEnergy = a_input.m_energyOut1;
3407
3408 double p_v = sqrt( a_input.m_energyOut1 * ( a_input.m_energyOut1 + 2. * a_productMass ) );
3409 if( product.m_isVelocity ) p_v *= MCGIDI_speedOfLight_cm_sec / ( a_input.m_energyOut1 + a_productMass );
3410
3411 product.m_pz_vz = p_v * a_input.m_mu;
3412 p_v *= sqrt( 1. - a_input.m_mu * a_input.m_mu );
3413 product.m_px_vx = p_v * sin( a_input.m_phi );
3414 product.m_py_vy = p_v * cos( a_input.m_phi ); }
3415 else if( a_input.m_sampledType == SampledType::firstTwoBody ) {
3416 product.m_kineticEnergy = a_input.m_energyOut1;
3417 product.m_px_vx = a_input.m_px_vx1;
3418 product.m_py_vy = a_input.m_py_vy1;
3419 product.m_pz_vz = a_input.m_pz_vz1;
3420 a_input.m_sampledType = SampledType::secondTwoBody; }
3421 else if( a_input.m_sampledType == SampledType::secondTwoBody ) {
3422 product.m_kineticEnergy = a_input.m_energyOut2;
3423 product.m_px_vx = a_input.m_px_vx2;
3424 product.m_py_vy = a_input.m_py_vy2;
3425 product.m_pz_vz = a_input.m_pz_vz2; }
3426 else if( a_input.m_sampledType == SampledType::photon ) {
3427 product.m_kineticEnergy = a_input.m_energyOut1;
3428
3429 double pz_vz_factor = a_input.m_energyOut1;
3430 if( product.m_isVelocity ) pz_vz_factor = MCGIDI_speedOfLight_cm_sec;
3431 product.m_pz_vz = a_input.m_mu * pz_vz_factor;
3432
3433 double v_perp = sqrt( 1.0 - a_input.m_mu * a_input.m_mu ) * pz_vz_factor;
3434 product.m_px_vx = cos( a_input.m_phi ) * v_perp;
3435 product.m_py_vy = sin( a_input.m_phi ) * v_perp; }
3436 else {
3437 product.m_kineticEnergy = a_input.m_energyOut2;
3438 product.m_px_vx = a_input.m_px_vx2;
3439 product.m_py_vy = a_input.m_py_vy2;
3440 product.m_pz_vz = a_input.m_pz_vz2;
3441 }
3442
3443 if( a_input.m_dataInTargetFrame && ( a_input.m_sampledType != SampledType::photon ) ) upScatterModelABoostParticle( a_input, a_rng, product );
3444
3445 a_push_back( product );
3446}
3447
3448#endif // End of MCGIDI_headerSource_hpp_included
epsilon
G4double epsilon(G4double density, G4double temperature)
Definition G4DNAElectronHoleRecombination.cc:98
LUPI_HOST_DEVICE
#define LUPI_HOST_DEVICE
Definition LUPI_declareMacro.hpp:83
LUPI_THROW
#define LUPI_THROW(arg)
Definition LUPI_declareMacro.hpp:84
LUPI_maybeUnused
#define LUPI_maybeUnused
Definition LUPI_defines.hpp:16
MCGIDI_speedOfLight_cm_sec
#define MCGIDI_speedOfLight_cm_sec
Definition MCGIDI.hpp:68
MCGIDI_particleBeta
#define MCGIDI_particleBeta(a_mass_unitOfEnergy, a_kineticEnergy)
Definition MCGIDI.hpp:71
MCGIDI_classicalElectronRadius
#define MCGIDI_classicalElectronRadius
Definition MCGIDI.hpp:69
MCGIDI_nullReaction
#define MCGIDI_nullReaction
Definition MCGIDI.hpp:64
MCGIDI_PRINTF
#define MCGIDI_PRINTF
Definition MCGIDI.hpp:40
crossSectionSumError
#define crossSectionSumError
Definition MCGIDI_headerSource.hpp:19
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
kinetics_COMKineticEnergy2LabEnergyAndMomentum
LUPI_HOST_DEVICE void kinetics_COMKineticEnergy2LabEnergyAndMomentum(double a_beta, double a_kinetic_com, double a_m3cc, double a_m4cc, MCGIDI::Sampling::Input &a_input)
Definition MCGIDI_headerSource.hpp:283
sampleBetaFromMaxwellian
LUPI_HOST_DEVICE double sampleBetaFromMaxwellian(RNG &&a_rng)
Definition MCGIDI_headerSource.hpp:2333
PoPI_electronMass_MeV_c2
#define PoPI_electronMass_MeV_c2
Definition PoPI.hpp:31
C1
#define C1
C3
#define C3
M_PI
#define M_PI
Definition SbMath.h:33
MCGIDI::DelayedNeutron::rate
LUPI_HOST_DEVICE double rate() const
Definition MCGIDI.hpp:1259
MCGIDI::DelayedNeutron::delayedNeutronIndex
LUPI_HOST_DEVICE int delayedNeutronIndex() const
Definition MCGIDI.hpp:1258
MCGIDI::DelayedNeutron::product
LUPI_HOST_DEVICE Product const & product() const
Definition MCGIDI.hpp:1260
MCGIDI::Distributions::AngularEnergyMC::angleBiasing
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
MCGIDI::Distributions::AngularEnergyMC::sample
LUPI_HOST_DEVICE void sample(double a_X, Sampling::Input &a_input, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:832
MCGIDI::Distributions::AngularTwoBody::sample
LUPI_HOST_DEVICE void sample(double a_X, Sampling::Input &a_input, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:366
MCGIDI::Distributions::AngularTwoBody::angleBiasing
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
MCGIDI::Distributions::Branching3d::sample
LUPI_HOST_DEVICE void sample(double a_X, Sampling::Input &a_input, RNG &&a_rng) const
MCGIDI::Distributions::Branching3d::angleBiasing
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
MCGIDI::Distributions::CoherentElasticTNSL::angleBiasing
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
MCGIDI::Distributions::CoherentElasticTNSL::sample
LUPI_HOST_DEVICE void sample(double a_energy, Sampling::Input &a_input, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:1423
MCGIDI::Distributions::CoherentPhotoAtomicScattering::sample
LUPI_HOST_DEVICE void sample(double a_X, Sampling::Input &a_input, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:975
MCGIDI::Distributions::CoherentPhotoAtomicScattering::angleBiasing
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
MCGIDI::Distributions::CoherentPhotoAtomicScattering::evaluateFormFactor
LUPI_HOST_DEVICE double evaluateFormFactor(double a_energyIn, double a_mu) const
Definition MCGIDI_distributions.cc:702
MCGIDI::Distributions::Distribution::angleBiasing
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
MCGIDI::Distributions::Distribution::targetMass
LUPI_HOST_DEVICE double targetMass() const
Definition MCGIDI_distributions.hpp:47
MCGIDI::Distributions::Distribution::sample
LUPI_HOST_DEVICE MCGIDI_VIRTUAL_FUNCTION void sample(double a_X, Sampling::Input &a_input, RNG &&a_rng) const MCGIDI_TRUE_VIRTUAL
MCGIDI::Distributions::Distribution::productFrame
LUPI_HOST_DEVICE GIDI::Frame productFrame() const
Definition MCGIDI_distributions.hpp:44
MCGIDI::Distributions::Distribution::productMass
LUPI_HOST_DEVICE double productMass() const
Definition MCGIDI_distributions.hpp:48
MCGIDI::Distributions::Distribution::projectileMass
LUPI_HOST_DEVICE double projectileMass() const
Definition MCGIDI_distributions.hpp:46
MCGIDI::Distributions::EnergyAngularMC::sample
LUPI_HOST_DEVICE void sample(double a_X, Sampling::Input &a_input, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:747
MCGIDI::Distributions::EnergyAngularMC::angleBiasing
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
MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::sample
LUPI_HOST_DEVICE void sample(double a_X, Sampling::Input &a_input, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:1124
MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::energyRatio
LUPI_HOST_DEVICE double energyRatio(double a_energyIn, double a_mu) const
Definition MCGIDI_distributions.cc:1038
MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::angleBiasing
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
MCGIDI::Distributions::IncoherentBoundToFreePhotoAtomicScattering::evaluateOccupationNumber
LUPI_HOST_DEVICE double evaluateOccupationNumber(double a_X, double a_mu) const
Definition MCGIDI_distributions.cc:1074
MCGIDI::Distributions::IncoherentElasticTNSL::angleBiasing
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
MCGIDI::Distributions::IncoherentElasticTNSL::sample
LUPI_HOST_DEVICE void sample(double a_energy, Sampling::Input &a_input, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:1494
MCGIDI::Distributions::IncoherentPhotoAtomicScatteringElectron::angleBiasing
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
MCGIDI::Distributions::IncoherentPhotoAtomicScatteringElectron::sample
LUPI_HOST_DEVICE void sample(double a_energy, Sampling::Input &a_input, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:1328
MCGIDI::Distributions::IncoherentPhotoAtomicScattering::sample
LUPI_HOST_DEVICE void sample(double a_X, Sampling::Input &a_input, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:1092
MCGIDI::Distributions::IncoherentPhotoAtomicScattering::angleBiasing
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
MCGIDI::Distributions::IncoherentPhotoAtomicScattering::evaluateScatteringFactor
LUPI_HOST_DEVICE double evaluateScatteringFactor(double a_X) const
Definition MCGIDI_distributions.cc:941
MCGIDI::Distributions::KalbachMann::angleBiasing
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
MCGIDI::Distributions::KalbachMann::sample
LUPI_HOST_DEVICE void sample(double a_X, Sampling::Input &a_input, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:877
MCGIDI::Distributions::PairProductionGamma::angleBiasing
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
MCGIDI::Distributions::PairProductionGamma::sample
LUPI_HOST_DEVICE void sample(double a_X, Sampling::Input &a_input, RNG &&a_rng) const
MCGIDI::Distributions::Uncorrelated::angleBiasing
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
MCGIDI::Distributions::Uncorrelated::sample
LUPI_HOST_DEVICE void sample(double a_X, Sampling::Input &a_input, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:630
MCGIDI::Distributions::Unspecified::angleBiasing
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
MCGIDI::Distributions::Unspecified::sample
LUPI_HOST_DEVICE void sample(double a_X, Sampling::Input &a_input, RNG &&a_rng) const
MCGIDI::Functions::Function1d::type
LUPI_HOST_DEVICE Function1dType type() const
Definition MCGIDI_functions.hpp:72
MCGIDI::Functions::Function1d::sampleBoundingInteger
LUPI_HOST_DEVICE MCGIDI_VIRTUAL_FUNCTION int sampleBoundingInteger(double a_x1, RNG &&a_rng) const
MCGIDI::Functions::Function1d::evaluate
LUPI_HOST_DEVICE MCGIDI_VIRTUAL_FUNCTION double evaluate(double a_x1) const MCGIDI_TRUE_VIRTUAL
Definition MCGIDI_functions.cc:197
MCGIDI::Functions::FunctionBase::interpolation
LUPI_HOST_DEVICE Interpolation interpolation() const
Definition MCGIDI_functions.hpp:50
MCGIDI::Functions::TerrellFissionNeutronMultiplicityModel::sampleBoundingInteger
LUPI_HOST_DEVICE int sampleBoundingInteger(double a_energy, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:1617
MCGIDI::GRIN_captureLevelProbability
Definition MCGIDI.hpp:1121
MCGIDI::GRIN_captureLevelProbability::sampleCaptureLevel
LUPI_HOST_DEVICE int sampleCaptureLevel(ProtareSingle const *a_protare, double a_energy, RNG &&a_rng)
Definition MCGIDI_headerSource.hpp:3208
MCGIDI::GRIN_captureToCompound
Definition MCGIDI.hpp:1090
MCGIDI::GRIN_captureToCompound::index
LUPI_HOST_DEVICE std::size_t index() const
Definition MCGIDI.hpp:1101
MCGIDI::GRIN_captureToCompound::sampleCaptureLevel
LUPI_HOST_DEVICE int sampleCaptureLevel(ProtareSingle const *a_protare, double a_energy, RNG &&a_rng, bool a_checkEnergy) const
Definition MCGIDI_headerSource.hpp:3180
MCGIDI::GRIN_capture::sampleProducts
LUPI_HOST_DEVICE bool sampleProducts(ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input, RNG &&a_rng, PUSHBACK &&a_push_back, Sampling::ProductHandler &a_products) const
Definition MCGIDI_headerSource.hpp:3069
MCGIDI::GRIN_inelasticForEnergy
Definition MCGIDI.hpp:1032
MCGIDI::GRIN_inelasticForEnergy::sampleLevelIndex
LUPI_HOST_DEVICE int sampleLevelIndex(double a_projectileEnergy, double a_random) const
Definition MCGIDI_GRIN.cc:152
MCGIDI::GRIN_inelastic::sampleProducts
LUPI_HOST_DEVICE bool sampleProducts(ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input, RNG &&a_rng, PUSHBACK &&a_push_back, Sampling::ProductHandler &a_products) const
Definition MCGIDI_headerSource.hpp:3119
MCGIDI::HeatedCrossSectionContinuousEnergy::reactionCrossSection2
LUPI_HOST_DEVICE double reactionCrossSection2(std::size_t a_reactionIndex, URR_protareInfos const &a_URR_protareInfos, int a_URR_index, double a_energy, std::size_t a_energyIndex, double a_energyFraction, bool a_sampling=false) const
Definition MCGIDI_heatedCrossSections.cc:717
MCGIDI::HeatedCrossSectionContinuousEnergy::evaluationInfo
LUPI_HOST_DEVICE std::size_t evaluationInfo(std::size_t a_hashIndex, double a_energy, double *a_energyFraction) const
Definition MCGIDI_heatedCrossSections.cc:626
MCGIDI::HeatedCrossSectionMultiGroup::reactionCrossSection
LUPI_HOST_DEVICE double reactionCrossSection(std::size_t a_reactionIndex, std::size_t a_hashIndex, bool a_sampling=false) const
Definition MCGIDI.hpp:858
MCGIDI::HeatedCrossSectionsContinuousEnergy::sampleReaction
LUPI_HOST_DEVICE std::size_t 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
Definition MCGIDI_headerSource.hpp:2187
MCGIDI::HeatedCrossSectionsMultiGroup::sampleReaction
LUPI_HOST_DEVICE std::size_t sampleReaction(std::size_t a_hashIndex, double a_temperature, double a_energy_in, double a_crossSection, RNG &&rng) const
Definition MCGIDI_headerSource.hpp:2267
MCGIDI::NuclideGammaBranchInfo::probability
LUPI_HOST_DEVICE double probability() const
Definition MCGIDI.hpp:952
MCGIDI::NuclideGammaBranchInfo::residualStateIndex
LUPI_HOST_DEVICE int residualStateIndex() const
Definition MCGIDI.hpp:955
MCGIDI::NuclideGammaBranchInfo::gammaEnergy
LUPI_HOST_DEVICE double gammaEnergy() const
Definition MCGIDI.hpp:954
MCGIDI::NuclideGammaBranchInfo::photonEmissionProbability
LUPI_HOST_DEVICE double photonEmissionProbability() const
Definition MCGIDI.hpp:953
MCGIDI::NuclideGammaBranchStateInfo::branchIndices
LUPI_HOST_DEVICE Vector< std::size_t > const & branchIndices() const
Definition MCGIDI.hpp:993
MCGIDI::NuclideGammaBranchStateInfo::nuclearLevelEnergy
LUPI_HOST_DEVICE double nuclearLevelEnergy() const
Definition MCGIDI.hpp:988
MCGIDI::NuclideGammaBranchStateInfo::intid
LUPI_HOST_DEVICE int intid() const
Definition MCGIDI.hpp:987
MCGIDI::NuclideGammaBranchStateInfo::nuclearLevelEnergyWidth
LUPI_HOST_DEVICE double nuclearLevelEnergyWidth() const
Definition MCGIDI.hpp:989
MCGIDI::OutputChannel::angleBiasing
LUPI_HOST_DEVICE void angleBiasing(Reaction const *a_reaction, int a_pid, double a_temperature, double a_energy_in, double a_mu_lab, double &a_probability, double &a_energy_out, RNG &&a_rng, double &a_cumulative_weight) const
Definition MCGIDI_headerSource.hpp:2469
MCGIDI::OutputChannel::angleBiasingViaIntid
LUPI_HOST_DEVICE void angleBiasingViaIntid(Reaction const *a_reaction, int a_intid, double a_temperature, double a_energy_in, double a_mu_lab, double &a_probability, double &a_energy_out, RNG &&a_rng, double &a_cumulative_weight) const
Definition MCGIDI_headerSource.hpp:2501
MCGIDI::OutputChannel::sampleProducts
LUPI_HOST_DEVICE void sampleProducts(ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input, RNG &&a_rng, PUSHBACK &&a_push_back, Sampling::ProductHandler &a_products) const
Definition MCGIDI_headerSource.hpp:2406
MCGIDI::OutputChannel::delayedNeutron
LUPI_HOST_DEVICE DelayedNeutron const * delayedNeutron(std::size_t a_index) const
Definition MCGIDI.hpp:1306
MCGIDI::Probabilities::Evaporation2d::sample
LUPI_HOST_DEVICE double sample(double a_x2, double a_rngValue, RNG &&a_rng) const
MCGIDI::Probabilities::GeneralEvaporation2d::sample
LUPI_HOST_DEVICE double sample(double a_x2, double a_rngValue, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:2040
MCGIDI::Probabilities::NBodyPhaseSpace2d::sample
LUPI_HOST_DEVICE double sample(double a_x2, double a_rngValue, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:2003
MCGIDI::Probabilities::ProbabilityBase1d::sample
LUPI_HOST_DEVICE MCGIDI_VIRTUAL_FUNCTION double sample(double a_rngValue, RNG &&a_rng) const MCGIDI_TRUE_VIRTUAL
MCGIDI::Probabilities::ProbabilityBase2d_d1
Definition MCGIDI_functions.hpp:420
MCGIDI::Probabilities::ProbabilityBase2d_d1::sample
LUPI_HOST_DEVICE double sample(double a_x2, double a_rngValue, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:1732
MCGIDI::Probabilities::ProbabilityBase2d_d1::sample2dOf3d
LUPI_HOST_DEVICE double sample2dOf3d(double a_x2, double a_rngValue, RNG &&a_rng, double *a_x1_1, double *a_x1_2) const
Definition MCGIDI_headerSource.hpp:1771
MCGIDI::Probabilities::ProbabilityBase2d_d2
Definition MCGIDI_functions.hpp:442
MCGIDI::Probabilities::ProbabilityBase2d_d2::sample2dOf3d
LUPI_HOST_DEVICE double sample2dOf3d(double a_x2, double a_rngValue, RNG &&a_rng, double *a_x1_1, double *a_x1_2) const
Definition MCGIDI_headerSource.hpp:1853
MCGIDI::Probabilities::ProbabilityBase2d_d2::sample
LUPI_HOST_DEVICE double sample(double a_x2, double a_rngValue, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:1798
MCGIDI::Probabilities::ProbabilityBase2d::sample
LUPI_HOST_DEVICE MCGIDI_VIRTUAL_FUNCTION double sample(double a_x2, double a_rngValue, RNG &&a_rng) const MCGIDI_TRUE_VIRTUAL
MCGIDI::Probabilities::ProbabilityBase2d::type
LUPI_HOST_DEVICE ProbabilityBase2dType type() const
Definition MCGIDI_functions.hpp:406
MCGIDI::Probabilities::ProbabilityBase3d::sample
LUPI_HOST_DEVICE MCGIDI_VIRTUAL_FUNCTION double sample(double a_x3, double a_x2_1, double a_x2_2, double a_rngValue, RNG &&a_rng) const MCGIDI_TRUE_VIRTUAL
MCGIDI::Probabilities::Recoil2d::sample
LUPI_HOST_DEVICE double sample(double a_x2, double a_rngValue, RNG &&a_rng) const
MCGIDI::Probabilities::Regions2d::sample
LUPI_HOST_DEVICE double sample(double a_x2, double a_rngValue, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:1968
MCGIDI::Probabilities::SimpleMaxwellianFission2d::sample
LUPI_HOST_DEVICE double sample(double a_x2, double a_rngValue, RNG &&a_rng) const
MCGIDI::Probabilities::Watt2d::sample
LUPI_HOST_DEVICE double sample(double a_x2, double a_rngValue, RNG &&a_rng) const
MCGIDI::Probabilities::WeightedFunctionals2d::sample
LUPI_HOST_DEVICE double sample(double a_x2, double a_rngValue, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:2096
MCGIDI::Probabilities::XYs2d::sample2dOf3d
LUPI_HOST_DEVICE double sample2dOf3d(double a_x2, double a_rngValue, RNG &&a_rng, double *a_x1_1, double *a_x1_2) const
Definition MCGIDI_headerSource.hpp:1923
MCGIDI::Probabilities::XYs2d::sample
LUPI_HOST_DEVICE double sample(double a_x2, double a_rngValue, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:1870
MCGIDI::Probabilities::XYs3d::sample
LUPI_HOST_DEVICE double sample(double a_x3, double a_x2_1, double a_x2_2, double a_rngValue, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:2128
MCGIDI::Probabilities::Xs_pdf_cdf1d
Definition MCGIDI_functions.hpp:373
MCGIDI::Probabilities::Xs_pdf_cdf1d::sample
LUPI_HOST_DEVICE double sample(double a_rngValue, RNG &&a_rng) const
MCGIDI::Product::userParticleIndex
LUPI_HOST_DEVICE int userParticleIndex() const
Definition MCGIDI.hpp:1202
MCGIDI::Product::sampleProducts
LUPI_HOST_DEVICE void sampleProducts(ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input, RNG &&a_rng, PUSHBACK &&a_push_back, Sampling::ProductHandler &a_products) const
Definition MCGIDI_headerSource.hpp:2531
MCGIDI::Product::index
LUPI_HOST_DEVICE int index() const
Definition MCGIDI.hpp:1201
MCGIDI::Product::angleBiasing
LUPI_HOST_DEVICE void angleBiasing(Reaction const *a_reaction, int a_pid, double a_temperature, double a_energy_in, double a_mu_lab, double &a_probability, double &a_energy_out, RNG &&a_rng, double &a_cumulative_weight) const
Definition MCGIDI_headerSource.hpp:2602
MCGIDI::Product::sampleFinalState
LUPI_HOST_DEVICE void sampleFinalState(ProtareSingle const *a_protare, double a_projectileEnergy, Sampling::Input &a_input, RNG &&a_rng, PUSHBACK &&a_push_back, Sampling::ProductHandler &a_products) const
Definition MCGIDI_headerSource.hpp:2574
MCGIDI::Product::intid
LUPI_HOST_DEVICE int intid() const
Definition MCGIDI.hpp:1200
MCGIDI::Product::mass
LUPI_HOST_DEVICE double mass() const
Definition MCGIDI.hpp:1208
MCGIDI::Product::angleBiasingViaIntid
LUPI_HOST_DEVICE void angleBiasingViaIntid(Reaction const *a_reaction, int a_intid, double a_temperature, double a_energy_in, double a_mu_lab, double &a_probability, double &a_energy_out, RNG &&a_rng, double &a_cumulative_weight) const
Definition MCGIDI_headerSource.hpp:2650
MCGIDI::ProtareComposite::sampleReaction
LUPI_HOST_DEVICE std::size_t sampleReaction(Sampling::Input &a_input, URR_protareInfos const &a_URR_protareInfos, std::size_t a_hashIndex, double a_crossSection, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:2888
MCGIDI::ProtareSingle::sampleReaction
LUPI_HOST_DEVICE std::size_t sampleReaction(Sampling::Input &a_input, URR_protareInfos const &a_URR_protareInfos, std::size_t a_hashIndex, double a_crossSection, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:2843
MCGIDI::ProtareSingle::sampleBranchingGammas
LUPI_HOST_DEVICE void sampleBranchingGammas(Sampling::Input &a_input, double a_projectileEnergy, int a_initialStateIndex, RNG &&a_rng, PUSHBACK &&push_back, Sampling::ProductHandler &a_products) const
Definition MCGIDI_headerSource.hpp:2718
MCGIDI::ProtareSingle::inURR
LUPI_HOST_DEVICE bool inURR(double a_energy) const
Definition MCGIDI_protare.cc:1799
MCGIDI::ProtareSingle::crossSection
LUPI_HOST_DEVICE double crossSection(URR_protareInfos const &a_URR_protareInfos, std::size_t a_hashIndex, double a_temperature, double a_energy, bool a_sampling=false) const
Definition MCGIDI_protare.cc:1818
MCGIDI::ProtareSingle::reactionCrossSection
LUPI_HOST_DEVICE double reactionCrossSection(std::size_t a_reactionIndex, URR_protareInfos const &a_URR_protareInfos, std::size_t a_hashIndex, double a_temperature, double a_energy, bool a_sampling=false) const
Definition MCGIDI_protare.cc:1857
MCGIDI::ProtareSingle::sampleTargetBetaForUpscatterModelA
LUPI_HOST_DEVICE bool sampleTargetBetaForUpscatterModelA(Sampling::Input &a_input, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:2772
MCGIDI::ProtareSingle::nuclideGammaBranchStateInfos
LUPI_HOST_DEVICE const Vector< NuclideGammaBranchStateInfo * > & nuclideGammaBranchStateInfos() const
Definition MCGIDI.hpp:1658
MCGIDI::ProtareSingle::URR_index
LUPI_HOST_DEVICE int URR_index() const
Definition MCGIDI.hpp:1703
MCGIDI::ProtareSingle::upscatterModelASupported
LUPI_HOST_DEVICE bool upscatterModelASupported() const
Definition MCGIDI.hpp:1735
MCGIDI::ProtareTNSL::sampleReaction
LUPI_HOST_DEVICE std::size_t sampleReaction(Sampling::Input &a_input, URR_protareInfos const &a_URR_protareInfos, std::size_t a_hashIndex, double a_crossSection, RNG &&a_rng) const
Definition MCGIDI_headerSource.hpp:2929
MCGIDI::Protare::projectileIntid
LUPI_HOST_DEVICE int projectileIntid() const
Definition MCGIDI.hpp:1517
MCGIDI::Protare::protareType
LUPI_HOST_DEVICE ProtareType protareType() const
Definition MCGIDI.hpp:1514
MCGIDI::Protare::ProtareTNSL
friend ProtareTNSL
Definition MCGIDI.hpp:1601
MCGIDI::Protare::projectileIndex
LUPI_HOST_DEVICE int projectileIndex() const
Definition MCGIDI.hpp:1518
MCGIDI::Protare::userPhotonIndex
LUPI_HOST_DEVICE int userPhotonIndex() const
Definition MCGIDI.hpp:1531
MCGIDI::Protare::ProtareComposite
friend ProtareComposite
Definition MCGIDI.hpp:1600
MCGIDI::Protare::ProtareSingle
friend ProtareSingle
Definition MCGIDI.hpp:1599
MCGIDI::Protare::projectileUserIndex
LUPI_HOST_DEVICE int projectileUserIndex() const
Definition MCGIDI.hpp:1519
MCGIDI::Protare::photonIndex
LUPI_HOST_DEVICE int photonIndex() const
Definition MCGIDI.hpp:1530
MCGIDI::Protare::numberOfProtares
MCGIDI_VIRTUAL_FUNCTION LUPI_HOST_DEVICE std::size_t numberOfProtares() const MCGIDI_TRUE_VIRTUAL
Definition MCGIDI_protare.cc:465
MCGIDI::Protare::sampleReaction
MCGIDI_VIRTUAL_FUNCTION LUPI_HOST_DEVICE std::size_t sampleReaction(Sampling::Input &a_input, URR_protareInfos const &a_URR_protareInfos, std::size_t a_hashIndex, double a_crossSection, RNG &&a_rng) const MCGIDI_TRUE_VIRTUAL
MCGIDI::Protare::projectileMass
LUPI_HOST_DEVICE double projectileMass() const
Definition MCGIDI.hpp:1520
MCGIDI::Protare::protare
MCGIDI_VIRTUAL_FUNCTION LUPI_HOST_DEVICE ProtareSingle const * protare(std::size_t a_index) const MCGIDI_TRUE_VIRTUAL
Definition MCGIDI_protare.cc:492
MCGIDI::Reaction::sampleNullProducts
static LUPI_HOST_DEVICE void sampleNullProducts(Protare const &a_protare, double a_projectileEnergy, Sampling::Input &a_input, RNG &&a_rng, PUSHBACK &&a_push_back, Sampling::ProductHandler &a_products)
Definition MCGIDI_headerSource.hpp:3240
MCGIDI::Reaction::angleBiasing
LUPI_HOST_DEVICE double angleBiasing(int a_pid, double a_temperature, double a_energy_in, double a_mu_lab, double &a_energy_out, RNG &&a_rng, double *a_cumulative_weight=nullptr, bool a_checkOrphanProducts=true) const
Definition MCGIDI_headerSource.hpp:3274
MCGIDI::Reaction::reactionIndex
LUPI_HOST_DEVICE std::size_t reactionIndex() const
Definition MCGIDI.hpp:1400
MCGIDI::Reaction::product
LUPI_HOST_DEVICE Product const * product(std::size_t a_index) const
Definition MCGIDI.hpp:1422
MCGIDI::Reaction::protareSingle
LUPI_HOST_DEVICE ProtareSingle const * protareSingle() const
Definition MCGIDI.hpp:1399
MCGIDI::Reaction::sampleProducts
LUPI_HOST_DEVICE void sampleProducts(Protare const *a_protare, Sampling::Input &a_input, RNG &&a_rng, PUSHBACK &&a_push_back, Sampling::ProductHandler &a_products, bool a_checkOrphanProducts=true) const
Definition MCGIDI_headerSource.hpp:2966
MCGIDI::Reaction::angleBiasingViaIntid
LUPI_HOST_DEVICE double 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=nullptr, bool a_checkOrphanProducts=true) const
Definition MCGIDI_headerSource.hpp:3326
MCGIDI::Sampling::Input::wantVelocity
LUPI_HOST_DEVICE bool wantVelocity() const
Definition MCGIDI_sampling.hpp:200
MCGIDI::Sampling::Input::temperature
LUPI_HOST_DEVICE double temperature() const
Definition MCGIDI_sampling.hpp:201
MCGIDI::Sampling::Input::modelEnergy
LUPI_HOST_DEVICE double modelEnergy() const
Definition MCGIDI_sampling.hpp:207
MCGIDI::Sampling::Input::m_upscatterModel
Upscatter::Model m_upscatterModel
Definition MCGIDI_sampling.hpp:162
MCGIDI::Sampling::Input::setSampledType
LUPI_HOST_DEVICE void setSampledType(SampledType a_sampledType)
Definition MCGIDI_sampling.hpp:210
MCGIDI::Sampling::Input::energy
LUPI_HOST_DEVICE double energy() const
Definition MCGIDI_sampling.hpp:202
MCGIDI::Sampling::ProductHandler::add
LUPI_HOST_DEVICE void add(double a_projectileEnergy, int a_productIntid, int a_productIndex, int a_userProductIndex, double a_productMass, Input &a_input, RNG &&a_rng, PUSHBACK &&push_back, bool isPhoton)
Definition MCGIDI_headerSource.hpp:3364
MCGIDI::URR_protareInfo
Definition MCGIDI.hpp:420
MCGIDI::URR_protareInfo::m_inURR
bool m_inURR
Definition MCGIDI.hpp:423
MCGIDI::URR_protareInfo::m_rng_Value
double m_rng_Value
Definition MCGIDI.hpp:424
MCGIDI::URR_protareInfo::inURR
LUPI_HOST_DEVICE bool inURR() const
Definition MCGIDI.hpp:441
MCGIDI::URR_protareInfos
Definition MCGIDI.hpp:451
MCGIDI::URR_protareInfos::updateProtare
LUPI_HOST_DEVICE void updateProtare(MCGIDI::Protare const *a_protare, double a_energy, RNG &&a_rng)
Definition MCGIDI_headerSource.hpp:33
GIDI::Frame::centerOfMass
@ centerOfMass
Definition GIDI.hpp:146
LUPI::Misc::doubleToString3
std::string doubleToString3(char const *a_format, double a_value, bool a_reduceBits=false)
Definition LUPI_misc.cc:327
MCGIDI::upScatterModelABoostParticle
LUPI_HOST_DEVICE void upScatterModelABoostParticle(Sampling::Input &a_input, RNG &&a_rng, Sampling::Product &a_product)
Definition MCGIDI_headerSource.hpp:2358
MCGIDI::particleKineticEnergy
LUPI_HOST_DEVICE double particleKineticEnergy(double a_mass_unitOfEnergy, double a_particleBeta)
Definition MCGIDI_misc.cc:127
MCGIDI::binarySearchVector
LUPI_HOST_DEVICE int binarySearchVector(double a_x, Vector< double > const &a_Xs, bool a_boundIndex=false)
Definition MCGIDI.hpp:318
MCGIDI::particleKineticEnergyFromBeta2
LUPI_HOST_DEVICE double particleKineticEnergyFromBeta2(double a_mass_unitOfEnergy, double a_particleBeta2)
Definition MCGIDI_misc.cc:144
MCGIDI::muCOM_From_muLab
LUPI_HOST_DEVICE int muCOM_From_muLab(double a_muLab, double a_boostBeta, double a_productBeta, double &a_muPlus, double &a_JacobianPlus, double &a_muMinus, double &a_JacobianMinus)
Definition MCGIDI_misc.cc:185
PoPI::Intids::photon
static int constexpr photon
Definition PoPI.hpp:178
PoPI::Intids::neutron
static int constexpr neutron
Definition PoPI.hpp:177