Generator.hh

//
// This file is part of MARLEY (Model of Argon Reaction Low Energy Yields)
//
// MARLEY is free software: you can redistribute it and/or modify it under the
// terms of version 3 of the GNU General Public License as published by the
// Free Software Foundation.
//
// For the full text of the license please see COPYING or
// visit http://opensource.org/licenses/GPL-3.0
//
// Please respect the MCnet academic usage guidelines. See GUIDELINES
// or visit https://www.montecarlonet.org/GUIDELINES for details.
 
#pragma once
#include <limits>
#include <memory>
#include <random>
#include <sstream>
#include <vector>
 
#include "marley/GammaStrengthFunctionModel.hh"
#include "marley/NeutrinoSource.hh"
#include "marley/NuclearReaction.hh"
#include "marley/LevelDensityModel.hh"
#include "marley/OpticalModel.hh"
#include "marley/Parity.hh"
#include "marley/ProjectileDirectionRotator.hh"
#include "marley/RotationMatrix.hh"
#include "marley/StructureDatabase.hh"
#include "marley/Target.hh"
#include "marley/marley_utils.hh"
 
namespace marley {
 
  class ChebyshevInterpolatingFunction;
  class JSONConfig;
 
  class Generator {
 
    // Allow the JSONConfig class to access the private
    // constructors of Generator
    friend JSONConfig;
 
    public:
 
      Generator();
 
      marley::Event create_event();
 
      inline uint_fast64_t get_seed() const;
 
      void reseed(uint_fast64_t seed);
 
      void seed_using_state_string(const std::string& state_string);
 
      std::string get_state_string() const;
 
      double uniform_random_double(double min, double max, bool inclusive);
 
      double rejection_sample(const std::function<double(double)>& f,
        double xmin, double xmax, double& fmax, double safety_factor = 1.01,
        double max_search_tolerance = DEFAULT_REJECTION_SAMPLING_TOLERANCE_);
 
      double inverse_transform_sample(
        const marley::ChebyshevInterpolatingFunction& cdf,
        double xmin, double xmax, double bisection_tolerance = 1e-12);
 
      double inverse_transform_sample(const std::function<double(double)>& f,
        double xmin, double xmax, double bisection_tolerance = 1e-12);
 
      marley::StructureDatabase& get_structure_db();
 
      inline const std::vector< std::unique_ptr<marley::Reaction> >&
        get_reactions() const;
 
      void add_reaction(std::unique_ptr<marley::Reaction> reaction);
 
      void clear_reactions();
 
      marley::Reaction& sample_reaction(double& E);
 
      double E_pdf(double E);
 
      const marley::NeutrinoSource& get_source() const;
 
      void set_source(std::unique_ptr<marley::NeutrinoSource> source);
 
      void set_target(std::unique_ptr<marley::Target> target);
 
      const marley::Target& get_target() const;
 
      template <class RandomNumberDistribution>
        inline auto sample_from_distribution(RandomNumberDistribution& rnd)
        -> decltype( std::declval<RandomNumberDistribution&>().operator()(
        std::declval<std::mt19937_64&>()) )
      {
        return rnd(rand_gen_);
      }
 
      template <class RandomNumberDistribution, typename ParamType>
        inline auto sample_from_distribution(RandomNumberDistribution& rnd,
        const ParamType& params) -> decltype(
        std::declval<RandomNumberDistribution&>().operator()(
        std::declval<std::mt19937_64&>(), std::declval<const ParamType&>() ) )
      {
        return rnd(rand_gen_, params);
      }
 
      void set_neutrino_direction(const std::array<double, 3>& dir_vec);
 
      inline const std::array<double, 3>& neutrino_direction();
 
      void set_weight_flux(bool should_we_weight );
 
      inline void set_do_deexcitations( bool do_them );
 
      double flux_averaged_total_xs() const;
 
      double total_xs(int pdg_a, double KEa, int pdg_atom) const;
 
      double total_xs(int pdg_a, double KEa) const;
 
      marley::Event create_event( int pdg_a, double KEa, int pdg_atom,
        const std::array<double, 3>& dir_vec );
 
      inline marley::ProjectileDirectionRotator& get_rotator()
        { return rotator_; }
 
    private:
 
      Generator(uint_fast64_t seed);
 
      void normalize_E_pdf();
 
      void print_logo();
 
      uint_fast64_t seed_;
 
      std::mt19937_64 rand_gen_;
 
      static constexpr double DEFAULT_REJECTION_SAMPLING_TOLERANCE_ = 1e-8;
 
      double norm_ = 1.;
 
      std::unique_ptr<marley::NeutrinoSource> source_;
 
      std::unique_ptr<marley::Target> target_;
 
      std::unique_ptr<marley::StructureDatabase> structure_db_;
 
      std::vector< std::unique_ptr<marley::Reaction> > reactions_;
 
      std::vector<double> total_xs_values_;
 
      std::discrete_distribution<size_t> r_index_dist_;
 
      bool weight_flux_ = true;
 
      double E_pdf_max_;
 
      double E_PDF_MAX_DEFAULT_ = marley_utils::UNKNOWN_MAX;
 
      inline void set_default_E_pdf_max( double def_max ) {
        E_PDF_MAX_DEFAULT_ = def_max;
      }
 
      bool dont_normalize_E_pdf_ = false;
 
      marley::ProjectileDirectionRotator rotator_;
 
      bool do_deexcitations_ = true;
 
      double total_xs(int pdg_a, double KEa, int pdg_atom,
        std::vector<size_t>* index_vec, std::vector<double>* xsec_vec) const;
  };
 
  // Inline function definitions
  inline uint_fast64_t Generator::get_seed() const { return seed_; }
 
  inline const std::vector<std::unique_ptr<marley::Reaction> >&
    Generator::get_reactions() const { return reactions_; }
 
  inline const std::array<double, 3>& Generator::neutrino_direction()
    { return rotator_.projectile_direction(); }
 
  inline void Generator::set_do_deexcitations( bool do_them )
    { do_deexcitations_ = do_them; }
}