marley::Generator Class Reference

The MARLEY Event generator. More...

#include <Generator.hh>

Public Member Functions
	Generator ()
	Create a Generator using default settings.
void	add_reaction (std::unique_ptr< marley::Reaction > reaction)
	Take ownership of a new Reaction.
void	add_state_to_event (HepMC3::GenEvent &ev) const
	Attach the current random number generator state to the input event as a string attribute.
void	assign_run_info (HepMC3::GenEvent &event) const
	Associates the owned GenRunInfo object with the input event.
void	clear_reactions ()
	Clear the vector of Reaction objects owned by this Generator.
std::shared_ptr< HepMC3::GenEvent >	create_event (bool attach_state=false)
	Create an Event using the NeutrinoSource, Target, Reaction, and StructureDatabase objects owned by this Generator.
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)
	Creates an event object for a fixed projectile species, kinetic energy, and atomic target.
double	E_pdf (double E)
	Probability density function that describes the distribution of reacting neutrino energies.
void	finish_event_metadata (HepMC3::GenEvent &ev, bool attach_state=false)
	Add final pieces of metadata (e.g., the RNG state) to an otherwise complete event.
double	flux_averaged_total_xs () const
	Computes the flux-averaged total cross section for all enabled neutrino reactions, taking target atom fractions into account as appropriate.
const std::vector< std::unique_ptr< marley::Reaction > > &	get_reactions () const
	Get a const reference to the vector of Reaction objects owned by this Generator.
marley::ProjectileDirectionRotator &	get_rotator ()
	Provides access to the owned ProjectileDirectionRotator.
uint_fast64_t	get_seed () const
	Get the seed used to initialize this Generator.
const marley::NeutrinoSource &	get_source () const
	Get a const reference to the NeutrinoSource owned by this Generator.
std::string	get_state_string () const
	Get a string that represents the current internal state of this Generator.
marley::StructureDatabase &	get_structure_db ()
	Get a reference to the StructureDatabase owned by this Generator.
const marley::Target &	get_target () const
	Get a const reference to the Target owned by this Generator.
marley::Weighter &	get_weighter ()
	Get a non-const reference to the owned Weighter object.
const marley::Weighter &	get_weighter () const
	Get a const reference to the owned Weighter object.
double	inverse_transform_sample (const marley::InterpolatingFunction &cdf, double xmin, double xmax, double bisection_tolerance=1e-12)
	Sample from a given 1D cumulative density function cdf(x) on the interval [xmin, xmax] using bisection.
double	inverse_transform_sample (const std::function< double(double)> &f, double xmin, double xmax, double bisection_tolerance=1e-12)
	Sample from a given 1D probability density function f(x) on the interval [xmin, xmax] using an inverse transform technique.
const std::string &	json_config () const
const std::array< double, 3 > &	neutrino_direction ()
	Gets the direction of the incident neutrinos that is used when generating events.
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_)
	Sample from a given 1D probability density function f(x) on the interval [xmin, xmax] using a simple rejection method.
void	reseed (uint_fast64_t seed)
	Reseeds the Generator.
const std::shared_ptr< HepMC3::GenRunInfo > &	run_info () const
double	sample_decay_time (double partial_width)
	Sample a random decay time given a partial decay width.
template<class RandomNumberDistribution>
auto	sample_from_distribution (RandomNumberDistribution &rnd) -> decltype(std::declval< RandomNumberDistribution & >().operator()(std::declval< std::mt19937_64 & >()))
	Sample from an arbitrary probability distribution (defined here as any object that implements an operator()(std::mt19937_64&) function) @detail This template function is based on https://stackoverflow.com/a/9154394/4081973.
template<class RandomNumberDistribution, typename ParamType>
auto	sample_from_distribution (RandomNumberDistribution &rnd, const ParamType &params) -> decltype(std::declval< RandomNumberDistribution & >().operator()(std::declval< std::mt19937_64 & >(), std::declval< const ParamType & >()))
	Sample from an arbitrary probability distribution (defined here as any object that implements an operator()(std::mt19937_64&, const ParamType&) function) using the parameters params.
marley::Reaction &	sample_reaction (double &E)
	Sample a Reaction and an energy for the reacting neutrino.
void	seed_using_state_string (const std::string &state_string)
	Use a string to set this Generator's internal state.
void	set_do_deexcitations (bool do_them)
	Sets the value of the do_deexcitations flag.
void	set_neutrino_direction (const std::array< double, 3 > &dir_vec)
	Sets the direction of the incident neutrinos to use when generating events.
void	set_run_info (const std::shared_ptr< HepMC3::GenRunInfo > &run_info)
void	set_source (std::unique_ptr< marley::NeutrinoSource > source)
	Take ownership of a new NeutrinoSource, replacing any existing source owned by this Generator.
void	set_target (std::unique_ptr< marley::Target > target)
	Take ownership of a new Target, replacing any existing target owned by this Generator.
void	set_up_run_info ()
	Initializes the owned GenRunInfo object that will be used to associate run metadata with the output events.
void	set_weight_flux (bool should_we_weight)
	Sets the value of the weight_flux flag.
void	set_weighter (std::shared_ptr< marley::Weighter > w)
	Replace the owned Weighter object.
double	total_xs (int pdg_a, double KEa) const
	Computes the abundance-weighted total cross section at fixed energy for all configured reactions.
double	total_xs (int pdg_a, double KEa, int pdg_atom) const
	Computes the total cross section at fixed energy for all configured reactions involving a particular target atom.
double	uniform_random_double (double min, double max, bool inclusive)
	Sample a random number uniformly on either [min, max) or [min, max].

Static Public Member Functions
static void	add_state_to_event (HepMC3::GenEvent &ev, const std::string &state)
	Attach a user-supplied random number generator state string to the input event as a string attribute.

Detailed Description

The MARLEY Event generator.

Definition at line 54 of file Generator.hh.

Constructor & Destructor Documentation

Generator()

marley::Generator::Generator

(

)

Create a Generator using default settings.

Definition at line 45 of file Generator.cc.

  : seed_( std::chrono::system_clock::now().time_since_epoch().count() ),
  source_( new marley::MonoNeutrinoSource ),
  structure_db_( new marley::StructureDatabase )
{
  print_logo();
  reseed( seed_ );
}

References reseed().

Member Function Documentation

add_reaction()

void marley::Generator::add_reaction

(

std::unique_ptr< marley::Reaction >

reaction

)

Take ownership of a new Reaction.

Parameters

reaction

A pointer to the new Reaction to use

Definition at line 479 of file Generator.cc.

{
  // If we're passed a nullptr, then don't bother to do anything
  if ( reaction ) {
 
    // Transfer ownership to a new unique_ptr in the reactions vector, leaving
    // the original empty
    reactions_.push_back( std::move(reaction) );
 
    // Add a new entry in the reaction cross sections vector
    total_xs_values_.push_back( 0. );
 
    // TODO: consider adding a check to see whether source_ is non-null.
    // Right now, this shouldn't be possible, but an explicit check might
    // be good.
 
    // Update the neutrino energy probability density function by including the
    // cross section for the new reaction
    normalize_E_pdf();
  }
}

add_state_to_event() [1/2]

void marley::Generator::add_state_to_event

(

HepMC3::GenEvent &

ev

)

const

Attach the current random number generator state to the input event as a string attribute.

Definition at line 880 of file Generator.cc.

                                                                   {
 
  // Query the random number generator for its current state string
  std::string state = this->get_state_string();
 
  // Attach the random number generator state string to the event as a
  // string attribute
  this->add_state_to_event( ev, state );
}

References add_state_to_event(), and get_state_string().

Referenced by add_state_to_event(), marley::CommandHandler::cmd_generate(), and finish_event_metadata().

add_state_to_event() [2/2]

void marley::Generator::add_state_to_event ( HepMC3::GenEvent & ev, const std::string & state )

static

Attach a user-supplied random number generator state string to the input event as a string attribute.

Definition at line 890 of file Generator.cc.

{
  ev.add_attribute( "MARLEY.GeneratorState",
    std::make_shared< HepMC3::StringAttribute >( state )
  );
}

References HepMC3::GenEvent::add_attribute().

assign_run_info()

void marley::Generator::assign_run_info

(

HepMC3::GenEvent &

event

)

const

Associates the owned GenRunInfo object with the input event.

Definition at line 854 of file Generator.cc.

                                                                   {
  event.set_run_info( run_info_ );
}

Referenced by finish_event_metadata().

clear_reactions()

void marley::Generator::clear_reactions

(

)

Clear the vector of Reaction objects owned by this Generator.

Definition at line 502 of file Generator.cc.

                                      {
  reactions_.clear();
  total_xs_values_.clear();
  // Reset the normalization factor to 1. We don't need it until we define
  // one or more new reactions.
  norm_ = 1.;
}

create_event() [1/2]

std::shared_ptr< HepMC3::GenEvent > marley::Generator::create_event

(

bool

attach_state = false

)

Create an Event using the NeutrinoSource, Target, Reaction, and StructureDatabase objects owned by this Generator.

Parameters

attach_state

Whether to attach the random number generator state to the event as a string attribute

Definition at line 77 of file Generator.cc.

{
  // (0) Initialize the run information if it has not been set up yet
  if ( !run_info_ ) {
    this->set_up_run_info();
  }
 
  // (1) Select a reacting neutrino energy and reaction using the
  // flux-weighted total cross section(s)
  double E_nu;
  marley::Reaction& r = sample_reaction( E_nu );
 
  // (2) Create the prompt two-two scattering event using the
  // sampled reaction object
  int pdg_a = source_->get_pid();
  std::shared_ptr< HepMC3::GenEvent > ev = r.create_event( pdg_a, E_nu, *this );
 
  // E.C.2 and E.C.3
  // Save total and reaction cross sections as metadata in the event.
  // NOTE: the expected value for NuHepMC3 E.C.2 is the total cross section for
  // the selected projectile and target, so we loop over relevant reactions
  // manually to avoid imposing target fraction weighting to this particular
  // calculation.
  double totXS = 0.;
  for ( const auto& temp_r : reactions_ ) {
 
    // Skip reactions which involve a different target atom from the selected
    // one
    if ( temp_r->atomic_target().pdg() != r.atomic_target().pdg() ) continue;
 
    // Skip reactions which involve a different projectile
    if ( temp_r->pdg_a() != r.pdg_a() ) continue;
 
    double temp_xsec = temp_r->total_xs( pdg_a, E_nu );
    if ( temp_xsec > 0. ) {
      totXS += temp_xsec;
    }
 
  }
 
  // Convert to picobarns
  totXS *= marley_utils::hbar_c2 * marley_utils::fm2_to_picobarn;
 
  ev->add_attribute( "tot_xs",
    std::make_shared< HepMC3::DoubleAttribute >( totXS )
  );
 
  double procXS = r.total_xs( pdg_a, E_nu ) * marley_utils::hbar_c2
    * marley_utils::fm2_to_picobarn;
 
  ev->add_attribute( "proc_xs",
    std::make_shared< HepMC3::DoubleAttribute >( procXS )
  );
 
  // (3) If needed, de-excite the final-state residue
  if ( do_deexcitations_ ) {
    marley::NucleusDecayer nd;
    nd.process_event( *ev, *this );
  }
 
  // (4) If needed, rotate the event to match the desired projectile direction
  rotator_.process_event( *ev, *this );
 
  // (5) Finish adding metadata to the event object
  this->finish_event_metadata( *ev, attach_state );
 
  // Return the completed event object
  return ev;
}

References HepMC3::GenEvent::add_attribute(), marley::Reaction::atomic_target(), marley::Reaction::create_event(), finish_event_metadata(), marley::TargetAtom::pdg(), marley::Reaction::pdg_a(), marley::NucleusDecayer::process_event(), sample_reaction(), set_up_run_info(), and marley::Reaction::total_xs().

create_event() [2/2]

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

Creates an event object for a fixed projectile species, kinetic energy, and atomic target.

If no energetically-accessible reaction is available for the given input parameters, then a marley::Error will be thrown.

Note

This function is not used as part of the normal MARLEY workflow. It serves as part of an API that enables MARLEY to be interfaced with an external flux driver.

Parameters

pdg_a	The PDG code for the projectile
KEa	The kinetic energy of the projectile (MeV)
pdg_atom	The nuclear PDG code for the atomic target
dir_vec	Direction three-vector of the projectile
attach_state	Whether to attach the random number generator state to the event as a string attribute

Definition at line 667 of file Generator.cc.

{
  // (0) Initialize the run information if it has not been set up yet
  if ( !run_info_ ) {
    this->set_up_run_info();
  }
 
  // (1) Sample a reaction mode from all configured reactions that can handle
  // the given initial-state parameters
  std::vector<size_t> indices;
  std::vector<double> xsecs;
  double tot_xsec = this->total_xs( pdg_a, KEa, pdg_atom, &indices, &xsecs );
 
  if ( xsecs.empty() || tot_xsec <= 0. ) throw marley::Error(
    "Cannot create an event for a projectile with kinetic energy = "
    + std::to_string(KEa) + " MeV and PDG code " + std::to_string(pdg_a)
    + " striking an atom with PDG code " + std::to_string(pdg_atom)
    + ". The total cross section vanishes for all configured reactions" );
 
  // The total cross section values and indices in the full reactions_ vector
  // have already been loaded into temporary vectors, so we can immediately use
  // those to sample a reaction using a discrete distribution.
  std::discrete_distribution< size_t > react_dist( xsecs.begin(), xsecs.end() );
  size_t sampled_index = react_dist( rand_gen_ );
  auto& r = reactions_.at( indices.at(sampled_index) );
 
  // (2) Create the prompt two-two scattering event using the sampled reaction
  // object
  std::shared_ptr< HepMC3::GenEvent > ev = r->create_event( pdg_a, KEa, *this );
 
  // E.C.2 and E.C.3
  // Save total and reaction cross sections as metadata in the event
  double totXS = tot_xsec * marley_utils::hbar_c2
    * marley_utils::fm2_to_picobarn;
 
  ev->add_attribute( "tot_xs",
    std::make_shared< HepMC3::DoubleAttribute >( totXS )
  );
 
  double procXS = xsecs.at( sampled_index ) * marley_utils::hbar_c2
    * marley_utils::fm2_to_picobarn;
 
  ev->add_attribute( "proc_xs",
    std::make_shared< HepMC3::DoubleAttribute >( procXS )
  );
 
  // Do the usual post-processing
 
  // (3) If needed, de-excite the final-state residue
  if ( do_deexcitations_ ) {
    marley::NucleusDecayer nd;
    nd.process_event( *ev, *this );
  }
 
  // (4) If needed, rotate the event to match the desired projectile direction
  // Set the incident neutrino direction for this event
  static marley::ProjectileDirectionRotator my_rotator;
  my_rotator.set_projectile_direction( dir_vec );
 
  // Rotate the coordinate system of the event if needed
  my_rotator.process_event( *ev, *this );
 
  // (5) Finish adding metadata to the event object
  this->finish_event_metadata( *ev, attach_state );
 
  // Return the completed event object
  return ev;
}

References HepMC3::GenEvent::add_attribute(), finish_event_metadata(), marley::NucleusDecayer::process_event(), marley::ProjectileDirectionRotator::process_event(), set_up_run_info(), and total_xs().

E_pdf()

double marley::Generator::E_pdf

(

double

E

)

Probability density function that describes the distribution of reacting neutrino energies.

This function computes the cross-section weighted neutrino flux (normalized to unity between source_->E_min and source_->E_max) including cross-section contributions from all Reactions owned by this Generator. For the distribution of incident neutrino energies, use marley::NeutrinoSource::pdf()

Parameters

E	Total energy of the reacting neutrino

Returns

Probability density (MeV ^-1)

Definition at line 335 of file Generator.cc.

                                        {
 
  // Initialize the return value to zero
  double pdf = 0.;
 
  // Sum all of the reaction total cross sections, saving
  // each individual value along the way. Take weighting
  // by atom fraction in the target material into account.
  for ( size_t j = 0, s = reactions_.size(); j < s; ++j ) {
 
    // Get the current reaction
    const auto& react = reactions_.at( j );
 
    // Compute the total cross section for the current reaction for a single
    // target atom
    double tot_xs = react->total_xs( source_->get_pid(), E );
 
    // If the target_ member has not been initialized, don't bother doing any
    // weighting by atom fraction (equivalent to a weight of unity for all
    // target atoms)
    if ( target_ ) {
      // If it has been configured, then apply the appropriate atom fraction
      // weight from the target as appropriate.
      tot_xs *= target_->atom_fraction( react->atomic_target() );
    }
 
    // Cache the atom-fraction-weighted total cross section for sampling a
    // reaction mode later
    total_xs_values_.at( j ) = tot_xs;
 
    // Add the weighted total cross section value to the total
    pdf += tot_xs;
  }
 
  // Normally, we want to fold the flux with the reaction cross section(s)
  // in order to obtain the distribution of reacting neutrino energies
  if ( weight_flux_ ) {
    // Multiply the total cross section by the neutrino spectrum
    // from the source object to get the (unnormalized) PDF
    // for sampling reacting neutrino energies.
    pdf *= source_->pdf( E );
  }
  else {
    // If the user has specifically requested it, don't weight the
    // energy PDF by the cross section(s), as long as at least one of them
    // is non-vanishing
    if ( pdf <= 0. ) return 0.;
    pdf = source_->pdf( E );
  }
 
  //  Divide by the normalization factor (computed when this source
  //  was made available to the Generator) to obtain the normalized PDF.
  return pdf / norm_;
}

Referenced by sample_reaction().

finish_event_metadata()

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

Add final pieces of metadata (e.g., the RNG state) to an otherwise complete event.

Parameters

attach_state

Whether to attach the random number generator state to the event as a string attribute

Definition at line 827 of file Generator.cc.

{
  // Associate the owned run information with the event
  this->assign_run_info( ev );
 
  // Calculate any needed weight(s) for the event
  weighter_->process_event( ev, *this );
 
  // Add the generator state after the event was completed as a string
  // attribute. This allows resuming an interrupted job from where it left off.
  if ( attach_state ) this->add_state_to_event( ev );
 
  // E.R.5
  // TODO: revisit spatial position when MARLEY is interfaced with a
  // detector geometry simulation
  const std::vector< double > lab_pos = { 0., 0., 0. };
  ev.add_attribute( "lab_pos",
    std::make_shared< HepMC3::VectorDoubleAttribute >( lab_pos )
  );
 
}

References HepMC3::GenEvent::add_attribute(), add_state_to_event(), and assign_run_info().

Referenced by marley::CommandHandler::cmd_decay(), create_event(), and create_event().

flux_averaged_total_xs()

double marley::Generator::flux_averaged_total_xs

(

)

const

Computes the flux-averaged total cross section for all enabled neutrino reactions, taking target atom fractions into account as appropriate.

If flux weighting is disabled (via a call to set_weight_flux()) then this function will return zero

Returns

Total cross section (MeV ^-2)

Definition at line 587 of file Generator.cc.

                                                     {
  // If we've disabled weighting the neutrino energy PDF
  // by the total cross section, just return zero
  if ( !weight_flux_ ) return 0.;
 
  double avg_total_xs = 0.;
 
  // For a monoenergetic source, don't bother to do the full
  // integral
  double Emin = source_->get_Emin();
  double Emax = source_->get_Emax();
  if ( Emin == Emax ) {
    avg_total_xs = norm_ / source_->pdf( Emin );
  }
  else {
    double source_norm = marley_utils::num_integrate(
      [this](double Ev) -> double { return this->source_->pdf(Ev); },
      Emin, Emax );
 
    // Use the precomputed integral of the reacting neutrino energy PDF
    avg_total_xs = norm_ / source_norm;
  }
  return avg_total_xs;
}

Referenced by set_up_run_info().

get_reactions()

const std::vector< std::unique_ptr< marley::Reaction > > & marley::Generator::get_reactions ( ) const

inline

Get a const reference to the vector of Reaction objects owned by this Generator.

Definition at line 468 of file Generator.hh.

{ return reactions_; }

Referenced by marley::StrengthVariationWeightCalculator::create_instances(), and marley::StrengthVariationWeightCalculator::ensure_initialized().

get_rotator()

marley::ProjectileDirectionRotator & marley::Generator::get_rotator ( )

inline

Provides access to the owned ProjectileDirectionRotator.

Definition at line 287 of file Generator.hh.

        { return rotator_; }

get_seed()

uint_fast64_t marley::Generator::get_seed ( ) const

inline

Get the seed used to initialize this Generator.

Definition at line 465 of file Generator.hh.

{ return seed_; }

get_source()

const marley::NeutrinoSource & marley::Generator::get_source

(

)

const

Get a const reference to the NeutrinoSource owned by this Generator.

Throws a marley::Error if this Generator does not own a NeutrinoSource object.

Definition at line 449 of file Generator.cc.

                                                              {
  if ( source_ ) return *source_;
  else throw marley::Error( "Error in marley::Generator::get_source()."
    " The member variable source_ == nullptr." );
}

get_state_string()

std::string marley::Generator::get_state_string

(

)

const

Get a string that represents the current internal state of this Generator.

Definition at line 168 of file Generator.cc.

                                                  {
  std::stringstream ss;
  ss << rand_gen_;
  return ss.str();
}

Referenced by add_state_to_event().

get_structure_db()

marley::StructureDatabase & marley::Generator::get_structure_db

(

)

Get a reference to the StructureDatabase owned by this Generator.

Definition at line 510 of file Generator.cc.

                                                           {
  if ( structure_db_ ) return *structure_db_;
  else throw marley::Error( "Error in marley::Generator::get_structure_db()."
    " The member variable structure_db_ == nullptr." );
}

Referenced by marley::CommandHandler::cmd_decay(), marley::DiscreteNuclearReaction::create_event(), and marley::NucleusDecayer::process_event().

get_target()

const marley::Target & marley::Generator::get_target

(

)

const

Get a const reference to the Target owned by this Generator.

Throws a marley::Error if this Generator does not own a Target object.

Definition at line 455 of file Generator.cc.

                                                      {
  if ( target_ ) return *target_;
  else throw marley::Error( "Error in marley::Generator::get_target()."
    " The member variable target_ == nullptr." );
}

get_weighter() [1/2]

marley::Weighter & marley::Generator::get_weighter ( )

inline

Get a non-const reference to the owned Weighter object.

Definition at line 317 of file Generator.hh.

{ return *weighter_; }

get_weighter() [2/2]

const marley::Weighter & marley::Generator::get_weighter ( ) const

inline

Get a const reference to the owned Weighter object.

Definition at line 313 of file Generator.hh.

        { return *weighter_; }

inverse_transform_sample() [1/2]

double marley::Generator::inverse_transform_sample	(	const marley::InterpolatingFunction &	cdf,
		double	xmin,
		double	xmax,
		double	bisection_tolerance = 1e-12 )

Sample from a given 1D cumulative density function cdf(x) on the interval [xmin, xmax] using bisection.

Parameters

cdf	Cumulative density function to use for sampling
xmin	Lower bound of the sampling interval
xmax	Upper bound of the sampling interval

Returns

Sampled value of x

Definition at line 552 of file Generator.cc.

{
  // Sample a probability value uniformly on [0, 1]
  double prob = uniform_random_double(0., 1., true);
 
  // If we chose an endpoint, we're done, so just return the appropriate one
  if ( prob == 0. ) return xmin;
  else if ( prob == 1. ) return xmax;
 
  // A properly normalized CDF should evaluate to unity at x = xmax. We enforce
  // this here so that the user doesn't have to do it in advance.
  double norm = cdf.evaluate( xmax );
 
  // Find the x value corresponding to the sampled probability via bisection
  // (slow but robust)
  double a = xmin;
  double b = xmax;
  while ( (b - a) > bisection_tolerance ) {
    double midpoint = (a + b) / 2.;
    double mid_cdf = cdf.evaluate( midpoint ) / norm;
    // If the CDF at the midpoint exactly matches our sampled
    // probability, we're done. Just return the midpoint.
    if ( mid_cdf == prob ) return midpoint;
    // Otherwise, shrink the bisection interval and try again
    else if ( mid_cdf > prob ) b = midpoint;
    else a = midpoint; // mid_cdf < prob
  }
 
  // Return the midpoint of the bisection interval as our sampled x value
  double x = (a + b) / 2.;
  return x;
}

References marley::InterpolatingFunction::evaluate(), and uniform_random_double().

Referenced by inverse_transform_sample().

inverse_transform_sample() [2/2]

double marley::Generator::inverse_transform_sample	(	const std::function< double(double)> &	f,
		double	xmin,
		double	xmax,
		double	bisection_tolerance = 1e-12 )

Sample from a given 1D probability density function f(x) on the interval [xmin, xmax] using an inverse transform technique.

Parameters

f	Probability density function to use for sampling
xmin	Lower bound of the sampling interval
xmax	Upper bound of the sampling interval

Returns

Sampled value of x

Todo

Remove hard-coded number of points here

Definition at line 536 of file Generator.cc.

{
  // Build an approximate CDF corresponding to the integral of the input PDF.
  // Use a polynomial approximant at Chebyshev points to do it.
  marley::ChebyshevInterpolatingFunction func(f, xmin, xmax,
    DEFAULT_N_CHEBYSHEV);
  auto cdf = func.cdf();
 
  // Now that we have a CDF to use for sampling, delegate the rest of the
  // action to the overloaded version of this function.
  return this->inverse_transform_sample(cdf, xmin, xmax, bisection_tolerance);
}

References inverse_transform_sample().

json_config()

const std::string & marley::Generator::json_config ( ) const

inline

Definition at line 290 of file Generator.hh.

{ return json_config_; }

neutrino_direction()

const std::array< double, 3 > & marley::Generator::neutrino_direction ( )

inline

Gets the direction of the incident neutrinos that is used when generating events.

Definition at line 470 of file Generator.hh.

    { return rotator_.projectile_direction(); }

rejection_sample()

double marley::Generator::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_ )

Sample from a given 1D probability density function f(x) on the interval [xmin, xmax] using a simple rejection method.

Parameters

f	Probability density function to use for sampling
xmin	Lower bound of the sampling interval
xmax	Upper bound of the sampling interval
max_search_tolerance	Tolerance to use when finding the maximum of f(x) using Brent's method

Returns

Sampled value of x

The rejection method used by this function consists of the following steps:

1. Find the maximum of the function f(x) on [xmin, xmax]
2. Sample an x value uniformly over f(x)'s domain
3. Sample a y value uniformly over [0, max(f(x))]
4. If y <= f(x), accept the sampled x value
5. If y > f(x), reject the sampled x value, and return to step 2 to try again

Note that f(x) does not need to be normalized, but its range must be nonnegative. In the first step, an iterative method (Brent's method) is used to find the maximum of f(x). The iterations will continue until two successive iterations agree within max_search_tolerance on the location of the maximum of f(x). To avoid problems with functions that yield double values that are small compared to a typical value of max_search_tolerance (say, max_search_tolerance = 1e-8, while many neutrino cross sections of interest for MARLEY are less than 1e-40 cm^2), MARLEY normalizes all probability density functions to unity before using rejection sampling.

Definition at line 281 of file Generator.cc.

{
  // If we were passed the value marley_utils::UNKNOWN_MAX for fmax, then this
  // signals that we need to search for the function maximum ourselves.
  // Otherwise, we'll assume that the value passed over is good.
  if ( fmax == marley_utils::UNKNOWN_MAX  ) {
    // This variable will be loaded with the value of x
    // that corresponds to the maximum of f(x).
    // We don't actually use this, but currently it's
    // a required parameter of marley_utils::maximize
    double x_at_max;
 
    // Maximize the function and multiply by a safety factor just
    // in case we didn't quite find the exact peak
    fmax = marley_utils::maximize( f, xmin, xmax, max_search_tolerance,
      x_at_max ) * safety_factor;
  }
 
  MARLEY_LOG( TRACE, "physics.generator.sampling" ) << "rejection_sample:"
    " xmin = " << xmin << ", xmax = " << xmax << ", initial fmax = " << fmax;
 
  double x, y, val;
 
  do {
    // Sample x value uniformly from [xmin, xmax]
    x = uniform_random_double( xmin, xmax, true );
 
    // Sample y uniformly from [0, fmax]
    y = uniform_random_double( 0., fmax, true );
 
    val = f( x );
    MARLEY_LOG( TRACE, "physics.generator.sampling" ) << "rejection_sample:"
      " trial x = " << x << ", y = " << y << ", f(x) = " << val;
    if ( val > fmax ) {
 
      MARLEY_LOG( WARN, "physics.generator.sampling" ) << "PDF value f(x) = "
        << val << " at x = " << x << " exceeded the estimated maximum"
        << " fmax = " << fmax << " during rejection sampling.";
 
      fmax = val * safety_factor;
      MARLEY_LOG( WARN, "physics.generator.sampling" )
        << "A new estimate fmax = " << val * safety_factor
        << " will now be adopted.";
    }
  }
  // Keep sampling until you get a y value less than f(x)
  // (the probability density function evaluated at the sampled value of x)
  while ( y > val );
 
  return x;
}

References uniform_random_double().

Referenced by marley::ElectronReaction::create_event(), marley::DiscreteNuclearReaction::sample_cos_theta_c_cm(), marley::NeutrinoSource::sample_incident_neutrino(), and sample_reaction().

reseed()

void marley::Generator::reseed

(

uint_fast64_t

seed

)

Reseeds the Generator.

Definition at line 156 of file Generator.cc.

                                                 {
  // This is an attempt to do a decent job of seeding the random number
  // generator, but optimally accomplishing this can be tricky (see, for
  // example, http://www.pcg-random.org/posts/cpp-seeding-surprises.html)
  seed_ = seed;
  std::seed_seq seed_sequence{ seed_ };
  rand_gen_.seed( seed_sequence );
 
  MARLEY_LOG( NOTICE, "physics.generator" ) << "Seeded random number generator with "
    << seed_;
}

Referenced by Generator().

run_info()

const std::shared_ptr< HepMC3::GenRunInfo > & marley::Generator::run_info ( ) const

inline

Definition at line 292 of file Generator.hh.

        { return run_info_; }

sample_decay_time()

double marley::Generator::sample_decay_time

(

double

partial_width

)

Sample a random decay time given a partial decay width.

Parameters

partial_width

Partial decay width (MeV) for the decay process of interest

Returns

Time for a decay via this process to occur (MeV ^-1)

Definition at line 859 of file Generator.cc.

                                                                {
  if ( partial_width <= 0. ) throw marley::Error( "Non-positive partial decay"
    " width passed to marley::Generator::sample_decay_time()" );
 
  // Mean lifetime (1/MeV)
  double tau = 1. / partial_width;
 
  // Find the double value that comes immediately after zero. This allows us
  // exclude zero and sample uniformly on (0, 1]. Including zero would lead to
  // the possibility of an infinite decay time. See http://tinyurl.com/n3ocg3p
  // for more information.
  double min_to_use = std::nextafter( 0., std::numeric_limits<double>::max() );
 
  // Choose a random number uniformly on (0, 1].
  double r = this->uniform_random_double( min_to_use, 1., true );
 
  // Decay time in MeV^{-1}
  double t = -tau * std::log( r );
  return t;
}

References uniform_random_double().

sample_from_distribution() [1/2]

template<class RandomNumberDistribution>

auto marley::Generator::sample_from_distribution ( RandomNumberDistribution & rnd ) -> decltype( std::declval<RandomNumberDistribution&>().operator()( std::declval<std::mt19937_64&>()) )

inline

Sample from an arbitrary probability distribution (defined here as any object that implements an operator()(std::mt19937_64&) function) @detail This template function is based on https://stackoverflow.com/a/9154394/4081973.

Definition at line 193 of file Generator.hh.

      {
        return rnd(rand_gen_);
      }

Referenced by marley::CommandHandler::cmd_decay(), marley::ContinuumNuclearReaction::create_event(), marley::DiscreteNuclearReaction::create_event(), marley::HauserFeshbachDecay::sample_exit_channel(), and marley::Level::sample_gamma().

sample_from_distribution() [2/2]

template<class RandomNumberDistribution, typename ParamType>

auto marley::Generator::sample_from_distribution ( RandomNumberDistribution & rnd, const ParamType & params ) -> decltype( std::declval<RandomNumberDistribution&>().operator()( std::declval<std::mt19937_64&>(), std::declval<const ParamType&>() ) )

inline

Sample from an arbitrary probability distribution (defined here as any object that implements an operator()(std::mt19937_64&, const ParamType&) function) using the parameters params.

Definition at line 205 of file Generator.hh.

      {
        return rnd(rand_gen_, params);
      }

sample_reaction()

marley::Reaction & marley::Generator::sample_reaction

(

double &

E

)

Sample a Reaction and an energy for the reacting neutrino.

Parameters

[out]

E

Total energy of the neutrino undergoing the reaction

Returns

Reference to the sampled Reaction owned by this Generator

Definition at line 390 of file Generator.cc.

                                                            {
  if ( reactions_.empty() ) throw marley::Error( "Cannot sample"
    " a reaction in marley::Generator::sample_reaction(). The vector of"
    " marley::Reaction objects owned by this generator is empty." );
 
  // Store the "old" value of E_pdf_max_, i.e., the one we had before calling
  // rejection_sample(). This will be used to check for problems.
  double old_max = E_pdf_max_;
 
  // TODO: protect against source_ changing E_min or E_max after you compute
  // the normalization factor norm_ in marley::Generator::init()
  E = rejection_sample( [this](double E_nu)
    -> double { return this->E_pdf(E_nu); }, source_->get_Emin(),
    source_->get_Emax(), E_pdf_max_ );
 
  // If the value of max changed after the call to rejection_sample() and the
  // old value wasn't UNKNOWN_MAX, then the rejection sampling routine must
  // have encountered a PDF value that was larger than our estimated maximum.
  // Alert the user about this and advise them to change the configuration
  // appropriately to avoid a biased reacting neutrino energy distribution.
  static bool issued_long_error_message = false;
  if ( old_max != marley_utils::UNKNOWN_MAX
    && old_max != E_pdf_max_ )
  {
    if ( !issued_long_error_message ) {
      MARLEY_LOG( ERROR, "physics.generator.sampling" )
        << "Estimation of the maximum PDF value failed when"
        << " using a rejection method to sample reacting neutrino energies.\n"
        << "This may occur when, e.g., an incident neutrino flux"
        << " is used that includes multiple sharp peaks.\n"
        << "To avoid biasing the energy distribution, please rerun the"
        << " simulation after adding the following line to the MARLEY job"
        << " configuration file:\n"
        << "    energy_pdf_max: " << E_pdf_max_ << ",\n"
        << "If this error message persists after raising energy_pdf_max to a"
        << " relatively high value, please contact the MARLEY developers for"
        << " troubleshooting help.";
      issued_long_error_message = true;
    }
    else {
      MARLEY_LOG( ERROR, "physics.generator.sampling" )
        << "The maximum PDF value for sampling reacting"
        << " neutrino energies was exceeded again. The new estimated maximum is"
        << "\n    energy_pdf_max: " << E_pdf_max_ << ',';
    }
  }
 
  // The atom-fraction-weighted total cross section values have already been
  // updated by the final call to E_pdf() during rejection sampling, so we can
  // now sample a reaction type using our discrete distribution object.
  std::discrete_distribution<size_t>::param_type
    params( total_xs_values_.begin(), total_xs_values_.end() );
  size_t r_index = r_index_dist_( rand_gen_, params );
  auto& sampled_reaction = *reactions_.at( r_index );
  MARLEY_LOG( DEBUG, "physics.generator" ) << "Sampled reaction: "
    << sampled_reaction.get_description() << " at E_nu = " << E << " MeV";
  return sampled_reaction;
}

References E_pdf(), and rejection_sample().

Referenced by create_event().

seed_using_state_string()

void marley::Generator::seed_using_state_string

(

const std::string &

state_string

)

Use a string to set this Generator's internal state.

This function is typically used to restore a Generator to a state saved using get_state_string().

Definition at line 148 of file Generator.cc.

{
  // TODO: add error handling here (check that state_string is valid)
  std::stringstream strstr( state_string );
  strstr >> rand_gen_;
}

set_do_deexcitations()

void marley::Generator::set_do_deexcitations ( bool do_them )

inline

Sets the value of the do_deexcitations flag.

This should almost always be true. Use only if you know what you are doing.

Definition at line 473 of file Generator.hh.

    { do_deexcitations_ = do_them; }

set_neutrino_direction()

void marley::Generator::set_neutrino_direction

(

const std::array< double, 3 > &

dir_vec

)

Sets the direction of the incident neutrinos to use when generating events.

Parameters

dir_vec

Vector that points in the direction of the incident neutrinos

Note

The dir_vec passed to this function does not need to be normalized, but it must have at least one nonzero element or a marley::Error will be thrown.

Definition at line 516 of file Generator.cc.

{
  rotator_.set_projectile_direction( dir_vec );
 
  const auto& normalized_dir_vec = rotator_.projectile_direction();
 
  // Print a log message announcing the change of direction
  std::string dir_msg("Incident neutrino direction: (");
  for ( size_t i = 0; i < 3; ++i ) {
    dir_msg += std::to_string( normalized_dir_vec[i] );
    if ( i < 2 ) dir_msg += ", ";
  }
  MARLEY_LOG( INFO, "physics.generator" ) << dir_msg << ')';
}

set_run_info()

void marley::Generator::set_run_info ( const std::shared_ptr< HepMC3::GenRunInfo > & run_info )

inline

Definition at line 298 of file Generator.hh.

                   { run_info_ = run_info; }

set_source()

void marley::Generator::set_source

(

std::unique_ptr< marley::NeutrinoSource >

source

)

Take ownership of a new NeutrinoSource, replacing any existing source owned by this Generator.

Parameters

source

A pointer to the new NeutrinoSource to use

Definition at line 461 of file Generator.cc.

{
  // If we're passed a nullptr, then don't bother to do anything
  if ( source ) {
    // Transfer ownership of the neutrino source the generator, leaving
    // the std::unique_ptr passed to this function null afterwards.
    source_.reset( source.release() );
 
    // Don't bother to renormalize if there are no reactions defined yet
    if ( reactions_.empty() ) return;
 
    // Update the neutrino energy probability density function based on the
    // new source spectrum
    this->normalize_E_pdf();
  }
}

set_target()

void marley::Generator::set_target

(

std::unique_ptr< marley::Target >

target

)

Take ownership of a new Target, replacing any existing target owned by this Generator.

Parameters

source

A pointer to the new Target to use

Definition at line 612 of file Generator.cc.

{
  // If we're passed a nullptr, then don't bother to do anything
  if ( target ) {
    // Transfer ownership of the neutrino target to the generator, leaving
    // the std::unique_ptr passed to this function null afterwards.
    target_.reset( target.release() );
 
    // Don't bother to renormalize if there are no reactions defined yet
    if ( reactions_.empty() ) return;
 
    // Update the neutrino energy probability density function based on the
    // new target composition
    this->normalize_E_pdf();
  }
}

set_up_run_info()

void marley::Generator::set_up_run_info

(

)

Initializes the owned GenRunInfo object that will be used to associate run metadata with the output events.

Definition at line 769 of file Generator.cc.

                                      {
 
  // G.R.1
  run_info_ = std::make_shared< HepMC3::GenRunInfo >();
 
  // G.R.2
  run_info_->add_attribute( "NuHepMC.Version.Major",
    std::make_shared< HepMC3::IntAttribute >(
      marley_hepmc3::NUHEPMC_MAJOR_VERSION )
  );
 
  run_info_->add_attribute( "NuHepMC.Version.Minor",
    std::make_shared< HepMC3::IntAttribute >(
      marley_hepmc3::NUHEPMC_MINOR_VERSION )
  );
 
  run_info_->add_attribute( "NuHepMC.Version.Patch",
    std::make_shared< HepMC3::IntAttribute >(
      marley_hepmc3::NUHEPMC_PATCH_VERSION )
  );
 
  // G.R.3
  run_info_->tools().emplace_back(
    HepMC3::GenRunInfo::ToolInfo{ "MARLEY", MARLEY_VERSION,
      MARLEY_GIT_REVISION }
  );
 
  // G.R.8
  marley_hepmc3::prepare_process_metadata( *run_info_ );
 
  // G.R.9
  marley_hepmc3::prepare_vertex_status_metadata( *run_info_ );
 
  // G.R.10
  marley_hepmc3::prepare_particle_status_metadata( *run_info_ );
 
  // G.R.7
  const auto& wgt_names = weighter_->get_weight_names();
  run_info_->set_weight_names( wgt_names );
 
  // G.R.11
  marley_hepmc3::prepare_non_standard_pdg_code_metadata( *run_info_ );
 
  // G.R.4, G.R.6, G.C.2, G.C.3
  double avg_xsec = this->flux_averaged_total_xs(); // MeV^{-2}
  marley_hepmc3::apply_nuhepmc_runinfo_conventions( *run_info_, avg_xsec );
 
  // Save the information needed to restore an interrupted MARLEY job
  run_info_->add_attribute( "MARLEY.RNGseed",
    std::make_shared< HepMC3::StringAttribute >( std::to_string(seed_) )
  );
 
  run_info_->add_attribute( "MARLEY.JSONconfig",
    std::make_shared< HepMC3::StringAttribute >( json_config_ )
  );
 
}

References flux_averaged_total_xs().

Referenced by marley::CommandHandler::cmd_decay(), create_event(), and create_event().

set_weight_flux()

void marley::Generator::set_weight_flux

(

bool

should_we_weight

)

Sets the value of the weight_flux flag.

This is potentially dangerous. Use only if you know what you are doing.

Definition at line 532 of file Generator.cc.

                                                             {
  weight_flux_ = should_we_weight;
}

set_weighter()

void marley::Generator::set_weighter ( std::shared_ptr< marley::Weighter > w )

inline

Replace the owned Weighter object.

Definition at line 320 of file Generator.hh.

        { weighter_ = w; }

Referenced by marley::OutputFile::merge_reweight_provenance_weights().

total_xs() [1/2]

double marley::Generator::total_xs	(	int	pdg_a,
		double	KEa ) const

Computes the abundance-weighted total cross section at fixed energy for all configured reactions.

Atom fractions in the owned Target are used to perform the weighting by nuclide abundance

Note

This function is not used as part of the normal MARLEY workflow. It exposes the abundance-weighted total cross section for use by the "marley xsec" subcommand

Parameters

pdg_a	The PDG code for the projectile
KEa	The kinetic energy of the projectile (MeV)

Returns

Abundance-weighted total cross section (MeV ^-2 / atom)

Definition at line 740 of file Generator.cc.

                                                              {
 
  // Initialize the return value to zero
  double tot_xsec = 0.;
 
  // Sum all of the reaction total cross sections. Take weighting by atom
  // fraction in the target material into account.
  for ( const auto& react : reactions_ ) {
 
    // Compute the total cross section for the current reaction for a single
    // target atom
    double xsec = react->total_xs( pdg_a, KEa );
 
    // If the target_ member has not been initialized, don't bother doing any
    // weighting by atom fraction (equivalent to a weight of unity for all
    // target atoms)
    if ( target_ ) {
      // If it has been configured, then apply the appropriate atom fraction
      // weight from the target as appropriate.
      xsec *= target_->atom_fraction( react->atomic_target() );
    }
 
    // Add the weighted total cross section value to the total
    tot_xsec += xsec;
  }
 
  return tot_xsec;
}

total_xs() [2/2]

double marley::Generator::total_xs	(	int	pdg_a,
		double	KEa,
		int	pdg_atom ) const

Computes the total cross section at fixed energy for all configured reactions involving a particular target atom.

Atom fractions in the owned Target are ignored by this function.

Note

This function is not used as part of the normal MARLEY workflow. It serves as part of an API that enables MARLEY to be interfaced with an external flux driver.

Parameters

pdg_a	The PDG code for the projectile
KEa	The kinetic energy of the projectile (MeV)
pdg_atom	The nuclear PDG code for the atomic target

Returns

Total cross section (MeV ^-2 / atom)

Definition at line 629 of file Generator.cc.

{
  return this->total_xs( pdg_a, KEa, pdg_atom, nullptr, nullptr );
}

References total_xs().

Referenced by marley::CommandHandler::cmd_xsec(), create_event(), and total_xs().

uniform_random_double()

double marley::Generator::uniform_random_double	(	double	min,
		double	max,
		bool	inclusive )

Sample a random number uniformly on either [min, max) or [min, max].

Parameters

min	Lower bound of the sampling interval
max	Upper bound of the sampling interval
inclusive	Whether the upper bound should be included (true) or excluded (false) from the possible sampling outcomes

Definition at line 235 of file Generator.cc.

{
  // Defaults to sampling from [0,1). We will always
  // explicitly supply the upper and lower bounds to
  // this distribution, so we won't worry about the
  // default setting.
  static std::uniform_real_distribution<double> udist;
 
  double max_to_use;
 
  if ( inclusive ) { // sample from [min, max]
 
    // Find the double value that comes immediately after max. This allows us
    // to sample uniformly on [min, max] rather than [min,max). This trick
    // comes from http://tinyurl.com/n3ocg3p.
    max_to_use = std::nextafter( max, std::numeric_limits<double>::max() );
  }
  else { // sample from [min, max)
    max_to_use = max;
  }
 
  std::uniform_real_distribution<double>::param_type params( min, max_to_use );
 
  // Sample a random double from this distribution
  return udist( rand_gen_, params );
}

Referenced by marley::CommandHandler::cmd_decay(), marley::ContinuumNuclearReaction::create_event(), marley::DiscreteNuclearReaction::create_event(), marley::ElectronReaction::create_event(), marley::DecayScheme::do_cascade(), inverse_transform_sample(), marley::ExitChannel::prepare_products(), rejection_sample(), and sample_decay_time().

---

The documentation for this class was generated from the following files:

• /home/sjg/Desktop/marley/include/marley/Generator.hh
• /home/sjg/Desktop/marley/src/Generator.cc