A neutrino-nucleus reaction whose cross section is calculated according to the allowed approximation. More...

#include <DiscreteNuclearReaction.hh>

Inheritance diagram for marley::DiscreteNuclearReaction:

Public Member Functions
	DiscreteNuclearReaction (ProcessType pt, int pdg_a, int pdg_b, int pdg_c, int pdg_d, int q_d, const std::shared_ptr< std::vector< marley::MatrixElement > > &mat_els, CoulombCorrector::CoulombMode mode, const JSON &ff_config, 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.

double	diff_xs (const marley::MatrixElement &mat_el, double KEa, double cos_theta_c_cm, double &beta_c_cm, bool check_max_E_level) const
	Differential cross section \(d\sigma/d\cos\theta_{c}^{\mathrm{CM}}\) (MeV^-2) evaluated in the center-of-momentum frame for a transition to a particular final nuclear level.

virtual double	diff_xs (int pdg_a, double KEa, double cos_theta_c_cm) const
	Differential cross section \(d\sigma/d\cos\theta_{c}^{\mathrm{CM}}\) (MeV^-2) summed over all kinematically-allowed final nuclear levels.

const std::vector< marley::MatrixElement > &	matrix_elements () const
	Allows access to the owned vector of MatrixElement objects.

void	set_decay_scheme (marley::DecayScheme *scheme)
	Sets the DecayScheme object to use for sampling excited levels in the residue.

virtual double	total_xs (const marley::MatrixElement &me, double KEa, double &beta_c_cm, bool check_max_E_level=true) const
	Total cross section (MeV^-2) for a given final nuclear level.

virtual double	total_xs (int pdg_a, double KEa) const override
	Total reaction cross section (MeV^-2) including all kinematically-allowed final nuclear levels.

virtual double	total_xs (int pdg_a, double KEa, marley::MatrixElement::TransitionType mat_el_type) const
	Total reaction cross section (MeV^-2) including all kinematically-allowed nuclear transitions of a particular type.

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.

Protected Member Functions
double	sample_cos_theta_c_cm (const marley::MatrixElement &matrix_el, double KEa, double beta_c_cm, marley::Generator &gen) const
	Samples a polar angle cosine for the ejectile using the relevant portion of the reaction nuclear matrix element.

virtual void	set_description () override
	Creates the description string based on the PDG code values for the initial and final particles.

double	summed_xs_helper (int pdg_a, double KEa, double cos_theta_c_cm, std::vector< double > *level_xsecs, bool differential) const

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
bool	allowed_approx_ = false
	Flag that indicates whether to include aditional terms beyond the q->0 limit.

CoulombCorrector	coulomb_corrector_
	Object that handles calculations of Coulomb correction factors.

std::shared_ptr< std::vector< marley::MatrixElement > >	matrix_elements_
	Matrix elements representing all of the possible nuclear transitions that may be caused by this reaction.

std::shared_ptr< NuclearFormFactor >	nuclear_ff_
	Object that handles calculations of the nuclear form factor.

NucleonFormFactors	nucleon_form_factors_
	Object that handles calculations of nucleon form factors.

std::pair< std::vector< int >, std::vector< double > >	nucleon_radii_
	Pair of vectors containg the mean radius of the nucleon wavefunctions and the degeneracies for each nuclear level.

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.

Additional Inherited Members
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...

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)

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

Detailed Description

A neutrino-nucleus reaction whose cross section is calculated according to the allowed approximation.

Definition at line 42 of file DiscreteNuclearReaction.hh.

Constructor & Destructor Documentation

◆ DiscreteNuclearReaction()

marley::DiscreteNuclearReaction::DiscreteNuclearReaction	(	ProcessType	pt,
		int	pdg_a,
		int	pdg_b,
		int	pdg_c,
		int	pdg_d,
		int	q_d,
		const std::shared_ptr< std::vector< marley::MatrixElement > > &	mat_els,
		CoulombCorrector::CoulombMode	mode,
		const JSON &	ff_config,
		const std::string &	source_file )

Parameters

pt	Type of scattering process represented by this Reaction
pdg_a	Projectile PDG code
pdg_b	Target PDG code
pdg_c	Ejectile PDG code
pdg_d	Residue PDG code
q_d	Charge of the residue after the prompt 2->2 interaction represented by this DiscreteNuclearReaction object
mat_els	A vector of MatrixElement objects that should be used to compute cross sections for this DiscreteNuclearReaction
mode	Indicates the method to use when computing Coulomb corrections for the reaction cross section
ff_config	JSON object containing information needed to configure the nucleon and nuclear form factor models

Definition at line 51 of file DiscreteNuclearReaction.cc.

  : marley::NuclearReaction( pt, pdg_a, pdg_b, pdg_c, pdg_d, q_d,
      source_file ),
  matrix_elements_( mat_els ), coulomb_corrector_( pdg_c, pdg_d, mode ),
  nucleon_form_factors_( ff_config )
{
  // Target nucleus proton number
  int Zb = marley_utils::get_particle_Z( pdg_b_ );
  // Target nucleus nucleon number
  int Ab = marley_utils::get_particle_A( pdg_b_ );
 
  // Choose the model to use for the nuclear form factor
  nuclear_ff_ = marley::NuclearFormFactor::create( Zb, Ab, ff_config );
 
  // Determine whether we're working within the complete allowed approximation
  // by checking the form factor JSON configuration. If the user has requested
  // this, then set the corresponding flag to true.
  allowed_approx_ = marley::JSONConfig
    ::check_for_allowed_approximation( ff_config );
}

References marley::NuclearReaction::NuclearReaction(), allowed_approx_, coulomb_corrector_, marley::NuclearFormFactor::create(), matrix_elements_, nuclear_ff_, nucleon_form_factors_, marley::Reaction::pdg_a(), marley::Reaction::pdg_b(), marley::Reaction::pdg_b_, and marley::Reaction::source_file().

Member Function Documentation

◆ create_event()

std::shared_ptr< HepMC3::GenEvent > marley::DiscreteNuclearReaction::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_.

Todo: Add more error checks to DiscreteNuclearReaction::create_event() as necessary

Todo: Revisit the equipartition assumption made here in favor of something better motivated.

Implements marley::Reaction.

Definition at line 78 of file DiscreteNuclearReaction.cc.

{
  // 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." );
 
  // Create an empty vector of sampling weights (partial total cross
  // sections to each kinematically accessible final level)
  std::vector<double> level_weights;
 
  // Create a discrete distribution object for level sampling.
  // Its default constructor creates a single weight of 1.
  // We will always explicitly give it weights to use when sampling
  // levels, so we won't worry about its default behavior.
  static std::discrete_distribution< size_t > ldist;
 
  // Compute the total cross section for a transition to each individual
  // nuclear level, and save the results in the level_weights vector (which
  // will be cleared by summed_xs_helper() before being loaded with the cross
  // sections). The summed_xs_helper() method can also be used for differential
  // (d\sigma/d\cos\theta_c^{CM}) cross sections, so supply a dummy
  // cos_theta_c_cm value and request total cross sections by setting the last
  // argument to false.
  double dummy = 0.;
  double sum_of_xsecs = summed_xs_helper( pdg_a, KEa, dummy,
    &level_weights, false );
 
  // Note that the elements in matrix_elements_ are given in order of
  // increasing excitation energy (this is currently enforced by the reaction
  // data format and is checked during parsing). This ensures that we can
  // sample a matrix element index from level_weights (which is populated in
  // the same order by summed_xs_helper()) and have it refer to the correct
  // object.
 
  // Complain if none of the levels we have data for are kinematically allowed
  if ( level_weights.empty() ) {
    throw marley::Error( "Could not create this event. The DecayScheme object"
      " associated with this reaction does not contain data for any"
      " kinematically accessible levels for a projectile kinetic energy of "
      + std::to_string( KEa ) + " MeV (max E_level = "
      + std::to_string( max_level_energy(KEa) ) + " MeV)." );
  }
 
  // Complain if the total cross section (the sum of all partial level 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 kinematically"
      " accessible levels for a projectile kinetic energy of "
      + std::to_string( KEa ) + " MeV (max E_level = "
      + std::to_string( max_level_energy(KEa) )
      + " MeV) have vanishing matrix elements." );
  }
 
  // Create a list of parameters used to supply the weights to our discrete
  // level sampling distribution
  std::discrete_distribution<size_t>::param_type params( level_weights.begin(),
    level_weights.end() );
 
  // Sample a matrix_element using our discrete distribution and the
  // current set of weights
  size_t me_index = gen.sample_from_distribution( ldist, params );
 
  const auto& sampled_matrix_el = matrix_elements_->at( me_index );
 
  // Get the energy of the selected level.
  double E_level = sampled_matrix_el.level_energy();
 
  // Update the residue mass based on its excitation energy for the current
  // event
  md_ = md_gs_ + E_level;
 
  // Compute Mandelstam s, the ejectile's CM frame total energy, the magnitude
  // of the ejectile's CM frame 3-momentum, and the residue's CM frame total
  // energy.
  double s, Ec_cm, pc_cm, Ed_cm;
  two_two_scatter( KEa, s, Ec_cm, pc_cm, Ed_cm );
 
  // Determine the CM frame velocity of the ejectile
  double beta_c_cm = pc_cm / Ec_cm;
 
  // Sample a CM frame scattering cosine for the ejectile.
  double cos_theta_c_cm = sample_cos_theta_c_cm( sampled_matrix_el, KEa,
    beta_c_cm, gen );
 
  // Sample a CM frame azimuthal scattering angle (phi) uniformly on [0, 2*pi).
  // We can do this because the matrix elements are azimuthally invariant
  double phi_c_cm = 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.
  int twoJ = BOGUS_TWO_J_VALUE;
  marley::Parity P; // defaults to positive parity
 
  // Get access to the nuclear structure database owned by the Generator
  auto& sdb = gen.get_structure_db();
 
  // Retrieve the ground-state spin-parity of the initial nucleus
  int twoJ_gs;
  marley::Parity P_gs;
 
  sdb.get_gs_spin_parity( pdg_b_, twoJ_gs, P_gs );
 
  // For transitions to discrete nuclear levels, all we need to do is retrieve
  // these values directly from the Level object
  const marley::Level* final_lev = sampled_matrix_el.level();
  if ( final_lev ) {
    twoJ = final_lev->twoJ();
    P = final_lev->parity();
  }
  // For transitions to the continuum, we rely on the spin-parity selection
  // rules to determine suitable values of twoJ and P. In cases where more than
  // one value is allowed, assume equipartition, and sample a spin-parity based
  // on the relative nuclear level densities at the excitation energy of
  // interest.
  else {
 
    // For a Fermi transition, the final spin-parity is always the same as the
    // initial one
    if ( sampled_matrix_el.type() == ME_Type::FERMI ) {
      twoJ = twoJ_gs;
      P = P_gs;
    }
    else if ( sampled_matrix_el.type() == ME_Type::GAMOW_TELLER ) {
      // For a Gamow-Teller transition, the final parity is the same as the
      // initial parity
      P = P_gs;
 
      // For a spin-zero initial state, take a shortcut: the final spin will
      // always be one.
      if ( twoJ_gs == 0 ) twoJ = 2;
      else {
 
        // For an initial state with a non-zero spin, make a vector storing
        // all of the spin values allowed by the GT selection rules.
        // Sample an allowed value assuming equipartition of spin. Use the
        // relative final nuclear level densities as sampling weights.
        std::vector<int> allowed_twoJs;
        std::vector<double> ld_weights;
 
        auto& ldm = sdb.get_level_density_model( pdg_d_ );
 
        for ( int myTwoJ = std::abs(twoJ_gs - 2); myTwoJ <= twoJ_gs + 2;
          myTwoJ += 2 )
        {
          allowed_twoJs.push_back( myTwoJ );
          ld_weights.push_back( ldm.level_density(E_level, myTwoJ, P) );
        }
 
        std::discrete_distribution<size_t> my_twoJ_dist( ld_weights.begin(),
          ld_weights.end() );
 
        size_t my_index = gen.sample_from_distribution( my_twoJ_dist );
        twoJ = allowed_twoJs.at( my_index );
      }
    }
    else throw marley::Error( "Unrecognized matrix element type encountered"
      " in marley::DiscreteNuclearReaction::create_event()" );
  }
 
  MARLEY_LOG( DEBUG, "physics.reaction" ) << "Sampled a "
    << sampled_matrix_el.type_str()
    << " transition from the " << marley::TargetAtom( pdg_b_ )
    << " ground state (with spin-parity " << static_cast<double>( twoJ_gs ) / 2.
    << P_gs << ") to the " << marley::TargetAtom( pdg_d_ )
    << " level with Ex = " << E_level << " MeV and spin-parity "
    << static_cast<double>( twoJ ) / 2. << P;
 
  // Create the preliminary event object (after 2-->2 scattering, but before
  // de-excitation of the residual nucleus)
  auto event = this->make_nuclear_event_object( KEa, pc_cm, cos_theta_c_cm,
    phi_c_cm, Ec_cm, Ed_cm, E_level, twoJ, P );
 
  // Store the matrix element index for use by weight calculators
  event->add_attribute( "me_index",
    std::make_shared< HepMC3::IntAttribute >(
      static_cast<int>( me_index ) ) );
 
  // Return the preliminary event object (to be processed later by the
  // NucleusDecayer class)
  return event;
}

References marley::Generator::get_structure_db(), marley::NuclearReaction::KEa_threshold_, marley::NuclearReaction::make_nuclear_event_object(), matrix_elements_, marley::NuclearReaction::max_level_energy(), marley::Reaction::md_, marley::NuclearReaction::md_gs_, marley::Level::parity(), marley::Reaction::pdg_a(), marley::Reaction::pdg_a_, marley::Reaction::pdg_b_, marley::Reaction::pdg_d_, sample_cos_theta_c_cm(), marley::Generator::sample_from_distribution(), summed_xs_helper(), marley::Reaction::two_two_scatter(), marley::Level::twoJ(), and marley::Generator::uniform_random_double().

◆ diff_xs() [1/2]

double marley::DiscreteNuclearReaction::diff_xs	(	const marley::MatrixElement &	mat_el,
		double	KEa,
		double	cos_theta_c_cm,
		double &	beta_c_cm,
		bool	check_max_E_level ) const

Differential cross section \(d\sigma/d\cos\theta_{c}^{\mathrm{CM}}\) (MeV^-2) evaluated in the center-of-momentum frame for a transition to a particular final nuclear level.

Parameters

mat_el	MatrixElement object describing the transition to the final nuclear level
KEa	Lab-frame projectile kinetic energy (MeV)
cos_theta_c_cm	Ejectile scattering cosine as measured in the CM frame

Definition at line 324 of file DiscreteNuclearReaction.cc.

{
  // Check that the scattering cosine is within the physically meaningful range
  if ( std::abs(cos_theta_c_cm) > 1. ) return 0.;
 
  // Don't bother to compute anything if the matrix element vanishes
  if ( mat_el.strength() == 0. ) return 0.;
 
  // Also don't proceed further if the reaction is below threshold
  // (equivalently, if the requested level excitation energy E_level exceeds
  // that maximum kinematically-allowed value). To avoid redundant checks of
  // the threshold, skip this check if check_max_E_level is set to false.
  if ( check_max_E_level ) {
    double max_E_level = max_level_energy( KEa );
    if ( mat_el.level_energy() > max_E_level ) return 0.;
  }
 
  // The final nuclear mass (before nuclear de-excitations) is the sum of the
  // ground state residue mass plus the excitation energy of the accessed level
  // This includes a correction to the ground-state atomic mass to account for
  // production of an ion by charged-current interactions (see the constructor
  // of the NuclearReaction class)
  double md2 = std::pow( md_gs_ + mat_el.level_energy(), 2 );
 
  // Compute Mandelstam s (the square of the total CM frame energy)
  double s = std::pow( ma_ + mb_, 2 ) + 2.*mb_*KEa;
  double sqrt_s = std::sqrt( s );
 
  // Compute some CM frame total energies and 3-momentum magnitudes
  double Eb_cm = ( s + mb_*mb_ - ma_*ma_ ) / ( 2. * sqrt_s );
 
  double Ea_cm = sqrt_s - Eb_cm;
  double pa_cm = marley_utils::real_sqrt( std::pow(Ea_cm, 2) - ma_*ma_ );
 
  double Ec_cm = ( s + mc_*mc_ - md2 ) / ( 2. * sqrt_s );
  double pc_cm = marley_utils::real_sqrt( std::pow(Ec_cm, 2) - mc_*mc_ );
 
  // Compute the CM frame value of the energy transfer (omega)
  double omega_cm = Ea_cm - Ec_cm;
 
  // Compute the CM frame value of the 3-momentum transfer magnitude (kappa)
  double kappa_cm = marley_utils::real_sqrt( std::pow( pa_cm, 2 )
    + std::pow( pc_cm, 2 ) - 2. * pa_cm * pc_cm * cos_theta_c_cm );
 
  // Negative square of the 4-momentum transfer
  double Q2 = kappa_cm*kappa_cm - omega_cm*omega_cm;
 
  // Compute the (dimensionless) speed of the ejectile in the CM frame
  beta_c_cm = pc_cm / Ec_cm;
 
  // CM frame total energy of the nuclear residue
  double Ed_cm = sqrt_s - Ec_cm;
 
  // Common factors for both CC and NC differential cross sections to a discrete
  // nuclear level in the CM frame
  double diff_xsec_prefactor = ( marley_utils::GF2 / ( 2 * marley_utils::pi ) )
    * ( Eb_cm * Ed_cm / s ) * Ec_cm * pc_cm;
 
  // Apply extra factors based on the current process type
  if ( process_type_ == ProcessType::NeutrinoCC_Discrete
    || process_type_ == ProcessType::AntiNeutrinoCC_Discrete )
  {
    // Dot product of the four-momenta of particles c and d
    double pc_dot_pd = Ed_cm*Ec_cm + std::pow( pc_cm, 2 );
 
    // Relative speed of particles c and d, computed with a manifestly
    // Lorentz-invariant expression
    double beta_rel_cd = marley_utils::real_sqrt(
      std::pow(pc_dot_pd, 2) - mc_*mc_*md2 ) / pc_dot_pd;
 
    // Calculate a Coulomb correction factor using either a Fermi function
    // or the effective momentum approximation
    double factor_C = coulomb_corrector_.coulomb_correction_factor(
      beta_rel_cd );
    diff_xsec_prefactor *= marley_utils::Vud2 * factor_C;
  }
  else if ( process_type_ == ProcessType::NC_Discrete )
  {
    // For NC, extra factors are only needed for Fermi transitions (which
    // correspond to CEvNS since they can only access the nuclear ground state)
    if ( mat_el.type() == ME_Type::FERMI ) {
      double Q_w = weak_nuclear_charge();
      diff_xsec_prefactor *= 0.25*std::pow( Q_w, 2 );
    }
  }
  else throw marley::Error( "Unrecognized or invalid process type encountered"
    " in marley::DiscreteNuclearReaction::diff_xs()" );
 
  // We're done with the overall factors. Now compute the lepton part of the
  // tensor contraction
  // @todo Reduce code duplication with a similar implementation of these (in
  // the lab frame) within the TabulatedXSec class. Note that there are
  // different conventions used there.
  int helicity = marley_utils::get_particle_helicity( pdg_a_ );
  double sin2_theta_c_cm = 1. - cos_theta_c_cm*cos_theta_c_cm;
  double vcc = 1. + beta_c_cm * cos_theta_c_cm;
  double vll = vcc - 2.*Ea_cm*Ec_cm*sin2_theta_c_cm
    * beta_c_cm*beta_c_cm/kappa_cm/kappa_cm;
  double vcl = -1.* ( omega_cm*vcc/kappa_cm + mc_*mc_/Ec_cm/kappa_cm );
  double vT = 1. - beta_c_cm*cos_theta_c_cm + Ea_cm*Ec_cm
    * beta_c_cm*beta_c_cm*sin2_theta_c_cm/kappa_cm/kappa_cm;
  double vTprime = helicity * ( (Ea_cm + Ec_cm)
    * (1. - beta_c_cm*cos_theta_c_cm)/kappa_cm - mc_*mc_/kappa_cm/Ec_cm );
  marley::LeptonFactors lf( vcc, vll, vcl, vT, vTprime );
 
  // Now compute the nuclear responses. Copy the energy transfer, 3-momentum
  // transfer magnitude, and Q^2 to effective values. These will be set to zero
  // if we're working in the allowed approximation.
  double omega_cm_eff = allowed_approx_ ? 0. : omega_cm;
  double kappa_cm_eff = allowed_approx_ ? 0. : kappa_cm;
  double Q2_eff = allowed_approx_ ? 0. : Q2;
 
  // Scale the transition strength by the squared nuclear form factor
  double strength_eff = mat_el.strength();
  strength_eff *= std::pow( nuclear_ff_->F( kappa_cm_eff ), 2 );
 
  // Scale the strength by the relevant nucleon form factor and divide by the
  // relevant coupling constant. Use the scaled value to compute the nuclear
  // responses relevant to the chosen transition type
  double rCC, rCL, rLL, rTvv, rTaa, rTprime;
  const double kM = kappa_cm_eff / marley_utils::m_nucleon;
  const double k2M = kappa_cm_eff * kappa_cm_eff / marley_utils::m_nucleon;
 
  if ( mat_el.type() == ME_Type::FERMI ) {
    double F1 = nucleon_form_factors_.F1( Q2_eff );
    strength_eff *= F1*F1 / marley_utils::g_V2;
 
    rCC = strength_eff;
    rCL = strength_eff * kM;
    rLL = strength_eff * kM * kM / 4.;
    rTvv = 0.;
    rTaa = 0.;
    rTprime = 0.;
  }
  else if ( mat_el.type() == ME_Type::GAMOW_TELLER ) {
    double FA = nucleon_form_factors_.FA( Q2_eff );
    strength_eff *= FA*FA / marley_utils::g_A2;
 
    double FP = nucleon_form_factors_.FP( Q2_eff );
    double F1 = nucleon_form_factors_.F1( Q2_eff );
    double F2 = nucleon_form_factors_.F2( Q2_eff );
 
    double FPA = FP / FA;
    double F12A = ( F1 + 2. * marley_utils::m_nucleon * F2 ) / FA;
 
    rCC = strength_eff * kM * kM / 12. * ( 1. - 2.*omega_cm_eff*FPA
      + omega_cm_eff*omega_cm_eff*FPA*FPA );
    rCL = strength_eff * kM / 3. * ( 1. - (omega_cm_eff + 0.5*k2M)*FPA
      + 0.5*omega_cm_eff*k2M*FPA*FPA );
    rLL = strength_eff * ( 1./3. - k2M/3.*FPA + k2M*k2M/12.*FPA*FPA );
    rTvv = 0.;
    rTaa = strength_eff * ( 2./3. + kM*kM/6.*F12A*F12A );
    rTprime = -1. * strength_eff * 2./3. * kM * F12A;
  }
  else throw marley::Error( "Unrecognized matrix element type encountered in"
    " marley::DiscreteNuclearReaction::diff_xs()" );
 
  marley::NuclearResponses nr( rCC, rLL, rCL, rTvv, rTaa, rTprime );
 
  double diff_xsec = diff_xsec_prefactor * ( lf * nr );
  return diff_xsec;
}

References allowed_approx_, marley::Reaction::AntiNeutrinoCC_Discrete, coulomb_corrector_, marley::MatrixElement::level_energy(), marley::Reaction::ma_, marley::NuclearReaction::max_level_energy(), marley::Reaction::mb_, marley::Reaction::mc_, marley::NuclearReaction::md_gs_, marley::Reaction::NC_Discrete, marley::Reaction::NeutrinoCC_Discrete, nuclear_ff_, nucleon_form_factors_, marley::Reaction::pdg_a_, marley::Reaction::process_type_, marley::MatrixElement::strength(), marley::MatrixElement::type(), and marley::NuclearReaction::weak_nuclear_charge().

Referenced by sample_cos_theta_c_cm(), summed_xs_helper(), and total_xs().

◆ diff_xs() [2/2]

double marley::DiscreteNuclearReaction::diff_xs	(	int	pdg_a,
		double	KEa,
		double	cos_theta_c_cm ) const

virtual

Differential cross section \(d\sigma/d\cos\theta_{c}^{\mathrm{CM}}\) (MeV^-2) summed over all kinematically-allowed final nuclear levels.

Parameters

pdg_a	PDG code for the projectile
KEa	Lab-frame kinetic energy (MeV) of the projectile
cos_theta_c_cm	CM frame scattering cosine of the ejectile

Note: This function returns 0. if pdg_a != pdg_a_.

Definition at line 315 of file DiscreteNuclearReaction.cc.

{
  return summed_xs_helper( pdg_a, KEa, cos_theta_c_cm, nullptr, true );
}

References marley::Reaction::pdg_a(), and summed_xs_helper().

◆ matrix_elements()

const std::vector< marley::MatrixElement > & marley::DiscreteNuclearReaction::matrix_elements ( ) const

inline

Allows access to the owned vector of MatrixElement objects.

Definition at line 130 of file DiscreteNuclearReaction.hh.

131 { return *matrix_elements_; }

References matrix_elements_.

Referenced by marley::StrengthVariationWeightCalculator::create_instances().

◆ sample_cos_theta_c_cm()

double marley::DiscreteNuclearReaction::sample_cos_theta_c_cm	(	const marley::MatrixElement &	matrix_el,
		double	KEa,
		double	beta_c_cm,
		marley::Generator &	gen ) const

protected

Samples a polar angle cosine for the ejectile using the relevant portion of the reaction nuclear matrix element.

Parameters

m_type	Integer representing the type of transition (0 = Fermi, 1 = Gamow-Teller)
beta_c_cm	Dimensionless speed of the ejectile in the CM frame
gen	Reference to the Generator to use for random sampling

Definition at line 605 of file DiscreteNuclearReaction.cc.

{
  // For now the max is unknown, so the rejection sampling algorithm will
  // calculate it "on the fly"
  double max = marley_utils::UNKNOWN_MAX;
 
  // For the allowed approximation, we know where it will be a priori
  if ( allowed_approx_ ) {
    if ( mat_el.type() == ME_Type::FERMI ) {
      max = this->diff_xs( mat_el, KEa, 1., beta_c_cm, false );
    }
    else if ( mat_el.type() == ME_Type::GAMOW_TELLER ) {
      max = this->diff_xs( mat_el, KEa, -1., beta_c_cm, false );
    }
    else throw marley::Error( "Unrecognized matrix element type encountered"
      " in marley::DiscreteNuclearReaction::sample_cos_theta_c_cm()" );
  }
 
  return gen.rejection_sample(
    [ &mat_el, KEa, &beta_c_cm, this ]( double cos_theta_cm ) -> double
    { return this->diff_xs( mat_el, KEa, cos_theta_cm, beta_c_cm, false ); },
    -1., 1., max );
}

References allowed_approx_, diff_xs(), marley::Generator::rejection_sample(), and marley::MatrixElement::type().

Referenced by create_event().

◆ set_description()

void marley::DiscreteNuclearReaction::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 633 of file DiscreteNuclearReaction.cc.

                                                    {
  // Initialize the description_ member variable with the basic information
  // first
  marley::NuclearReaction::set_description();
  // Now decide whether this reaction can access excited levels in the daughter
  // nucleus
  bool has_excited_state = false;
  for ( const auto& me : *matrix_elements_ ) {
    if ( me.level_energy() > 0. ) {
      has_excited_state = true;
      break;
    }
  }
  if ( has_excited_state ) description_ += '*';
  else description_ += " (g.s.)";
}

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

◆ summed_xs_helper()

double marley::DiscreteNuclearReaction::summed_xs_helper	(	int	pdg_a,
		double	KEa,
		double	cos_theta_c_cm,
		std::vector< double > *	level_xsecs,
		bool	differential ) const

protected

Helper function for total_xs and summed_diff_xs()

Parameters

pdg_a	PDG code for the projectile
KEa	Lab-frame projectile kinetic energy (MeV)
cos_theta_c_cm	CM frame scattering cosine for the ejectile
level_xsecs	If this pointer is not nullptr, then the std::vector<double> that it points to will be cleared and loaded with the partial cross sections to each individual nuclear level included in the sum. This feature is helpful when computing weights for sampling a nuclear transition in create_event().
differential	Whether the total cross section (false) or the differential cross section (true, \(d\sigma/d\cos\theta_{c}^{\mathrm{CM}}\)) should be summed over all kinematically accessible nuclear levels by this function

Returns: The requested cross section (MeV^-2)

Todo: Consider whether you should use an exception if pdg_a != pdg_a_

Definition at line 520 of file DiscreteNuclearReaction.cc.

{
  // Check that the projectile supplied to this event is correct. If not,
  // return a total cross section of zero since this reaction is not available
  // for the given projectile.
  if ( pdg_a != pdg_a_ ) return 0.;
 
  // If we're evaluating a differential cross section, check that the
  // scattering cosine is within the physically meaningful range. If it's
  // not, then just return 0.
  if ( differential && std::abs(cos_theta_c_cm) > 1. ) return 0.;
 
  // If the projectile kinetic energy is zero (or negative), then
  // just return zero.
  if ( KEa <= 0. ) return 0.;
 
  // If we've been passed a vector to load with the partial cross sections
  // to each nuclear level, then clear it before storing them
  if ( level_xsecs ) level_xsecs->clear();
 
  double max_E_level = max_level_energy( KEa );
  double xsec = 0.;
  for ( const auto& mat_el : *matrix_elements_ ) {
 
    // Get the excitation energy for the current level
    double level_energy = mat_el.level_energy();
 
    // Exit the loop early if you reach a level with an energy that's too high
    if ( level_energy > max_E_level ) break;
 
    // Check whether the matrix element (B(F) + B(GT)) is nonvanishing for the
    // current level. If it is, just set the weight equal to zero rather than
    // computing the total xs.
    double partial_xsec = 0.;
    if ( mat_el.strength() != 0. ) {
 
      // Set the check_max_E_level flag to false when calculating the total
      // cross section for this level (we've already verified that the
      // current level is kinematically accessible in the check against
      // max_E_level above)
      double beta_c_cm = 0.;
 
      // Compute either the total or differential (d\sigma / d\cos\theta_{CM})
      // cross section as requested
      if ( differential ) {
        partial_xsec = diff_xs( mat_el, KEa, cos_theta_c_cm, beta_c_cm, false );
      } else {
        partial_xsec = total_xs( mat_el, KEa, beta_c_cm, false );
      }
 
      if ( std::isnan(partial_xsec) ) {
        MARLEY_LOG( WARN, "physics.reaction.xsec" )
          << "Partial cross section for reaction "
          << description_ << " gave NaN result.";
        MARLEY_LOG( DEBUG, "physics.reaction.xsec" )
          << "Parameters were level energy = "
          << mat_el.level_energy() << " MeV, projectile kinetic energy = "
          << KEa << " MeV, and reduced matrix element = " << mat_el.strength()
          << ". Differential was set to " << differential << ".";
        MARLEY_LOG( DEBUG, "physics.reaction.xsec" )
          << "The partial cross section to this level"
          << " will be set to zero.";
        partial_xsec = 0.;
      }
 
      xsec += partial_xsec;
 
    } // mat_el.strength() != 0.
 
    // Store the partial cross section to the current individual nuclear
    // level if needed (i.e., if level_xsecs is not nullptr). This is
    // done even when the current matrix element is zero so that the
    // partial cross sections for each transition match the ordering
    // in the vector of matrix elements.
    if ( level_xsecs ) level_xsecs->push_back( partial_xsec );
 
  } // loop over matrix elements
 
  return xsec;
}

References marley::Reaction::description_, diff_xs(), matrix_elements_, marley::NuclearReaction::max_level_energy(), marley::Reaction::pdg_a(), marley::Reaction::pdg_a_, and total_xs().

Referenced by create_event(), diff_xs(), total_xs(), and total_xs().

◆ total_xs() [1/3]

double marley::DiscreteNuclearReaction::total_xs	(	const marley::MatrixElement &	me,
		double	KEa,
		double &	beta_c_cm,
		bool	check_max_E_level = true ) const

virtual

Total cross section (MeV^-2) for a given final nuclear level.

Parameters

	me	MatrixElement object describing the transition to the final nuclear level
	KEa	Lab-frame projectile kinetic energy (MeV)
[out]	beta_c_cm	After this function is called, this variable will contain the (dimensionless) speed of the ejectile as measured in the center of momentum frame
	check_max_E_level	Whether to check the reaction threshold via a call to max_level_energy(). If this argument is set to false, then the check will be skipped.

Definition at line 491 of file DiscreteNuclearReaction.cc.

{
  if ( allowed_approx_ ) return 2. * this->diff_xs( mat_el, KEa, 0., beta_c_cm,
    check_max_E_level );
 
  // Integrator object to integrate over the scattering angle
  static marley::Integrator integrator;
 
  // Don't bother to compute anything if the matrix element vanishes
  if ( mat_el.strength() == 0. ) return 0.;
 
  // Integrate the differential cross section over cos_theta_c_cm
  double total_xsec = integrator.num_integrate(
    [ &mat_el, KEa, &beta_c_cm, check_max_E_level, this ](
      double cos_theta_cm ) -> double
    { return this->diff_xs( mat_el, KEa, cos_theta_cm, beta_c_cm,
      check_max_E_level ); }, -1., 1.
  );
 
  MARLEY_LOG( DEBUG, "physics.reaction.xsec" ) << "total xsec " << description_
    << " to level with energy " << mat_el.level_energy() << " MeV is "
    << total_xsec << " MeV^(-2).";
 
  return total_xsec;
}

References allowed_approx_, marley::Reaction::description_, diff_xs(), marley::MatrixElement::level_energy(), marley::Integrator::num_integrate(), and marley::MatrixElement::strength().

◆ total_xs() [2/3]

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

overridevirtual

Total reaction cross section (MeV^-2) including all kinematically-allowed final nuclear levels.

Parameters

pdg_a	PDG code for the projectile
KEa	Lab-frame kinetic energy (MeV) of the projectile

Returns: Reaction total cross section (MeV^-2)

Note: This function returns 0. if pdg_a != pdg_a_.

Implements marley::Reaction.

Definition at line 281 of file DiscreteNuclearReaction.cc.

{
  double dummy_cos_theta = 0.;
  return summed_xs_helper( pdg_a, KEa, dummy_cos_theta, nullptr, false );
}

References marley::Reaction::pdg_a(), and summed_xs_helper().

Referenced by summed_xs_helper().

◆ total_xs() [3/3]

double marley::DiscreteNuclearReaction::total_xs	(	int	pdg_a,
		double	KEa,
		marley::MatrixElement::TransitionType	mat_el_type ) const

virtual

Total reaction cross section (MeV^-2) including all kinematically-allowed nuclear transitions of a particular type.

Parameters

pdg_a	PDG code for the projectile
KEa	Lab-frame kinetic energy (MeV) of the projectile
mat_el_type	marley::MatrixElement::TransitionType for the kind of nuclear transition of interest

Returns: Reaction total cross section (MeV^-2)

Note: This function returns 0. if pdg_a != pdg_a_.

Definition at line 291 of file DiscreteNuclearReaction.cc.

{
  double dummy_cos_theta = 0.;
  std::vector< double > level_xsecs;
  // Accumulate partial cross section values for all of the accessible
  // nuclear transitions
  this->summed_xs_helper( pdg_a, KEa, dummy_cos_theta, &level_xsecs, false );
 
  // Sum up the contributions from transitions that match the requested
  // matrix element type, then return the result
  double xsec = 0.;
  for ( size_t j = 0u; j < level_xsecs.size(); ++j ) {
    const auto& ml = matrix_elements_->at( j );
    if ( ml.type() == mat_el_type ) {
      xsec += level_xsecs.at( j );
    }
  }
  return xsec;
}