marley::ContinuumNuclearReaction Class Reference

Generates inclusive scattering events and computes cross sections using tabulated nuclear responses. More...

#include <ContinuumNuclearReaction.hh>

Inheritance diagram for marley::ContinuumNuclearReaction:

Public Types
enum class	SubContinuumMode { IGNORE , MIRROR , ACCUMULATE }
Public Types inherited from marley::Reaction
enum	DataFormat { DiscreteStrengths = 0 , MultipoleResponses = 1 }
	Enumerated type describing the file format for reaction data. More...
enum	ProcessType {   Unknown = -1 , NeutrinoCC_Discrete = 0 , AntiNeutrinoCC_Discrete = 1 , NC_Discrete = 2 ,   NuElectronElastic = 3 , NeutrinoCC_Continuum = 4 , AntiNeutrinoCC_Continuum = 5 , NC_Continuum = 6 ,   StandaloneDecay = 7 }
	Enumerated type describing the kind of scattering process represented by a Reaction. More...

Public Member Functions
	ContinuumNuclearReaction (Reaction::ProcessType pt, int pdg_a, int pdg_b, int pdg_c, int pdg_d, int q_d, const std::shared_ptr< TabulatedXSec > &txsec, const std::string &source_file)
virtual std::shared_ptr< HepMC3::GenEvent >	create_event (int particle_id_a, double KEa, marley::Generator &gen) const override
	Create an event object for this reaction.
const TabulatedXSec &	get_tabulated_xsec () const
virtual double	total_xs (int pdg_a, double KEa) const override
	Compute the reaction's total cross section (MeV -2)
Public Member Functions inherited from marley::NuclearReaction
	NuclearReaction (ProcessType pt, int pdg_a, int pdg_b, int pdg_c, int pdg_d, int q_d, const std::string &source_file)
virtual marley::TargetAtom	atomic_target () const override final
	Returns the target atom involved in this reaction.
double	max_level_energy (double KEa) const
	Get the maximum possible excitation energy (MeV) of the final-state residue that is kinematically allowed.
double	threshold_kinetic_energy () const override
	Get the minimum lab-frame kinetic energy (MeV) of the projectile that allows this reaction to proceed via a transition to the residue's ground state.
double	weak_nuclear_charge () const
Public Member Functions inherited from marley::Reaction
	Reaction (const std::string &source_file)
	Construct a Reaction with the resolved path of the data file.
const std::string &	get_description () const
	Get a string that contains the formula for this reaction.
int	pdg_a () const
	Get the projectile PDG code.
int	pdg_b () const
	Get the target PDG code.
ProcessType	process_type () const
	Get the process type for this reaction.
const std::string &	source_file () const
	Get the resolved path of the reaction data file used to construct this Reaction.

Static Public Member Functions
static void	set_sub_continuum_mode (SubContinuumMode scm)
	Sets the approach to handling sub-continuum cross-section strength.
static std::string	string_from_sub_continuum_mode (const SubContinuumMode sc_mode)
static SubContinuumMode	sub_continuum_mode ()
	Gets the approach to handling sub-continuum cross-section strength.
static SubContinuumMode	sub_continuum_mode_from_string (const std::string &sc_mode_str)
Static Public Member Functions inherited from marley::Reaction
static int	get_ejectile_pdg (int pdg_a, ProcessType proc_type)
static void	get_residue_pdg_and_charge (ProcessType proc_type, int pdg_b, int &pdg_d, int &q_d)
static std::vector< std::unique_ptr< Reaction > >	load_from_file (const std::string &filename, StructureDatabase &db, CoulombCorrector::CoulombMode coulomb_mode, const JSON &ff_config)
static std::string	proc_type_to_string (const ProcessType &pt)

Protected Member Functions
double	get_Ec_from_Ex (const double Ex, const double cos_theta, const double KEa, double *jacobian=nullptr) const
	Helper function for reassign_sub_continuum() that solves for the outgoing lepton total energy that corresponds to the input kinematic variables.
bool	reassign_sub_continuum (double &w, const double ctl, const double KEa) const
	Helper function for create_event() that potentially reassigns the value of the energy transfer.
virtual void	set_description () override
	Creates the description string based on the PDG code values for the initial and final particles.
Protected Member Functions inherited from marley::NuclearReaction
virtual std::shared_ptr< HepMC3::GenEvent >	make_nuclear_event_object (double KEa, const std::shared_ptr< HepMC3::GenParticle > &ejectile, const std::shared_ptr< HepMC3::GenParticle > &residue, double E_level, int twoJ, const marley::Parity &P) const
	Helper function that makes a complete event object for a nuclear reaction.
virtual std::shared_ptr< HepMC3::GenEvent >	make_nuclear_event_object (double KEa, double pc_cm, double cos_theta_c_cm, double phi_c_cm, double Ec_cm, double Ed_cm, double E_level, int twoJ, const marley::Parity &P) const
	Helper function that makes a complete event object for a nuclear reaction.
void	set_charge_attributes (std::shared_ptr< HepMC3::GenEvent > &event) const
	Helper function that sets the charges of the target and residue in an otherwise complete event record.
void	set_nuclear_residue_attributes (const std::shared_ptr< HepMC3::GenParticle > &residue, double E_level, int twoJ, const marley::Parity &P) const
	Helper function that adds the nuclear level attributes ( \( E_x \), \( 2J \), and parity) needed to keep track of the residue's de-excitation state to an otherwise complete event record.
Protected Member Functions inherited from marley::Reaction
virtual std::shared_ptr< HepMC3::GenEvent >	make_event_object (double KEa, const std::shared_ptr< HepMC3::GenParticle > &ejectile, const std::shared_ptr< HepMC3::GenParticle > &residue) const
	Helper function that makes an event object.
virtual std::shared_ptr< HepMC3::GenEvent >	make_event_object (double KEa, double pc_cm, double cos_theta_c_cm, double phi_c_cm, double Ec_cm, double Ed_cm, int residue_status) const
	Helper function that makes an event object.
void	two_two_scatter (double KEa, double &s, double &Ec_cm, double &pc_cm, double &Ed_cm) const
	Helper function that handles CM frame kinematics for the reaction.

Protected Attributes
std::shared_ptr< TabulatedXSec >	xsec_
	Helper object that handles cross section calculations.
Protected Attributes inherited from marley::NuclearReaction
int	Af_
	Residue mass number.
int	Ai_
	Target mass number.
double	KEa_threshold_
	Lab-frame kinetic energy of the projectile at threshold for this reaction (i.e., the residue is produced in its ground state, and all final-state particles are at rest in the CM frame)
double	md_gs_
	Ground state mass (MeV) of the residue.
int	q_d_
	Net charge of the residue (in units of the proton charge) following this reaction.
int	Zf_
	Residue atomic number.
int	Zi_
	Target atomic number.
Protected Attributes inherited from marley::Reaction
std::string	description_
	String that contains a formula describing the reaction.
double	ma_
	Projectile mass (MeV)
double	mb_
	Target mass (MeV)
double	mc_
	Ejectile mass (MeV)
double	md_
	Residue mass (MeV)
int	pdg_a_
	PDG code for the projectile.
int	pdg_b_
	PDG code for the target.
int	pdg_c_
	PDG code for the ejectile.
int	pdg_d_
	PDG code for the residue.
ProcessType	process_type_
	Type of scattering process (CC, NC) represented by this reaction.
const std::string	source_file_
	Resolved path of the reaction data file.

Static Protected Attributes
static SubContinuumMode	sc_mode_ = SubContinuumMode::ACCUMULATE
	Indicates the desired method for handling events with excitation energies originally sampled below the continuum threshold.
static const std::map< SubContinuumMode, std::string >	sc_mode_string_map_
	Helper map used for conversions between a SubContinuumMode value and a std::string.

Additional Inherited Members
Static Protected Member Functions inherited from marley::Reaction
static const std::vector< int > &	get_projectiles (ProcessType proc_type)

Detailed Description

Generates inclusive scattering events and computes cross sections using tabulated nuclear responses.

Definition at line 30 of file ContinuumNuclearReaction.hh.

Member Enumeration Documentation

SubContinuumMode

enum class marley::ContinuumNuclearReaction::SubContinuumMode

strong

Defines the approach to handling possible reassignment of cross-section strength that falls below the excitation energy threshold for the continuum

Definition at line 37 of file ContinuumNuclearReaction.hh.

{ IGNORE, MIRROR, ACCUMULATE };

Constructor & Destructor Documentation

ContinuumNuclearReaction()

marley::ContinuumNuclearReaction::ContinuumNuclearReaction	(	Reaction::ProcessType	pt,
		int	pdg_a,
		int	pdg_b,
		int	pdg_c,
		int	pdg_d,
		int	q_d,
		const std::shared_ptr< TabulatedXSec > &	txsec,
		const std::string &	source_file )

Definition at line 45 of file ContinuumNuclearReaction.cc.

  : marley::NuclearReaction( pt, pdg_a, pdg_b, pdg_c, pdg_d, q_d,
      source_file ),
  xsec_( txsec )
{
}

Member Function Documentation

create_event()

std::shared_ptr< HepMC3::GenEvent > marley::ContinuumNuclearReaction::create_event ( int pdg_a, double KEa, marley::Generator & gen ) const

overridevirtual

Create an event object for this reaction.

Parameters

pdg_a	PDG code for the incident projectile
KEa	Lab-frame kinetic energy of the projectile
gen	Reference to the Generator to use for random sampling

Note

Functions that override create_event() should throw an Error if pdg_a != pdg_a_.

Implements marley::Reaction.

Definition at line 62 of file ContinuumNuclearReaction.cc.

{
  // TODO: reduce code duplication here with DiscreteNuclearReaction using the
  // common base class NuclearReaction
 
  // Check that the projectile supplied to this event is correct. If not, alert
  // the user that this event does not use the requested projectile.
  if ( pdg_a != pdg_a_ ) throw marley::Error( "Could not create this event."
    " The requested projectile particle ID, " + std::to_string( pdg_a )
    + ", does not match the projectile particle ID, "
    + std::to_string( pdg_a_ ) + ", in the reaction dataset." );
 
  // Sample a final residue energy level. First, check to make sure the given
  // projectile energy is above threshold for this reaction.
  if ( KEa < KEa_threshold_ ) throw std::range_error( "Could"
    " not create this event. Projectile kinetic energy "
    + std::to_string( KEa ) + " MeV is below the threshold value "
    + std::to_string( KEa_threshold_ ) + " MeV." );
 
  // Select a specific multipole to use for the current event using the
  // individual total cross sections
  std::vector< double > multipole_weights;
  std::vector< marley::TabulatedXSec::MultipoleLabel > multipoles;
  std::vector< double > diff_max_values;
  const auto& table_map = xsec_->get_table_map();
  double sum_of_xsecs = 0.;
  for ( const auto& pair : table_map ) {
    const auto& ml = pair.first;
    double diff_max;
    double total_xsec = xsec_->integral( pdg_a_, KEa, ml, diff_max );
 
    sum_of_xsecs += total_xsec;
 
    multipole_weights.push_back( total_xsec );
    multipoles.push_back( ml );
    diff_max_values.push_back( diff_max );
  }
 
  // If there are no multipole weights, we can't go on. Complain if this
  // is the case.
  if ( multipole_weights.empty() ) {
    throw marley::Error( "Could not create this event. The TabulatedXSec"
      " object associated with this reaction does not own any nuclear response"
      " tables." );
  }
 
  // Complain if the total cross section (the sum of all partial cross
  // sections) is zero or negative (the latter is just to cover all
  // possibilities).
  if ( sum_of_xsecs <= 0. ) {
    throw marley::Error( "Could not create this event. All multipole total"
      " cross sections are nonpositive." );
  }
 
  // Create a discrete distribution based on the weights. This will be
  // used to choose a single multipole for the current event.
  std::discrete_distribution< size_t > multipole_dist(
    multipole_weights.begin(), multipole_weights.end() );
 
  // Sample a matrix_element using our discrete distribution and the
  // current set of weights
  size_t multipole_index = gen.sample_from_distribution( multipole_dist );
 
  // Label of the multipole chosen for this event
  const auto& sampled_ml = multipoles.at( multipole_index );
  // Maximum value of the differential cross section for this multipole.
  // This will be used for rejection sampling of inclusive kinematics below.
  double diff_max = diff_max_values.at( multipole_index );
 
  // ResponseTable object to use for computing the differential cross section
  // during kinematic sampling below
  const auto& rt = xsec_->get_table( sampled_ml );
 
  // Get the values of the energy transfer that correspond to the edges
  // of the table of nuclear responses. Note that the table is actually given
  // in terms of the effective energy transfer, which differs by delta_ias
  // from the actual energy transfer. We therefore apply a shift here
  // to correct for this.
  double table_wmin = rt.w_min() - xsec_->delta_ias();
  double table_wmax = rt.w_max() - xsec_->delta_ias();
 
  // Choose a reasonable sampling interval for the energy transfer
  double Ea = KEa + ma_; // Projectile total energy
 
  // Set the lower bound for the energy transfer to be either zero or
  // the lowest tabulated value (whichever is larger)
  double wmin = std::max( 0., table_wmin );
 
  // Set the upper bound for the energy transfer to be either the projectile
  // energy minus the final lepton mass or the highest tabulated value
  // (whichever is lower)
  double wmax = std::min( Ea - mc_, table_wmax );
 
  // Sample values for the energy transfer and scattering cosine using the
  // differential cross section for the chosen multipole.
  // Use a simple rejection sampling technique.
  double w, ctl, diff, y;
  int sampling_attempts = 0;
  bool recomputed_diff_max = false;
  do {
    // Occasionally the value of diff_max retrieved above can be extremely
    // overestimated when using an optimized version of the total cross
    // section calculation (which relies on interpolation). This typically
    // happens very close to threshold and leads to sampling getting stuck
    // due to a very low acceptance efficiency.
    //
    // To guard against this situation, when the number of sampling attempts
    // exceeds a large value, the estimate of the maximum differential cross
    // section diff_max is recalculated at exactly the input projectile kinetic
    // energy rather than relying on the precomputed value. The updated
    // estimate is then used in a new set of sampling attempts.
    if ( sampling_attempts > marley_utils::LARGE_NUMBER_OF_ITERATIONS ) {
 
      if ( recomputed_diff_max ) {
        // If we make it here, then we are still hitting a huge number of
        // iterations in this sampling loop despite recalculating the maximum
        // differential cross section. This suggests that we are stuck in an
        // infinite loop, so abort with an exception indicating the problem.
        throw marley::Error( "Reached maximum iteration count for rejection"
          " sampling in marley::ContinuumNuclearReaction::create_event()" );
      }
 
      // The value of diff_max is updated by this call to
      // marley::TabulatedXSec::compute_integral()
      xsec_->compute_integral( pdg_a_, KEa, sampled_ml, diff_max );
      sampling_attempts = 0;
      recomputed_diff_max = true;
    }
 
    w = gen.uniform_random_double( wmin, wmax, true );
    ctl = gen.uniform_random_double( -1., 1., true );
    diff = xsec_->diff_xsec( pdg_a_, KEa, w, ctl, sampled_ml );
    y = gen.uniform_random_double( 0., diff_max, true );
    ++sampling_attempts;
 
    // If reassignment is enabled and the excitation energy corresponding to the
    // sampled energy transfer w falls below the continuum, reassign the value
    // of the energy transfer to lie within the continuum. The scattering cosine
    // ctl and projectile kinetic energy are needed for the calculation but are
    // never altered. If reassignment is attempted but fails due to kinematic
    // limits, then force another iteration of this rejection sampling loop even
    // if the event would have been accepted without reassignment.
    //
    // NOTE: The reassignment operation is included in the condition of the
    // do-while loop for efficiency. There is no need to perform the
    // reassignment for w values that would be rejected anyway. Use of the
    // logical OR operation (||) will only evaluate the first condition if it
    // is false, thus skipping the reassignment when it is obviously
    // unnecessary.
  } while ( y > diff || !this->reassign_sub_continuum(w, ctl, KEa) );
 
  // Sample a lab-frame azimuthal scattering angle uniformly
  double phi_c = gen.uniform_random_double( 0., marley_utils::two_pi, false );
 
  // Load the initial residue twoJ and parity values into twoJ and P. These
  // variables are included in the event record and used by NucleusDecayer to
  // start the Hauser-Feshbach decay cascade.
  // NOTE: right now, these are taken directly from the multipole involved in
  // the current event. This is only valid for scattering on a 0+ target
  // nucleus
  // TODO: revisit this assumption and do something better
  int twoJ = 2 * sampled_ml.J_;
  marley::Parity P = sampled_ml.Pi_; // defaults to positive parity
 
  // Sine of the ejectile scattering angle
  double stl = marley_utils::real_sqrt( 1. - std::pow(ctl, 2) );
 
  // Calculate the full kinematics of the primary interaction based on the
  // lepton scattering cosine (ctl) and energy transfer (w) sampled above.
 
  // Determine the components of the ejectile's lab-frame 4-momentum
  double Ec = Ea - w;
  double pc = marley_utils::real_sqrt( Ec*Ec - mc_*mc_ );
  double pc_x = stl * std::cos( phi_c ) * pc;
  double pc_y = stl * std::sin( phi_c ) * pc;
  double pc_z = ctl * pc;
 
  // Determine the magnitude of the lab-frame 3-momentum of the projectile
  double pa = marley_utils::real_sqrt( Ea*Ea - ma_*ma_ );
 
  // Construct the lab-frame 4-momenta of the projectile, target, and ejectile
  HepMC3::FourVector pro_mom4( 0., 0., pa, Ea );
  HepMC3::FourVector tar_mom4( 0., 0., 0., mb_ );
  HepMC3::FourVector eje_mom4( pc_x, pc_y, pc_z, Ec );
 
  // Get the 4-momentum of the residue in the lab frame using conservation
  double Ed = pro_mom4.e() + tar_mom4.e() - eje_mom4.e();
  double pd_x = pro_mom4.px() + tar_mom4.px() - eje_mom4.px();
  double pd_y = pro_mom4.py() + tar_mom4.py() - eje_mom4.py();
  double pd_z = pro_mom4.pz() + tar_mom4.pz() - eje_mom4.pz();
 
  // Determine the residue mass from its 4-momentum
  md_ = marley_utils::real_sqrt( Ed*Ed - pd_x*pd_x - pd_y*pd_y - pd_z*pd_z );
 
  // The excitation energy is the mass difference between this mass and
  // the residue's ground-state mass
  double Ex = md_ - md_gs_;
 
  // Create particle objects representing the ejectile and residue
  auto ejectile = marley_hepmc3::make_particle( eje_mom4, pdg_c_,
    marley_hepmc3::NUHEPMC_FINAL_STATE_STATUS, mc_ );
 
  auto residue = marley_hepmc3::make_particle( pdg_d_, pd_x, pd_y, pd_z, Ed,
    marley_hepmc3::NUHEPMC_UNDECAYED_RESIDUE_STATUS, md_ );
 
  // Make the event object (this also sets the charge and nuclear level
  // attributes)
  auto event = this->make_nuclear_event_object( KEa, ejectile, residue, Ex,
    twoJ, P );
 
  return event;
}

References HepMC3::FourVector::e(), marley::NuclearReaction::KEa_threshold_, marley::Reaction::ma_, marley::NuclearReaction::make_nuclear_event_object(), marley::Reaction::mb_, marley::Reaction::mc_, marley::Reaction::md_, marley::NuclearReaction::md_gs_, marley::Reaction::pdg_a(), marley::Reaction::pdg_a_, marley::Reaction::pdg_c_, marley::Reaction::pdg_d_, HepMC3::FourVector::px(), HepMC3::FourVector::py(), HepMC3::FourVector::pz(), reassign_sub_continuum(), marley::Generator::sample_from_distribution(), marley::Generator::uniform_random_double(), and xsec_.

get_Ec_from_Ex()

double marley::ContinuumNuclearReaction::get_Ec_from_Ex ( const double Ex, const double cos_theta, const double KEa, double * jacobian = nullptr ) const

protected

Helper function for reassign_sub_continuum() that solves for the outgoing lepton total energy that corresponds to the input kinematic variables.

Parameters

[in]	Ex	Nuclear excitation energy (MeV)
[in]	cos_theta	Lepton scattering cosine
[in]	KEa	Projectile kinetic energy
[out]	jacobian	If this argument is not nullptr, then the target double will be filled with the value of the Jacobian needed to convert from \( d\sigma/dE_\ell \) to \( d\sigma/dE_x \).

Definition at line 366 of file ContinuumNuclearReaction.cc.

{
  // Mass of the final-state ion (including excitation energy)
  double md = md_gs_ + Ex;
 
  // Total energy of the projectile
  double Ea = KEa + ma_;
 
  // Total energy of the two-body system in the lab frame
  double Etot = Ea + mb_;
 
  // Construct helper variables
  double pa = marley_utils::real_sqrt( Ea*Ea - ma_*ma_ );
  double help = md*md - mc_*mc_ + pa*pa - Etot*Etot;
  double other_help = 4.*pa*pa*cos_theta*cos_theta;
 
  // Quadratic coefficients (a*Ec^2 + b*Ec + c == 0)
  double a = 4.*Etot*Etot - other_help;
  double b = 4.*Etot*help;
  double c = help*help + other_help*mc_*mc_;
 
  // Get both solutions to the quadratic equation
  double sol_plus, sol_minus;
  marley_utils::solve_quadratic_equation( a, b, c,
    sol_plus, sol_minus );
 
  // Now for a trick: due to the way we derived the results above,
  // the two solutions correspond to positive (sol_plus) and
  // negative (sol_minus) values of cos_theta, with the two
  // solutions exactly equal when cos_theta == 0. Choose the
  // appropriate one to return here.
  double Ec = sol_plus;
  if ( cos_theta < 0. ) Ec = sol_minus;
 
  if ( jacobian ) {
    double pc = marley_utils::real_sqrt( Ec*Ec - mc_*mc_ );
    *jacobian = md / ( Ea + mb_ - pa*Ec*cos_theta/pc );
  }
 
  return Ec;
}

References marley::Reaction::ma_, marley::Reaction::mb_, marley::Reaction::mc_, and marley::NuclearReaction::md_gs_.

Referenced by reassign_sub_continuum().

get_tabulated_xsec()

const TabulatedXSec & marley::ContinuumNuclearReaction::get_tabulated_xsec ( ) const

inline

Definition at line 49 of file ContinuumNuclearReaction.hh.

        { return *xsec_; }

reassign_sub_continuum()

bool marley::ContinuumNuclearReaction::reassign_sub_continuum ( double & w, const double ctl, const double KEa ) const

protected

Helper function for create_event() that potentially reassigns the value of the energy transfer.

If reassignment is enabled, moves cross-section strength falling below the unbound threshold to the continuum by updating the value of the energy transfer. Returns true if the reassignment was successful and is kinematically allowed. Returns false otherwise.

Parameters

[in,out]	w	Energy transfer (MeV, not shifted by \( \Delta_\mathrm{IAS} \))
[in]	ctl	Lepton scattering cosine
[in]	KEa	Projectile kinetic energy (MeV)

Returns

Returns true if the reassignment procedure was successful (or skipped because it was unnecessary or disabled). Returns false otherwise.

Definition at line 290 of file ContinuumNuclearReaction.cc.

{
  // If reassignment of the sub-continuum strength is disabled, then this
  // function returns immediately without doing anything
  if ( sc_mode_ == SubContinuumMode::IGNORE ) return true;
 
  // Lab-frame total energy and 3-momentum of the projectile
  double Ea = KEa + ma_;
  double pa = marley_utils::real_sqrt( Ea*Ea - ma_*ma_ );
 
  // Lab-frame total energy and 3-momentum of the ejectile
  double Ec = Ea - w;
  double pc = marley_utils::real_sqrt( Ec*Ec - mc_*mc_ );
 
  // Squared magnitude of the 3-momentum transfer
  double kappa2 = pa*pa + pc*pc - 2.*pa*pc*ctl;
 
  // Total energy and excitation energy of the residue
  double Ed = mb_ + w;
  double Ex = marley_utils::real_sqrt( Ed*Ed - kappa2 ) - md_gs_;
 
  // Get the "unbound threshold" used to determine the start of the continuum
  const auto& mt = marley::MassTable::Instance();
  double unbound_threshold = mt.unbound_threshold( pdg_d_ );
 
  // If the excitation energy corresponding to the sampled energy transfer
  // is already within the continuum, no special action is needed. Just return
  // without making any modifications.
  if ( Ex >= unbound_threshold ) return true;
 
  // If we've made it here, then we need to reassign the excitation energy
  // and compute a new value of the energy transfer. First choose the new
  // excitation energy based on the recipe selected by the user configuration.
  if ( sc_mode_ == SubContinuumMode::ACCUMULATE ) {
    // For the "accumulate" option, just update the excitation energy to be
    // exactly at the unbound threshold
    Ex = unbound_threshold;
 
    MARLEY_LOG( DEBUG, "physics.reaction" ) << "Excitation energy " << Ex
      << " MeV is below the unbound threshold " << unbound_threshold
      << " MeV. Sampling exactly at the unbound threshold.";
  }
  else if ( sc_mode_ == SubContinuumMode::MIRROR ) {
    // For the "mirror" option, "reflect" the original excitation energy to the
    // upper side of the unbound threshold so it is the same distance above
    // as it was originally below.
    Ex = 2.*unbound_threshold - Ex;
 
    MARLEY_LOG( DEBUG, "physics.reaction" ) << "Excitation energy " << Ex
      << " MeV is below the unbound threshold " << unbound_threshold
      << " MeV. Mirroring the energy transfer around the unbound"
      << " threshold.";
  }
  else {
    throw marley::Error( "Unrecognized sub-continuum mode encountered"
      " in marley::ContinuumNuclearReaction::reassign_sub_continuum()" );
    return false;
  }
 
  // Solve for the new outgoing lepton total energy given the updated
  // excitation energy value. Update the value of Ec with the solution.
  Ec = this->get_Ec_from_Ex( Ex, ctl, KEa );
 
  // Now update the energy transfer accordingly
  w = Ea - Ec;
 
  // If the total energy of the outgoing lepton is now below its rest mass,
  // then the reassignment procedure failed due to the kinematic threshold.
  // Indicate this failure in the return value.
  if ( Ec < mc_ ) return false;
 
  // Otherwise, everything worked out, so indicate success.
  return true;
}

References get_Ec_from_Ex(), marley::MassTable::Instance(), marley::Reaction::ma_, marley::Reaction::mb_, marley::Reaction::mc_, marley::NuclearReaction::md_gs_, marley::Reaction::pdg_d_, and sc_mode_.

Referenced by create_event().

set_description()

void marley::ContinuumNuclearReaction::set_description ( )

overrideprotectedvirtual

Creates the description string based on the PDG code values for the initial and final particles.

Reimplemented from marley::NuclearReaction.

Definition at line 281 of file ContinuumNuclearReaction.cc.

                                                     {
  marley::NuclearReaction::set_description();
  description_ += '*';
}

References marley::Reaction::description_, and marley::NuclearReaction::set_description().

set_sub_continuum_mode()

static void marley::ContinuumNuclearReaction::set_sub_continuum_mode ( SubContinuumMode scm )

inlinestatic

Sets the approach to handling sub-continuum cross-section strength.

Definition at line 56 of file ContinuumNuclearReaction.hh.

        { sc_mode_ = scm; }

References sc_mode_.

string_from_sub_continuum_mode()

std::string marley::ContinuumNuclearReaction::string_from_sub_continuum_mode ( const SubContinuumMode sc_mode )

static

Definition at line 421 of file ContinuumNuclearReaction.cc.

{
  auto it = sc_mode_string_map_.find( mode );
  if ( it != sc_mode_string_map_.end() ) return it->second;
  else throw marley::Error( "Unrecognized sub-continuum mode value encountered"
    " in marley::ContinuumNuclearReaction::string_from_sub_continuum_mode()" );
}

sub_continuum_mode()

static SubContinuumMode marley::ContinuumNuclearReaction::sub_continuum_mode ( )

inlinestatic

Gets the approach to handling sub-continuum cross-section strength.

Definition at line 53 of file ContinuumNuclearReaction.hh.

{ return sc_mode_; }

References sc_mode_.

sub_continuum_mode_from_string()

SubContinuumMode marley::ContinuumNuclearReaction::sub_continuum_mode_from_string ( const std::string & sc_mode_str )

static

Definition at line 410 of file ContinuumNuclearReaction.cc.

{
  for ( const auto& pair : sc_mode_string_map_ ) {
    if ( str == pair.second ) return pair.first;
  }
  throw marley::Error( "The string \"" + str + "\" was not recognized"
    " as a valid sub-continuum mode setting" );
}

total_xs()

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

overridevirtual

Compute the reaction's total cross section (MeV ^-2)

Parameters

pdg_a	Projectile's PDG code
KEa	Lab-frame kinetic energy of the incident projectile

Returns

Reaction total cross section (MeV ^-2)

Note

Functions that override total_xs() should always return zero if pdg_a != pdg_a_.

Implements marley::Reaction.

Definition at line 55 of file ContinuumNuclearReaction.cc.

{
  if ( pdg_a != pdg_a_ ) return 0.;
  return xsec_->integral( pdg_a, KEa );
}

References marley::Reaction::pdg_a(), marley::Reaction::pdg_a_, and xsec_.

Member Data Documentation

sc_mode_

SubContinuumMode marley::ContinuumNuclearReaction::sc_mode_ = SubContinuumMode::ACCUMULATE

staticprotected

Indicates the desired method for handling events with excitation energies originally sampled below the continuum threshold.

Definition at line 102 of file ContinuumNuclearReaction.hh.

Referenced by reassign_sub_continuum(), set_sub_continuum_mode(), and sub_continuum_mode().

sc_mode_string_map_

const std::map< SubContinuumMode, std::string > marley::ContinuumNuclearReaction::sc_mode_string_map_

staticprotected

Initial value:

= {
    { SubContinuumMode::IGNORE, "ignore" },
    { SubContinuumMode::MIRROR, "mirror" },
    { SubContinuumMode::ACCUMULATE, "accumulate" }
}

Helper map used for conversions between a SubContinuumMode value and a std::string.

Definition at line 39 of file ContinuumNuclearReaction.hh.

xsec_

std::shared_ptr< TabulatedXSec > marley::ContinuumNuclearReaction::xsec_

protected

Helper object that handles cross section calculations.

Definition at line 98 of file ContinuumNuclearReaction.hh.

Referenced by create_event(), and total_xs().

---

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

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