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
 
// Standard library includes
#include <limits>
#include <memory>
#include <random>
#include <sstream>
#include <vector>
 
// MARLEY includes
#include "marley/GammaStrengthFunctionModel.hh"
#include "marley/InterpolatingFunction.hh"
#include "marley/NeutrinoSource.hh"
#include "marley/NuclearReaction.hh"
#include "marley/LevelDensityModel.hh"
#include "marley/NuclearFormFactor.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 HepMC3 {
  class GenEvent;
  class GenRunInfo;
}
 
namespace marley {
 
  class ChebyshevInterpolatingFunction;
  class JSON;
  class JSONConfig;
  class Weighter;

  class Generator {
 
    // Allow the JSONConfig class to access the private
    // constructors of Generator
    friend JSONConfig;
 
    public:

      Generator();

      std::shared_ptr< HepMC3::GenEvent > create_event(
        bool attach_state = false );

      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::InterpolatingFunction& 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 );

      double sample_decay_time( double partial_width );

      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;

      std::shared_ptr< HepMC3::GenEvent > create_event( int pdg_a, double KEa,
        int pdg_atom, const std::array<double, 3>& dir_vec,
        bool attach_state = false );

      inline marley::ProjectileDirectionRotator& get_rotator()
        { return rotator_; }
 
      inline const std::string& json_config() const { return json_config_; }
 
      inline const std::shared_ptr< HepMC3::GenRunInfo >& run_info() const
        { return run_info_; }

      void assign_run_info( HepMC3::GenEvent& event ) const;
 
      inline void set_run_info( const std::shared_ptr< HepMC3::GenRunInfo >&
        run_info ) { run_info_ = run_info; }

      void set_up_run_info();

      void finish_event_metadata( HepMC3::GenEvent& ev,
        bool attach_state = false );

      inline const marley::Weighter& get_weighter() const
        { return *weighter_; }

      inline marley::Weighter& get_weighter() { return *weighter_; }

      inline void set_weighter(std::shared_ptr<marley::Weighter> w)
        { weighter_ = w; }

      void add_state_to_event( HepMC3::GenEvent& ev ) const;

      static void add_state_to_event( HepMC3::GenEvent& ev,
        const std::string& state );
 
    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;

      std::shared_ptr< HepMC3::GenRunInfo > run_info_;

      std::string json_config_;

      double total_xs(int pdg_a, double KEa, int pdg_atom,
        std::vector<size_t>* index_vec, std::vector<double>* xsec_vec) const;

      void set_json_config( const marley::JSON& jc );

      std::shared_ptr< Weighter > weighter_;
  };
 
  // 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; }
 
}