include/Geant4/MCGIDI_sampling.hpp

0001 /*
0002 # <<BEGIN-copyright>>
0003 # Copyright 2019, Lawrence Livermore National Security, LLC.
0004 # This file is part of the gidiplus package (https://github.com/LLNL/gidiplus).
0005 # gidiplus is licensed under the MIT license (see https://opensource.org/licenses/MIT).
0006 # SPDX-License-Identifier: MIT
0007 # <<END-copyright>>
0008 */
0009
0010 #ifndef MCGIDI_sampling_hpp_included
0011 #define MCGIDI_sampling_hpp_included 1
0012
0013 #include <LUPI_declareMacro.hpp>
0014 #include <MCGIDI_vector.hpp>
0015 #include <MCGIDI_string.hpp>
0016
0017 namespace MCGIDI {
0018
0019 /*
0020 ============================================================
0021 ======================= DomainHash =========================
0022 ============================================================
0023 */
0024 class DomainHash {
0025
0026     private:
0027         std::size_t m_bins;                                                 /**< The number of bins for the hash. */
0028         double m_domainMin;                                                 /**< The minimum domain value for the hash. */
0029         double m_domainMax;                                                 /**< The maximum domain value for the hash. */
0030         double m_u_domainMin;                                               /**< The log of m_domainMin ). */
0031         double m_u_domainMax;                                               /**< The log of m_domainMax ). */
0032         double m_inverse_du;                                                /**< The value *m_bins* / ( *m_u_domainMax* - *m_u_domainMin* ). */
0033
0034     public:
0035         LUPI_HOST_DEVICE DomainHash( );
0036         LUPI_HOST_DEVICE DomainHash( std::size_t a_bins, double a_domainMin, double a_domainMax );
0037         LUPI_HOST_DEVICE DomainHash( DomainHash const &a_domainHash );
0038
0039         LUPI_HOST_DEVICE std::size_t bins( ) const { return( m_bins ); }                     /**< Returns the value of the **m_bins**. */
0040         LUPI_HOST_DEVICE double domainMin( ) const { return( m_domainMin ); }        /**< Returns the value of the **m_domainMax**. */
0041         LUPI_HOST_DEVICE double domainMax( ) const { return( m_domainMax ); }        /**< Returns the value of the **m_domainMax**. */
0042         LUPI_HOST_DEVICE double u_domainMin( ) const { return( m_u_domainMin ); }    /**< Returns the value of the **m_u_domainMin**. */
0043         LUPI_HOST_DEVICE double u_domainMax( ) const { return( m_u_domainMax ); }    /**< Returns the value of the **m_u_domainMax**. */
0044         LUPI_HOST_DEVICE double inverse_du( ) const { return( m_inverse_du ); }      /**< Returns the value of the **m_inverse_du**. */
0045
0046         LUPI_HOST_DEVICE std::size_t index( double a_domain ) const ;
0047         LUPI_HOST_DEVICE Vector<std::size_t > map( Vector<double> const &a_domainValues ) const ;
0048
0049         LUPI_HOST_DEVICE void serialize( LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode );
0050
0051         LUPI_HOST void print( bool a_printValues ) const ;
0052 };
0053
0054 namespace Sampling {
0055
0056 enum class SampledType { firstTwoBody, secondTwoBody, uncorrelatedBody, unspecified, photon };
0057 class ProductHandler;
0058
0059 LUPI_HOST_DEVICE std::size_t evaluationForHashIndex( std::size_t a_hashIndex, Vector<std::size_t> const &a_hashIndices, double a_energy,
0060                 Vector<double> const &a_energies, double *a_energyFraction );
0061
0062 namespace Upscatter {
0063
0064     enum class Model { none, A, B, BSnLimits, DBRC };
0065
0066 /*
0067 ============================================================
0068 ===================== ModelDBRC_data =======================
0069 ============================================================
0070 */
0071
0072 class ModelDBRC_data {
0073
0074     public:
0075         double m_neutronMass;                   /**< The mass of the neutron. */
0076         double m_targetMass;                    /**< The mass of the target. */
0077         Vector<double> m_energies;              /**< The energy grid for the cross section. */
0078         Vector<double> m_crossSections;         /**< The cross sections corresponding to the energy grid. */
0079         Vector<std::size_t> m_hashIndices;      /**< The indicies for the energy hash function. */
0080         MCGIDI::DomainHash m_domainHash;        /**< The hash "function". */
0081
0082     public:
0083         LUPI_HOST_DEVICE ModelDBRC_data( );
0084         LUPI_HOST ModelDBRC_data( double a_neutronMass, double a_targetMass, Vector<double> const &a_energies, Vector<double> const &a_crossSections,
0085                 DomainHash const &a_domainHash );
0086         LUPI_HOST_DEVICE ~ModelDBRC_data( );
0087
0088         LUPI_HOST_DEVICE double evaluate( double a_energy );
0089         LUPI_HOST_DEVICE double targetThermalSpeed( double a_temperature );
0090         LUPI_HOST_DEVICE double crossSectionMax( double a_energy, double a_targetThermalSpeed );
0091
0092         LUPI_HOST_DEVICE void serialize( LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode );
0093 };
0094
0095 LUPI_HOST_DEVICE ModelDBRC_data *serializeModelDBRC_data( LUPI::DataBuffer &a_buffer, LUPI::DataBuffer::Mode a_mode, ModelDBRC_data *a_modelDBRC_data );
0096
0097 }           // End of namespace Upscatter.
0098
0099 /*
0100 ============================================================
0101 ================ ClientRandomNumberGenerator ===============
0102 ============================================================
0103 */
0104 class ClientRandomNumberGenerator {
0105     private:
0106         double (*m_generator)( void * );                    /**< User supplied generator. */
0107         void *m_state;                                      /**< User supplied state. */
0108
0109     public:
0110         LUPI_HOST_DEVICE ClientRandomNumberGenerator( double (*a_generator)( void * ), void *a_state );
0111
0112         LUPI_HOST_DEVICE double (*generator( ))( void * ) { return( m_generator ); }
0113         LUPI_HOST_DEVICE void *state( ) { return( m_state ); }
0114         LUPI_HOST_DEVICE double Double( ) { return( m_generator( m_state ) ); }
0115
0116 // The following are deprecated.
0117         LUPI_HOST_DEVICE double (*rng( ))( void * ) { return( generator( ) ); }
0118         LUPI_HOST_DEVICE void *rngState( ) { return( state( ) ); }
0119         LUPI_HOST_DEVICE double dRng( ) { return( Double( ) ); }
0120 };
0121
0122 /*
0123 ============================================================
0124 =================== Client Code RNG Data ===================
0125 ============================================================
0126 */
0127 class ClientCodeRNGData : public ClientRandomNumberGenerator {
0128
0129     public:
0130         LUPI_HOST_DEVICE ClientCodeRNGData( double (*a_generator)( void * ), void *a_state );
0131 };
0132
0133 /*
0134 ============================================================
0135 =========================== Input ==========================
0136 ============================================================
0137 */
0138
0139 class Input {
0140
0141     friend ProtareSingle;
0142     friend Reaction;
0143     friend MCGIDI::Sampling::ProductHandler;
0144     friend GRIN_capture;
0145     friend GRIN_inelastic;
0146
0147     private:
0148         bool m_wantVelocity = true ;                        /**< See member m_isVelocity in class Product for meaning. This is user input. */
0149
0150         bool m_dataInTargetFrame = false;                   /**< **True** if the data are in the target's frame and **false** otherwise. */
0151         double m_modelTemperature = 0.0;                    /**< The temperature used when sampling product data. For example, in upscatter model A the projectile is boosted into a sampled target's frame and the modelled temperature is 0.0. */
0152         double m_modelEnergy = 0.0;                         /**< The projectile energy used when sampling product data (see comment for member **m_modelTemperature**. */
0153
0154         SampledType m_sampledType = SampledType::uncorrelatedBody;  /**< For internal use only. Set by distributions and used in the method **MCGIDI::Sampling::ProductHandler::add**. */
0155
0156                                                     // The next 2 members are set by the user via the setTemperatureAndEnergy method.
0157         double m_temperature = 0.0;                         /**< The temperature of the material. This member is set by the user. */
0158         double m_energy = 0.0;                              /**< The energy of the projectile. This member is set by the user. */
0159
0160     public:
0161
0162         Upscatter::Model m_upscatterModel = Upscatter::Model::none; /**< The upscatter model to use when sampling a target's velocity. */
0163
0164
0165                                                     // The rest of the members are set by MCGIDI methods.
0166                                                     // These five are used for upscatter model A and the last 4 also used by model B.
0167         double m_projectileBeta = 0.0;                      /**< The beta = speed / c of the projectile. */
0168         double m_muLab = 0.0;                               /**< The cosine of the angle between the projectile's and the sampled target's velocities. */
0169         double m_targetBeta = 0.0;                          /**< The beta = speed / c of the target. */
0170         double m_relativeBeta = 0.0;                        /**< The beta = speed / c of the relative speed between the projectile and the target. */
0171
0172         Reaction const *m_reaction = nullptr;               /**< The current reaction whose products are being sampled. */
0173
0174         double m_projectileMass = 0.0;                      /**< The mass of the projectile. */
0175         double m_targetMass = 0.0;                          /**< The mass of the target. */
0176
0177         GIDI::Frame m_frame = GIDI::Frame::lab;             /**< The frame the product data are returned in. */
0178         int m_numberOfDBRC_rejections = 0;                  /**< For the DBRC upscattering model, this is the number of rejections + 1 per product sample. */
0179
0180         double m_mu = 0.0;                                  /**< The sampled mu = cos( theta ) for the product. */
0181         double m_phi = 0.0;                                 /**< The sampled phi for the product. */
0182
0183         double m_energyOut1 = 0.0;                          /**< The sampled energy of the product. */
0184         double m_px_vx1 = 0.0;                              /**< Variable used for two-body sampling. */
0185         double m_py_vy1 = 0.0;                              /**< Variable used for two-body sampling. */
0186         double m_pz_vz1 = 0.0;                              /**< Variable used for two-body sampling. */
0187
0188         double m_energyOut2 = 0.0;                          /**< The sampled energy of the second product for a two-body interaction. */
0189         double m_px_vx2 = 0.0;                              /**< Variable used for two-body sampling. */
0190         double m_py_vy2 = 0.0;                              /**< Variable used for two-body sampling. */
0191         double m_pz_vz2 = 0.0;                              /**< Variable used for two-body sampling. */
0192
0193         int m_delayedNeutronIndex = -1;                     /**< If the product is a delayed neutron, this is its index. */
0194         double m_delayedNeutronDecayRate = 0.0;             /**< If the product is a delayed neutron, this is its decay rate. */
0195
0196         int m_GRIN_intermediateResidual = -1;               /**< For special GRIN product sampling, this is the GNDS intid of the intermediate residual. */
0197
0198         LUPI_HOST_DEVICE Input( bool a_wantVelocity, Upscatter::Model a_upscatterModel );
0199
0200         LUPI_HOST_DEVICE bool wantVelocity( ) const { return( m_wantVelocity ); }                       /**< Returns the value of the *m_wantVelocity* member. */
0201         LUPI_HOST_DEVICE double temperature( ) const { return( m_temperature ); }                       /**< Returns the value of the *m_temperature* member. */
0202         LUPI_HOST_DEVICE double energy( ) const { return( m_energy ); }                                 /**< Returns the value of the *m_energy* member. */
0203         LUPI_HOST_DEVICE void setTemperatureAndEnergy( double a_temperature, double a_energy );
0204
0205         LUPI_HOST_DEVICE bool dataInTargetFrame( ) const { return( m_dataInTargetFrame ); }             /**< Returns the value of the *m_dataInTargetFrame*. */
0206         LUPI_HOST_DEVICE double modelTemperature( ) const { return( m_modelTemperature ); }             /**< Returns the value of the *m_dataInTargetFrame* member. */
0207         LUPI_HOST_DEVICE double modelEnergy( ) const { return( m_modelEnergy ); }                       /**< Returns the value of the *m_modelEnergy* member. */
0208
0209         SampledType sampledType( ) const { return( m_sampledType ); }                                   /**< Returns the value of the *m_sampledType* member. */
0210         LUPI_HOST_DEVICE void setSampledType( SampledType a_sampledType ) { m_sampledType = a_sampledType; }             /**< Sets the member *m_sampledType* to *a_sampledType*. */
0211 };
0212
0213 /*
0214 ============================================================
0215 ========================== Product =========================
0216 ============================================================
0217 */
0218 class Product {
0219
0220     public:
0221         SampledType m_sampledType;
0222         bool m_isVelocity;                      /**< If true, m_px_vx, m_py_vy and m_pz_vz are velocities otherwise momenta. */
0223         int m_productIntid;                     /**< The intid of the sampled product. */
0224         int m_productIndex;                     /**< The index of the sampled product. */
0225         int m_userProductIndex;                 /**< The user particle index of the sampled product. */
0226         int m_numberOfDBRC_rejections;          /**< For the DBRC upscattering model, this is the number of rejections + 1 per product sample. */
0227         double m_productMass;                   /**< The mass of the sampled product. */
0228         double m_kineticEnergy;                 /**< The kinetic energy of the sampled product. */
0229         double m_px_vx;                         /**< The velocity or momentum along the x-axis of the sampled product. */
0230         double m_py_vy;                         /**< The velocity or momentum along the y-axis of the sampled product. */
0231         double m_pz_vz;                         /**< The velocity or momentum along the z-axis of the sampled product. The z-axis is along the direction of the projectile's velolcity. */
0232         int m_delayedNeutronIndex;              /**< If the product is a delayed neutron, this is its index. */
0233         double m_delayedNeutronDecayRate;       /**< If the product is a delayed neutron, this is its decay rate. */
0234         double m_birthTimeSec;                  /**< Some products, like delayed fission neutrons, are to appear (be born) later. This is the time in seconds that such a particle should be born since the interaction. */
0235 };
0236
0237 /*
0238 ============================================================
0239 ====================== ProductHandler ======================
0240 ============================================================
0241 */
0242 class ProductHandler {
0243
0244     public:
0245         LUPI_HOST_DEVICE ProductHandler( ) {}
0246         LUPI_HOST_DEVICE ~ProductHandler( ) {}
0247
0248         template <typename RNG, typename PUSHBACK>
0249         LUPI_HOST_DEVICE void add( double a_projectileEnergy, int a_productIntid, int a_productIndex, int a_userProductIndex, double a_productMass, Input &a_input,
0250                 RNG && a_rng, PUSHBACK && push_back, bool isPhoton );
0251 };
0252
0253 /*
0254 ============================================================
0255 ================ StdVectorProductHandler ===================
0256 ============================================================
0257 */
0258 #ifdef __CUDACC__
0259
0260 #define MCGIDI_CUDACC_numberOfProducts 1000
0261
0262 class StdVectorProductHandler : public ProductHandler {
0263
0264     private:
0265         std::size_t m_size;
0266         Product m_products[1024];
0267
0268     public:
0269         LUPI_HOST_DEVICE StdVectorProductHandler( ) : m_size( 0 ) { }
0270         LUPI_HOST_DEVICE ~StdVectorProductHandler( ) { }
0271
0272         LUPI_HOST_DEVICE std::size_t size( ) { return( m_size ); }
0273         LUPI_HOST_DEVICE Product &operator[]( long a_index ) { return( m_products[a_index] ); }
0274         LUPI_HOST_DEVICE void push_back( Product &a_product ) {
0275             if( m_size < MCGIDI_CUDACC_numberOfProducts ) {
0276                 m_products[m_size] = a_product;
0277                 ++m_size;
0278             }
0279         }
0280         LUPI_HOST_DEVICE void clear( ) { m_size = 0; }
0281 };
0282
0283 #else
0284 class StdVectorProductHandler : public ProductHandler {
0285
0286     private:
0287         std::vector<Product> m_products;            /**< The list of products sampled. */
0288
0289     public:
0290         LUPI_HOST_DEVICE StdVectorProductHandler( ) : m_products( ) { }
0291         LUPI_HOST_DEVICE ~StdVectorProductHandler( ) { }
0292
0293         LUPI_HOST_DEVICE std::size_t size( ) { return( m_products.size( ) ); }
0294         LUPI_HOST_DEVICE Product &operator[]( std::size_t a_index ) { return( m_products[a_index] ); }
0295         LUPI_HOST_DEVICE std::vector<Product> &products( ) { return( m_products ); }
0296         LUPI_HOST_DEVICE void push_back( Product &a_product ) { m_products.push_back( a_product ); }
0297         LUPI_HOST_DEVICE void clear( ) { m_products.clear( ); }
0298 };
0299 #endif
0300
0301 /*
0302 ============================================================
0303 ============== MCGIDIVectorProductHandler ==================
0304 ============================================================
0305 */
0306 class MCGIDIVectorProductHandler : public ProductHandler {
0307
0308     private:
0309         Vector<Product> m_products;             /**< The list of products sampled. */
0310
0311     public:
0312         LUPI_HOST_DEVICE MCGIDIVectorProductHandler( std::size_t a_size = 20 ) :
0313                 m_products( ) {
0314
0315             m_products.reserve( a_size );
0316         }
0317         LUPI_HOST_DEVICE ~MCGIDIVectorProductHandler( ) {}
0318
0319         LUPI_HOST_DEVICE std::size_t size( ) { return( m_products.size( ) ); }
0320         LUPI_HOST_DEVICE Product const &operator[]( std::size_t a_index ) const { return( m_products[a_index] ); }
0321         LUPI_HOST_DEVICE Vector<Product> const &products( ) const { return( m_products ); }
0322         LUPI_HOST_DEVICE void push_back( Product &a_product ) { m_products.push_back( a_product ); }
0323         LUPI_HOST_DEVICE void clear( ) { m_products.clear( ); }
0324 };
0325
0326 }       // End of namespace Sampling.
0327
0328 }       // End of namespace MCGIDI.
0329
0330 #endif      // End of MCGIDI_sampling_hpp_included