marley::DecayScheme Class Reference

Discrete level and γ-ray data for a specific nuclide. More...

#include <DecayScheme.hh>

Public Types
enum class	FileFormat { native , talys }
	The FileFormat type is used to tell the DecayScheme class which format to assume when parsing a discrete level data file. More...

Public Member Functions
	DecayScheme (int Z, int A)
	Create a DecayScheme without any levels.
	DecayScheme (int Z, int A, const std::string &filename, FileFormat format=FileFormat::talys)
	Create a DecayScheme using discrete level data from a file.
int	A () const
	Get the mass number.
marley::Level &	add_level (const marley::Level &level)
	Add a level to the DecayScheme.
void	do_cascade (marley::Level &initial_level, HepMC3::GenEvent &event, marley::Generator &gen, std::shared_ptr< HepMC3::GenParticle > &residue)
	Simulates nuclear de-excitation via γ-ray emission(s)
const std::vector< std::unique_ptr< marley::Level > > &	get_levels () const
	Get a const reference to the vector that holds the Level objects.
marley::Level *	get_pointer_to_closest_level (double E_level)
	Gets a pointer to the Level in the DecayScheme whose excitation energy is closest to E_level.
int	pdg () const
	Returns the nuclear PDG code corresponding to Z and A.
void	print (std::ostream &out=std::cout) const
	Print this DecayScheme object to a std::ostream.
void	print_latex_table (std::ostream &ostr=std::cout)
	Print LaTeX source code that gives a tabular representation of the DecayScheme object.
void	print_report (std::ostream &ostr=std::cout) const
	Print a human-readable text representation of the DecayScheme object.
void	read_from_stream (std::istream &in)
	Use a std::istream to initialize this DecayScheme object, replacing any previous data.
int	Z () const
	Get the atomic number.

Protected Member Functions
size_t	level_lower_bound_index (double Ex)
	Get the index of the first level whose energy is not less than Ex.

Protected Attributes
int	A_
	Mass number.
std::vector< std::unique_ptr< marley::Level > >	levels_
	Level objects owned by this DecayScheme.
int	Z_
	Atomic number.

Detailed Description

Discrete level and γ-ray data for a specific nuclide.

Definition at line 34 of file DecayScheme.hh.

Member Enumeration Documentation

FileFormat

enum class marley::DecayScheme::FileFormat

strong

The FileFormat type is used to tell the DecayScheme class which format to assume when parsing a discrete level data file.

Currently, only discrete level data in MARLEY's native format (which may be read from a std::istream to initialize a marley::DecayScheme object via the >> operator) or in the format used by the TALYS nuclear code are allowed.

Definition at line 45 of file DecayScheme.hh.

{ native, talys };

Constructor & Destructor Documentation

DecayScheme() [1/3]

marley::DecayScheme::DecayScheme ( )

inline

Definition at line 47 of file DecayScheme.hh.

: Z_(0), A_(0) {}

DecayScheme() [2/3]

marley::DecayScheme::DecayScheme	(	int	Z,
		int	A )

Create a DecayScheme without any levels.

Parameters

Z	Atomic number of the desired nuclide
A	Mass number of the desired nuclide

Definition at line 221 of file DecayScheme.cc.

                                             : Z_( Z ), A_( A )
{
}

References A(), A_, Z(), and Z_.

DecayScheme() [3/3]

marley::DecayScheme::DecayScheme	(	int	Z,
		int	A,
		const std::string &	filename,
		FileFormat	format = FileFormat::talys )

Create a DecayScheme using discrete level data from a file.

Parameters

Z	Atomic number of the desired nuclide
A	Mass number of the desired nuclide
filename	Name of the file containing the discrete level data
format	FileFormat specifier that indicates which nuclear data format is used in the file

Definition at line 225 of file DecayScheme.cc.

                                   : Z_( Z ), A_( A )
{
  this->parse( filename, ff );
}

References A(), A_, Z(), and Z_.

Member Function Documentation

A()

int marley::DecayScheme::A ( ) const

inline

Get the mass number.

Definition at line 155 of file DecayScheme.hh.

{ return A_; }

References A_.

Referenced by DecayScheme(), and DecayScheme().

add_level()

marley::Level & marley::DecayScheme::add_level

(

const marley::Level &

level

)

Add a level to the DecayScheme.

Returns

A reference to the newly-added Level object

Definition at line 444 of file DecayScheme.cc.

{
  // Compute the numerical index for where we will insert the new level
  size_t index = level_lower_bound_index( level.energy() );
 
  // Insert the new level into the decay scheme
  levels_.insert( levels_.begin() + index,
    std::make_unique< marley::Level >(level) );
 
  // Return a reference to the newly-added level
  return *levels_.at( index );
}

References marley::Level::energy(), level_lower_bound_index(), and levels_.

Referenced by read_from_stream().

do_cascade()

void marley::DecayScheme::do_cascade	(	marley::Level &	initial_level,
		HepMC3::GenEvent &	event,
		marley::Generator &	gen,
		std::shared_ptr< HepMC3::GenParticle > &	residue )

Simulates nuclear de-excitation via γ-ray emission(s)

Gamma-rays will be randomly emitted until the nucleus reaches its ground state.

Parameters

[in]	initial_level	Reference to the first level that will de-excite via γ-ray emission
[in,out]	event	Reference to an Event object that will store the emitted γs
[in]	gen	Reference to the Generator to use for random sampling
	qIon	Net charge of the atom or ion whose nucleus is de-exciting

Definition at line 80 of file DecayScheme.cc.

{
  int qIon = marley_hepmc3::get_particle_charge( *residue );
 
  MARLEY_LOG( DEBUG, "physics.deexcitation.gamma" )
    << "Beginning gamma cascade at level with energy "
    << initial_level.energy() << " MeV";
 
  marley::Level* p_current_level = &initial_level;
 
  const marley::MassTable& mt = marley::MassTable::Instance();
 
  // Initialize some variables used in the gamma cascade loop below
  bool cascade_finished = false;
  double gamma_branching_ratio = 0.;
  double level_total_width = 0.;
 
  while ( !cascade_finished ) {
 
    // Randomly select a gamma to produce while storing its branching ratio
    const marley::Gamma* p_gamma = p_current_level->sample_gamma( gen,
      &gamma_branching_ratio );
 
    if ( !p_gamma ) {
      MARLEY_LOG( DEBUG, "physics.deexcitation.gamma" )
        << "  this level does not have any gammas";
      cascade_finished = true;
    }
    else {
 
      // Get the total decay width (MeV) for the initial level
      level_total_width = marley_utils::hbar * marley_utils::log_2
        / p_current_level->half_life();
 
      // Update the current level now that the gamma has been emitted
      p_current_level = p_gamma->end_level();
      if ( !p_current_level ) {
        throw marley::Error( "This gamma does not have an end level."
          " Cannot continue cascade." );
      }
      MARLEY_LOG( DEBUG, "physics.deexcitation.gamma" )
        << std::setprecision(15) << std::scientific
        << "  emitted gamma with energy "
        << p_gamma->energy() << " MeV. New level has energy "
        << p_current_level->energy() << " MeV.";
 
      // Get the excitation energy of the end level. This will be added to
      // the ground state mass of the nucleus to determine its
      // post-gamma-emission mass.
      double Exf = p_current_level->energy();
 
      // Create new particle objects to represent the emitted gamma and
      // recoiling nucleus
      auto gamma = marley_hepmc3::make_particle( marley_utils::PHOTON,
        marley_hepmc3::NUHEPMC_FINAL_STATE_STATUS, 0.0 );
 
      int pdg = marley_utils::get_nucleus_pid( Z_, A_ );
 
      auto nucleus = marley_hepmc3::make_particle( pdg,
        marley_hepmc3::NUHEPMC_INTERMEDIATE_RESIDUE_STATUS,
        mt.get_atomic_mass(pdg) + Exf
        - qIon*mt.get_particle_mass(marley_utils::ELECTRON) );
 
      // Create a new binary decay vertex
      auto decay_vtx = std::make_shared< HepMC3::GenVertex >();
      decay_vtx->set_status( marley_hepmc3::NUHEPMC_GAMMA_DECAY_VERTEX );
 
      decay_vtx->add_particle_in( residue );
      decay_vtx->add_particle_out( gamma );
      decay_vtx->add_particle_out( nucleus );
 
      // Sample a decay time (MeV^{-1}) for emission of the chosen gamma-ray
      // and store this timing information in the new binary decay vertex
      double gamma_partial_width = gamma_branching_ratio * level_total_width;
      marley_hepmc3::store_decay_time( gamma_partial_width, gen, decay_vtx,
        residue );
 
      // Add the decay vertex to the event record
      event.add_vertex( decay_vtx );
 
      // We can set the charge attribute now that the daughter nucleus
      // belongs to the decay vertex (and thus the parent event)
      marley_hepmc3::set_particle_charge( *nucleus, qIon );
 
      // We can also now set the attributes representing the
      // excitation energy, spin, and parity of the daughter nucleus
      nucleus->add_attribute( "Ex",
        std::make_shared< HepMC3::DoubleAttribute >(Exf) );
      nucleus->add_attribute( "twoJ",
        std::make_shared< HepMC3::IntAttribute >( p_current_level->twoJ() )
      );
      nucleus->add_attribute( "parity", std::make_shared< HepMC3::IntAttribute >(
        static_cast< int >(p_current_level->parity()) )
      );
 
      // If the total decay width is finite (equivalently, there is a tabulated
      // value for the level half-life), then store it in an attribute
      // attached to the vertex
      if ( std::isfinite(level_total_width) ) {
        decay_vtx->add_attribute( "TotalWidth",
          std::make_shared<HepMC3::DoubleAttribute>(level_total_width)
        );
      }
 
      // Also store the selected gamma-ray's branching ratio
      decay_vtx->add_attribute( "GammaBranchingRatio",
        std::make_shared< HepMC3::DoubleAttribute>( gamma_branching_ratio )
      );
 
      // Sample a direction assuming that the gammas are emitted isotropically
      // in the nucleus's rest frame.
      // sample from [-1, 1]
      double gamma_cos_theta = gen.uniform_random_double( -1.0, 1.0, true );
      // sample from [0, 2*pi)
      double gamma_phi = gen.uniform_random_double( 0., 2.*marley_utils::pi,
        false );
 
      MARLEY_LOG( TRACE, "physics.deexcitation.gamma" )
        << "  sampled gamma direction: cos_theta = " << gamma_cos_theta
        << ", phi = " << gamma_phi << " rad";
 
      // Determine the final energies and momenta for the recoiling nucleus and
      // emitted gamma ray. Store them in the final state particle objects.
      marley_kinematics::two_body_decay( residue, gamma, nucleus,
        gamma_cos_theta, gamma_phi );
 
      // Update the residue for this event to take into account changes from
      // gamma ray emission
      residue.swap( nucleus );
    }
  }
 
  MARLEY_LOG( DEBUG, "physics.deexcitation.gamma" )
    << "Finished gamma cascade at level with energy "
    << p_current_level->energy();
 
  residue->set_status( marley_hepmc3::NUHEPMC_FINAL_STATE_STATUS );
}

References A_, marley::Gamma::end_level(), marley::Gamma::energy(), marley::Level::energy(), marley::MassTable::get_atomic_mass(), marley::MassTable::get_particle_mass(), marley::Level::half_life(), marley::MassTable::Instance(), marley::Level::parity(), pdg(), marley::Level::sample_gamma(), HepMC3::GenParticle::set_status(), marley::Level::twoJ(), marley::Generator::uniform_random_double(), and Z_.

Referenced by marley::NucleusDecayer::process_event().

get_levels()

const std::vector< std::unique_ptr< marley::Level > > & marley::DecayScheme::get_levels ( ) const

inline

Get a const reference to the vector that holds the Level objects.

Definition at line 158 of file DecayScheme.hh.

{ return levels_; }

References levels_.

get_pointer_to_closest_level()

marley::Level * marley::DecayScheme::get_pointer_to_closest_level

(

double

E_level

)

Gets a pointer to the Level in the DecayScheme whose excitation energy is closest to E_level.

Parameters

E_level

excitation energy (MeV)

Note

Returns nullptr if the DecayScheme doesn't own any Level objects

Definition at line 44 of file DecayScheme.cc.

{
  // If this decay scheme doesn't own any levels, return nullptr immediately
  size_t num_levels = levels_.size();
  if (num_levels == 0) return nullptr;
 
  // Search for the level whose energy is closest to the given value of E_level
  size_t e_index = level_lower_bound_index(E_level);
 
  if (e_index == num_levels) {
    // The given energy is greater than every level energy in our decay scheme.
    // We will therefore assume that the desired level is the highest level.
    // Its index is given by one less than the number of elements in the sorted
    // vector, so subtract one from our previous result.
    --e_index;
  }
  else if (e_index > 0) {
    // If the calculated index does not correspond to the first element, we
    // still need to check which of the two levels found (one on each side) is
    // really the closest. Do so and reassign the index if needed.
    if (std::abs(E_level - levels_.at(e_index)->energy())
      > std::abs(E_level - levels_.at(e_index - 1)->energy()))
    {
      --e_index;
    }
  }
 
  // Return a pointer to the selected level object
  return levels_.at(e_index).get();
}

References level_lower_bound_index(), and levels_.

Referenced by marley::NucleusDecayer::process_event().

level_lower_bound_index()

size_t marley::DecayScheme::level_lower_bound_index ( double Ex )

protected

Get the index of the first level whose energy is not less than Ex.

Parameters

Ex	Excitation energy (MeV)

Note

Returns levels_.size() if all levels have energies below Ex

Definition at line 435 of file DecayScheme.cc.

                                                             {
  const auto E_begin = marley::Level::make_energy_iterator( levels_.cbegin() );
  const auto E_end = marley::Level::make_energy_iterator( levels_.cend() );
 
  const auto closest_E_iter = std::lower_bound( E_begin, E_end, Ex );
  return std::distance( E_begin, closest_E_iter );
}

References levels_, and marley::Level::make_energy_iterator().

Referenced by add_level(), and get_pointer_to_closest_level().

pdg()

int marley::DecayScheme::pdg

(

)

const

Returns the nuclear PDG code corresponding to Z and A.

Definition at line 76 of file DecayScheme.cc.

                                 {
  return marley_utils::get_nucleus_pid( Z_, A_ );
}

References A_, and Z_.

Referenced by do_cascade().

print()

void marley::DecayScheme::print

(

std::ostream &

out = std::cout

)

const

Print this DecayScheme object to a std::ostream.

Definition at line 457 of file DecayScheme.cc.

                                                     {
 
  size_t num_levels = levels_.size();
 
  out << Z_ << ' ' << A_ << ' ' << num_levels << '\n';
 
  for ( const auto& lev : levels_ ) {
    out << "  " << lev->energy() << ' ' << lev->twoJ() << ' '
      << lev->parity() << ' ' << lev->gammas().size() << '\n';
    for ( const auto& g : lev->gammas() ) {
      out << "    " << g.energy() << ' ' << g.relative_intensity();
 
      const auto cit = std::find_if( levels_.cbegin(), levels_.cend(),
        [&g]( const std::unique_ptr< marley::Level >& l )
        -> bool { return l.get() == g.end_level(); } );
 
      int level_f_idx = -1;
      if ( cit != levels_.cend() ) {
        level_f_idx = std::distance(levels_.cbegin(), cit);
      }
      out << " " << level_f_idx << '\n';
    }
  }
}

References A_, levels_, and Z_.

print_latex_table()

void marley::DecayScheme::print_latex_table

(

std::ostream &

ostr = std::cout

)

Print LaTeX source code that gives a tabular representation of the DecayScheme object.

Definition at line 368 of file DecayScheme.cc.

                                                            {
 
  std::string nuc_id = marley_utils::nuc_id( Z_, A_ );
 
  std::string caption_beginning =
    std::string("{\\textbf{Levels") +
    " and $\\boldsymbol{\\gamma}$ transitions \n for " +
    "\\isotope[\\boldsymbol{" +
    marley_utils::trim_copy(nuc_id.substr(0,3)) +
    "}]{\\textbf{" + nuc_id.substr(3,1) +
    marley_utils::trim_copy(
      marley_utils::to_lowercase(nuc_id.substr(4,1))) +
    "}} \n";// from file " + filename + " ";
 
  ostr << marley_utils::latex_table_1;
 
  ostr << caption_beginning + "}}\\\\\n";
 
  ostr << marley_utils::latex_table_2;
 
  ostr << caption_beginning + " -- \\textit{continued}}} \\\\\n";
 
  ostr << marley_utils::latex_table_3;
 
  // Cycle through each of the levels owned by this decay scheme
  // object in order of increasing energy
  for (const auto& lev : levels_ ) {
 
    std::string sp = lev->spin_parity_string();
 
    ostr << lev->energy() << " & " << sp  << " & ";
 
    const auto& gammas = lev->gammas();
 
    // If there aren't any gammas for this level, finish writing
    // the current row of the table. Add extra space between this
    // level and the next one.
    if ( gammas.empty() ) {
      ostr << " &  &";
      // If this is the last row of the table, don't add extra space.
      if (lev == levels_.back()) ostr << '\n';
      else ostr << " \\\\ \\addlinespace[\\ExtraRowSpace]\n";
    }
 
    // Cycle through each of the gammas owned by the current level
    for ( const auto& g : gammas ) {
      // If this is not the first gamma, add empty columns
      // for the level energy and spin-parity
      if ( &g != &gammas.front() ) ostr << " & & ";
      // Output information about the current gamma
      ostr << g.energy() << " & " << g.relative_intensity()
        << " & " << g.end_level()->energy();
      // Add vertical space after the final gamma row. Also prevent page breaks
      // in the middle of a list of gammas by outputting a star at the end of
      // each row except the final gamma row.
      if ( &g == &gammas.back() ) {
    // Don't add the extra row space for the very last row in the table
    if ( lev == levels_.back() ) ostr << '\n';
        else ostr << " \\\\ \\addlinespace[\\ExtraRowSpace]" << '\n';
      }
      else ostr << " \\\\*" << '\n';
    }
  }
  ostr << marley_utils::latex_table_4 << '\n';
}

References A_, levels_, and Z_.

print_report()

void marley::DecayScheme::print_report

(

std::ostream &

ostr = std::cout

)

const

Print a human-readable text representation of the DecayScheme object.

Definition at line 341 of file DecayScheme.cc.

                                                             {
  // Cycle through each of the levels owned by this decay scheme
  // object in order of increasing energy
  for ( const auto& lev : levels_ ) {
    int twoj = lev->twoJ();
    std::string spin = std::to_string( twoj / 2 );
    // If 2*J is odd, then the level has half-integer spin
    if ( twoj % 2 ) spin += "/2";
    marley::Parity parity = lev->parity();
 
    ostr << "Level at " << lev->energy() << " MeV has spin-parity "
      << spin << parity << " and half-life " << lev->half_life() << " s\n";
 
    std::vector< marley::Gamma >& gammas = lev->gammas();
 
    // Cycle through each of the gammas owned by the current level
    // (according to the ENSDF specification, these will already be
    // sorted in order of increasing energy)
    for ( const auto& g : gammas ) {
      ostr << "  has a gamma with energy " << g.energy() << " MeV";
      ostr << " (transition to level at "
        << g.end_level()->energy() << " MeV)" << '\n';
      ostr << "    and relative intensity " << g.relative_intensity() << '\n';
    }
  }
}

References levels_.

read_from_stream()

void marley::DecayScheme::read_from_stream

(

std::istream &

in

)

Use a std::istream to initialize this DecayScheme object, replacing any previous data.

The expected format for data in the std::istream is the same as the output format in DecayScheme::print()

Definition at line 482 of file DecayScheme.cc.

                                                         {
 
  levels_.clear();
 
  int num_levels;
  in >> Z_ >> A_ >> num_levels;
 
  // If we had trouble parsing the decay scheme header, then
  // just return the stream without doing anything else.
  if ( !in ) return;
 
  double energy, ri, half_life;
  int two_j, num_gammas, level_f_idx;
  marley::Parity pi;
 
  for ( int i = 0; i < num_levels; ++i ) {
    in >> energy >> two_j >> pi >> num_gammas >> half_life;
 
    marley::Level& l = add_level( marley::Level(energy, two_j, pi, half_life) );
    for ( int j = 0; j < num_gammas; ++j ) {
      in >> energy >> ri >> level_f_idx;
      l.add_gamma( energy, ri, levels_.at(level_f_idx).get() );
    }
  }
 
  // Remove levels above the unbound threshold (these will always be
  // handled using a continuous level density treatment)
  const auto& mt = marley::MassTable::Instance();
  double unbound_Ex = mt.unbound_threshold( Z_, A_ );
 
  // TODO: replace with std::erase_if when updating to C++20
  auto iter_new_end = std::remove_if( levels_.begin(), levels_.end(),
    [ unbound_Ex ]( const std::unique_ptr< marley::Level >& lev ) -> bool {
      double lvl_Ex = lev->energy();
      bool unbound = ( lvl_Ex > unbound_Ex );
      return unbound;
    }
  );
  levels_.erase( iter_new_end, levels_.end() );
 
}

References A_, marley::Level::add_gamma(), add_level(), marley::MassTable::Instance(), levels_, and Z_.

Z()

int marley::DecayScheme::Z ( ) const

inline

Get the atomic number.

Definition at line 154 of file DecayScheme.hh.

{ return Z_; }

References Z_.

Referenced by DecayScheme(), and DecayScheme().

Member Data Documentation

A_

int marley::DecayScheme::A_

protected

Mass number.

Definition at line 122 of file DecayScheme.hh.

Referenced by DecayScheme(), DecayScheme(), A(), do_cascade(), pdg(), print(), print_latex_table(), and read_from_stream().

levels_

std::vector< std::unique_ptr<marley::Level> > marley::DecayScheme::levels_

protected

Level objects owned by this DecayScheme.

Definition at line 125 of file DecayScheme.hh.

Referenced by add_level(), get_levels(), get_pointer_to_closest_level(), level_lower_bound_index(), print(), print_latex_table(), print_report(), and read_from_stream().

Z_

int marley::DecayScheme::Z_

protected

Atomic number.

Definition at line 121 of file DecayScheme.hh.

Referenced by DecayScheme(), DecayScheme(), do_cascade(), pdg(), print(), print_latex_table(), read_from_stream(), and Z().

---

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

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