diff --git a/corsika/detail/framework/random/RNGManager.inl b/corsika/detail/framework/random/RNGManager.inl index 171a8ddcdfb3088e67322f3fdf8b9c20189f73f4..dec9c849372e3eaa67721d011c309a0ada42ec24 100644 --- a/corsika/detail/framework/random/RNGManager.inl +++ b/corsika/detail/framework/random/RNGManager.inl @@ -12,8 +12,8 @@ namespace corsika { - inline void RNGManager::registerRandomStream(std::string const& pStreamName) { - RNG rng; + inline void RNGManager::registerRandomStream(string_type const& pStreamName) { + prng_type rng; if (auto const& it = seeds_.find(pStreamName); it != seeds_.end()) { rng.seed(it->second); @@ -22,7 +22,7 @@ namespace corsika { rngs_[pStreamName] = std::move(rng); } - inline RNG& RNGManager::getRandomStream(std::string const& pStreamName) { + inline RNGManager::prng_type& RNGManager::getRandomStream(string_type const& pStreamName) { if (isRegistered(pStreamName)) { return rngs_.at(pStreamName); } else { // this stream name is not in the map @@ -31,7 +31,7 @@ namespace corsika { } - inline bool RNGManager::isRegistered(std::string const& pStreamName) const { + inline bool RNGManager::isRegistered(string_type const& pStreamName) const { return rngs_.count(pStreamName) > 0; } @@ -45,7 +45,7 @@ namespace corsika { return buffer; } - inline void RNGManager::seedAll(uint64_t vSeed) { + inline void RNGManager::seedAll(seed_type vSeed) { for (auto& entry : rngs_) { entry.second.seed(vSeed++); } } diff --git a/corsika/detail/modules/qgsjetII/qgsjet-II-04.inl b/corsika/detail/modules/qgsjetII/qgsjet-II-04.inl index 56578021db0093f8a3c922ffa10004607520c492..4efc9563212ab966575bced7443c34511937b1fc 100644 --- a/corsika/detail/modules/qgsjetII/qgsjet-II-04.inl +++ b/corsika/detail/modules/qgsjetII/qgsjet-II-04.inl @@ -29,7 +29,7 @@ datadir::datadir(const std::string& dir) { } double qgran_(int&) { - static corsika::RNG& rng = + static corsika::default_prng_type& rng = corsika::RNGManager::getInstance().GetRandomStream("qgran"); std::uniform_real_distribution<double> dist; diff --git a/corsika/framework/core/Cascade.hpp b/corsika/framework/core/Cascade.hpp index 27b665f6c0708e639ae7d8e9a2142bd4f85f6391..73dce0b79a92e1dede62a41c15f13c0ea55ca3ec 100644 --- a/corsika/framework/core/Cascade.hpp +++ b/corsika/framework/core/Cascade.hpp @@ -110,7 +110,7 @@ namespace corsika { TTracking& fTracking; TProcessList& fProcessSequence; TStack& fStack; - corsika::RNG& fRNG = corsika::RNGManager::getInstance().getRandomStream("cascade"); + corsika::default_prng_type& fRNG = corsika::RNGManager::getInstance().getRandomStream("cascade"); }; } // namespace corsika diff --git a/corsika/framework/random/ExponentialDistribution.hpp b/corsika/framework/random/ExponentialDistribution.hpp index 0dc7f4d575f18b6100042ee84e9442475395e0d2..2a2a1b74f3cf8ad7d751c32f15bab3f38c425a5f 100644 --- a/corsika/framework/random/ExponentialDistribution.hpp +++ b/corsika/framework/random/ExponentialDistribution.hpp @@ -12,26 +12,74 @@ n/* #include <random> namespace corsika { - // FIXME: This whole facility needs to re-designed. - // It is not parallel friendly neither polymorphic - // and the streaming management is prone to produce - // huge correlation between the streams template <class TQuantity> class ExponentialDistribution { - using RealType = typename TQuantity::value_type; - std::exponential_distribution<RealType> dist_{1.}; - TQuantity const beta_; + using RealType = typename TQuantity::value_type; + typedef std::exponential_distribution<RealType> distribution_type; + typedef TQuantity quantity_type; public: - ExponentialDistribution(TQuantity beta) - : beta_(beta) {} + ExponentialDistribution()=delete; + + ExponentialDistribution(quantity_type vBeta) + : beta_(vBeta) {} + + ExponentialDistribution(ExponentialDistribution<quantity_type>const& other): + beta_(other.getBeta()) + {} + + ExponentialDistribution<quantity_type>& + operator=(ExponentialDistribution<quantity_type>const& other){ + if( this == &other) return *this; + beta_ = other.getBeta(); + return *this; + } + + /** + * @fn quantity_type getBeta()const + * @brief Get parameter of exponential distribution \f[ \beta e^{-X}\f] + * @pre + * @post + * @return quantity_type + */ + quantity_type getBeta() const { + return beta_; + } + + /** + * @fn void setBeta(quantity_type) + * @brief Set parameter of exponential distribution \f[ \beta e^{-X}\f] + * + * @pre + * @post + * @param vBeta + */ + void setBeta(quantity_type vBeta) { + beta_ = vBeta; + } + + /** + * @fn quantity_type operator ()(Generator&) + * @brief Generate a random number distributed like \f[ \beta e^{-X}\f] + * + * @pre + * @post + * @tparam Generator + * @param g + * @return + */ template <class Generator> - TQuantity operator()(Generator& g) { + quantity_type operator()(Generator& g) { return beta_ * dist_(g); } + + private: + + distribution_type dist_{1.}; + quantity_type beta_; }; } // namespace corsika diff --git a/corsika/framework/random/RNGManager.hpp b/corsika/framework/random/RNGManager.hpp index 0c6ebfd9c67b3f1a64121ce611bf49e50eb98399..d78691fa9942a35f11e0d4a38c80bdc3b53a4a8d 100644 --- a/corsika/framework/random/RNGManager.hpp +++ b/corsika/framework/random/RNGManager.hpp @@ -8,70 +8,134 @@ #pragma once -#include <corsika/framework/utility/Singleton.hpp> -#include <corsika/framework/logging/Logging.h> #include <map> +#include <cstdint> #include <random> #include <string> +#include <corsika/framework/utility/Singleton.hpp> +#include <corsika/framework/logging/Logging.h> + + /*! * With this class modules can register streams of random numbers. */ namespace corsika { - // FIXME: This while facility needs to re-designed. - // It is not parallel friendly neither polymorphic - // and the streaming management is prone to produce - // huge correlation between the streams - - using RNG = std::mt19937; //!< the actual RNG type that will be used - class RNGManager : public corsika::Singleton<RNGManager> { friend class corsika::Singleton<RNGManager>; - std::map<std::string, RNG> rngs_; - std::map<std::string, std::seed_seq> seeds_; + public: - protected: - RNGManager() {} // why ? + typedef std::mt19937_64 prng_type; + typedef std::uint64_t seed_type; + typedef std::string string_type; + typedef std::map<std::string, prng_type> streams_type; + typedef std::map<std::string, std::seed_seq> seeds_type; + + + RNGManager( RNGManager const& ) = default; + + RNGManager& operator=(RNGManager const& ) = delete; - public: /*! * This function is to be called by a module requiring a random-number * stream during its initialization. * * \throws sth. when stream \a pModuleName is already registered */ - void registerRandomStream(std::string const& pStreamName); + inline void registerRandomStream(string_type const& vStreamName); /*! * returns the pre-stored stream of given name \a pStreamName if * available */ - RNG& getRandomStream(std::string const& pStreamName); + inline prng_type& getRandomStream(string_type const& vStreamName); /*! * Check whether a stream has been registered. */ - bool isRegistered(std::string const& pStreamName) const; + inline bool isRegistered(string_type const& vStreamName) const; /*! * dumps the names and states of all registered random-number streams * into a std::stringstream. */ - std::stringstream dumpState() const; + inline std::stringstream dumpState() const; /** * Set explicit seeds for all currently registered streams. The actual seed values * are incremented from \a vSeed. */ - void seedAll(uint64_t vSeed); + inline void seedAll(seed_type vSeed); + + /** + * Set seeds for all currently registered streams. + */ + inline void seedAll(void); //!< seed all currently registered streams with "real" randomness - void seedAll(void); //!< seed all currently registered streams with "real" randomness + /** + * @fn const streams_type getRngs&()const + * @brief Constant access to the streams. + * + * @pre + * @post + * @return RNGManager::streams_type + */ + inline const streams_type& getRngs() const { + return rngs_; + } + + /** + * @fn const seeds_type getSeeds&()const + * @brief Constant access to the seeds. + * + * @pre + * @post + * @return RNGManager::seeds_type + */ + inline const seeds_type& getSeeds() const { + return seeds_; + } + + /** + * @fn streams_type Rngs() + * @brief Non-constant access to the streams. + * + * @pre + * @post + * @return RNGManager::streams_type& + */ + inline streams_type Rngs() { + return rngs_; + } + + /** + * @fn seeds_type Seeds&() + * @brief Non-constant access to seeds. + * + * @pre + * @post + * @return RNGManager::seeds_type& + */ + inline seeds_type& Seeds() { + return seeds_; + } + + protected: + + RNGManager()=default; + + private: + + streams_type rngs_; + seeds_type seeds_; }; + typedef typename RNGManager::prng_type default_prng_type; + } // namespace corsika #include <corsika/detail/framework/random/RNGManager.inl> diff --git a/corsika/framework/random/UniformRealDistribution.hpp b/corsika/framework/random/UniformRealDistribution.hpp index 3ea701b0497bb0db1249f23003f25d8b09ce45a9..49ca93fc5c86c1750d9818b455fecd54f51bcbd7 100644 --- a/corsika/framework/random/UniformRealDistribution.hpp +++ b/corsika/framework/random/UniformRealDistribution.hpp @@ -13,30 +13,108 @@ namespace corsika { - // FIXME: This while facility needs to re-designed. - // It is not parallel friendly neither polymorphic - // and the streaming management is prone to produce - // huge correlation between the streams template <class TQuantity> class UniformRealDistribution { using RealType = typename TQuantity::value_type; - std::uniform_real_distribution<RealType> dist_{RealType(0.), RealType(1.)}; - - TQuantity const min_, max_; + typedef std::uniform_real_distribution<RealType> distribution_type; public: + + typedef TQuantity quantity_type; + + UniformRealDistribution() = delete; + UniformRealDistribution(TQuantity b) - : min_{TQuantity(phys::units::detail::magnitude_tag, 0)} + : min_{quantity_type(phys::units::detail::magnitude_tag, 0)} , max_(b) {} - UniformRealDistribution(TQuantity a, TQuantity b) - : min_(a) - , max_(b) {} + UniformRealDistribution(quantity_type vMin, quantity_type vMax) + : min_(vMin) + , max_(vMax) {} + + UniformRealDistribution(UniformRealDistribution<quantity_type> const& other): + min_(other.getMin()) + , max_(other.getMax()) + {} + + inline UniformRealDistribution<quantity_type>& + operator=(UniformRealDistribution<quantity_type> const& other){ + if( this == &other) return *this; + min_ = other.getMin(); + max_ = other.getMax(); + return *this; + } + + /** + * @fn quantity_type getMax()const + * @brief Get the upper limit. + * + * @pre + * @post + * @return quantity_type + */ + inline quantity_type getMax() const { + return max_; + } + + /** + * @fn void setMax(quantity_type) + * @brief Set the upper limit. + * + * @pre + * @post + * @param vMax + */ + inline void setMax(quantity_type vMax) { + max_ = vMax; + } + + /** + * @fn quantity_type getMin()const + * @brief Get the lower limit. + * + * @pre + * @post + * @return + */ + inline quantity_type getMin() const { + return min_; + } + + /** + * @fn void setMin(quantity_type) + * @brief Set the lower limit. + * + * @pre + * @post + * @param vMin + */ + inline void setMin(quantity_type vMin) { + min_ = vMin; + } + + /** + * @fn quantity_type operator ()(Generator&) + * @brief Generate a random numberin the range [min, max] + * + * @pre + * @post + * @tparam Generator + * @param g + * @return quantity_type + */ template <class Generator> - TQuantity operator()(Generator& g) { + inline quantity_type operator()(Generator& g) { return min_ + dist_(g) * (max_ - min_); } + + private: + + distribution_type dist_{RealType(0.), RealType(1.)}; + + quantity_type min_; + quantity_type max_; }; } // namespace corsika diff --git a/corsika/modules/HadronicElasticModel.hpp b/corsika/modules/HadronicElasticModel.hpp index 6ea5250749b52e0ed7d9ef3c83e1a847d5bdc130..2ec29ab46501eb50df997372d2001ee588d1ba79 100644 --- a/corsika/modules/HadronicElasticModel.hpp +++ b/corsika/modules/HadronicElasticModel.hpp @@ -39,7 +39,7 @@ namespace corsika::hadronic_elastic_model { using eV2 = decltype(square(electronvolt)); using inveV2 = decltype(1 / square(electronvolt)); - corsika::RNG& fRNG = + corsika::default_prng_type& fRNG = corsika::RNGManager::getInstance().getRandomStream("HadronicElasticModel"); inveV2 B(eV2 s) const; diff --git a/corsika/modules/pythia8/Interaction.hpp b/corsika/modules/pythia8/Interaction.hpp index d2639b9f9b182bbfe8a9ffea303ff01fd924f212..abaf5fdf777b35c891b09bceaa8ea94a76d5359a 100644 --- a/corsika/modules/pythia8/Interaction.hpp +++ b/corsika/modules/pythia8/Interaction.hpp @@ -55,7 +55,7 @@ namespace corsika::pythia8 { corsika::EProcessReturn DoInteraction(TProjectile&); private: - corsika::RNG& fRNG = corsika::RNGManager::getInstance().getRandomStream("pythia"); + corsika::default_prng_type& fRNG = corsika::RNGManager::getInstance().getRandomStream("pythia"); Pythia8::Pythia fPythia; Pythia8::SigmaTotal fSigma; const bool fInternalDecays = true; diff --git a/corsika/modules/pythia8/Random.hpp b/corsika/modules/pythia8/Random.hpp index 1cad31b2f4518fc9f86d6fc7d858487e209df028..9a0b6085dc739ebf0b7b323c8ca244d9911bc01e 100644 --- a/corsika/modules/pythia8/Random.hpp +++ b/corsika/modules/pythia8/Random.hpp @@ -19,7 +19,7 @@ namespace corsika::pythia8 { private: std::uniform_real_distribution<double> fDist; - corsika::RNG& fRNG = corsika::RNGManager::getInstance().getRandomStream("pythia"); + corsika::default_prng_type& fRNG = corsika::RNGManager::getInstance().getRandomStream("pythia"); }; } // namespace corsika::pythia8 diff --git a/corsika/modules/qgsjetII/Interaction.hpp b/corsika/modules/qgsjetII/Interaction.hpp index 4e0ed339782636a3373083a553e95305923de037..f272b8af14417550cfcfa1c0a65af0679559bc38 100644 --- a/corsika/modules/qgsjetII/Interaction.hpp +++ b/corsika/modules/qgsjetII/Interaction.hpp @@ -55,7 +55,7 @@ namespace corsika::qgsjetII { corsika::EProcessReturn DoInteraction(TProjectile&); private: - corsika::RNG& fRNG = corsika::RNGManager::getInstance().getRandomStream("qgran"); + corsika::default_prng_type& fRNG = corsika::RNGManager::getInstance().getRandomStream("qgran"); const int maxMassNumber_ = 208; }; diff --git a/corsika/modules/qgsjetII/Random.hpp b/corsika/modules/qgsjetII/Random.hpp index 5941b190e473f9a539354bf0a088f43d5a74204c..3ad727b3b673e65e35866fd0a8283f91835bd660 100644 --- a/corsika/modules/qgsjetII/Random.hpp +++ b/corsika/modules/qgsjetII/Random.hpp @@ -14,7 +14,7 @@ namespace qgsjetII { double rndm_interface() { - static corsika::RNG& rng = + static corsika::default_prng_type& rng = corsika::RNGManager::getInstance().getRandomStream("qgran"); std::uniform_real_distribution<double> dist; return dist(rng); diff --git a/corsika/modules/sibyll/Interaction.hpp b/corsika/modules/sibyll/Interaction.hpp index 05d3ad33dbe0c0634d59e073a1d4ed4cab16bb1b..4b7d5cac95bac3283fe43d115a20a6e0bac6cc76 100644 --- a/corsika/modules/sibyll/Interaction.hpp +++ b/corsika/modules/sibyll/Interaction.hpp @@ -65,7 +65,7 @@ namespace corsika::sibyll { corsika::EProcessReturn DoInteraction(TProjectile&); private: - corsika::RNG& RNG_ = corsika::RNGManager::getInstance().getRandomStream("s_rndm"); + corsika::default_prng_type& RNG_ = corsika::RNGManager::getInstance().getRandomStream("s_rndm"); // FOR NOW keep trackedParticles private, could be configurable std::vector<corsika::Code> const trackedParticles_ = { corsika::Code::PiPlus, corsika::Code::PiMinus, corsika::Code::Pi0, diff --git a/corsika/modules/sibyll/NuclearInteraction.hpp b/corsika/modules/sibyll/NuclearInteraction.hpp index d273d1bf827f59e64b42b5a6ae1d965fc7856bf2..a743cae72c161ec7cc69bc17fe63b6dd03ab4db9 100644 --- a/corsika/modules/sibyll/NuclearInteraction.hpp +++ b/corsika/modules/sibyll/NuclearInteraction.hpp @@ -56,7 +56,7 @@ namespace corsika::sibyll { TEnvironment const& environment_; corsika::sibyll::Interaction& hadronicInteraction_; std::map<corsika::Code, int> targetComponentsIndex_; - corsika::RNG& RNG_ = corsika::RNGManager::getInstance().getRandomStream("s_rndm"); + corsika::default_prng_type& RNG_ = corsika::RNGManager::getInstance().getRandomStream("s_rndm"); static constexpr unsigned int gNSample_ = 500; // number of samples in MC estimation of cross section static constexpr unsigned int gMaxNucleusAProjectile_ = 56; diff --git a/corsika/modules/sibyll/Random.hpp b/corsika/modules/sibyll/Random.hpp index 475d7e5c73f208a8d75c58af2f13ba4e967882c0..56483223210748a639c90298872af68610787344 100644 --- a/corsika/modules/sibyll/Random.hpp +++ b/corsika/modules/sibyll/Random.hpp @@ -14,7 +14,7 @@ namespace sibyll { double rndm_interface() { - static corsika::RNG& rng = + static corsika::default_prng_type& rng = corsika::RNGManager::getInstance().getRandomStream("s_rndm"); std::uniform_real_distribution<double> dist; return dist(rng); diff --git a/corsika/modules/urqmd/Random.hpp b/corsika/modules/urqmd/Random.hpp index 3a1d8e347be09231669311a2fe379bd78a4c3370..e14a457ba462a34ef4e58d13b4fb069eeba5a189 100644 --- a/corsika/modules/urqmd/Random.hpp +++ b/corsika/modules/urqmd/Random.hpp @@ -17,7 +17,7 @@ namespace urqmd { * the random number generator function of UrQMD */ double rndm_interface() { - static corsika::RNG& rng = + static corsika::default_prng_type& rng = corsika::RNGManager::getInstance().getRandomStream("UrQMD"); static std::uniform_real_distribution<double> dist; return dist(rng); diff --git a/corsika/modules/urqmd/UrQMD.hpp b/corsika/modules/urqmd/UrQMD.hpp index 42ad785c84e7fd1d4304c6e1b1d4277fac00c3af..c9f6a227f828cf75b3476975d9d1ab7c395c1cc6 100644 --- a/corsika/modules/urqmd/UrQMD.hpp +++ b/corsika/modules/urqmd/UrQMD.hpp @@ -38,7 +38,7 @@ namespace corsika::urqmd { private: static CrossSectionType GetCrossSection(corsika::Code, corsika::Code, HEPEnergyType, int); - corsika::RNG& fRNG = corsika::RNGManager::getInstance().getRandomStream("UrQMD"); + corsika::default_prng_type& fRNG = corsika::RNGManager::getInstance().getRandomStream("UrQMD"); std::uniform_int_distribution<int> fBooleanDist{0, 1}; };