marley::NucleusDecayer Class Reference

EventProcessor that handles nuclear de-excitations. More...

#include <NucleusDecayer.hh>

Inheritance diagram for marley::NucleusDecayer:

Public Member Functions
virtual void	process_event (HepMC3::GenEvent &event, marley::Generator &gen) override
	Processes an input GenEvent object.

Detailed Description

EventProcessor that handles nuclear de-excitations.

Definition at line 23 of file NucleusDecayer.hh.

Constructor & Destructor Documentation

NucleusDecayer()

marley::NucleusDecayer::NucleusDecayer ( )

inline

Definition at line 27 of file NucleusDecayer.hh.

{}

Member Function Documentation

process_event()

void marley::NucleusDecayer::process_event ( HepMC3::GenEvent & ev, marley::Generator & gen )

overridevirtual

Processes an input GenEvent object.

Parameters

[in,out]	ev	The GenEvent object to be processed
[in]	gen	If needed, the Generator object to use during processing (e.g., for obtaining random numbers)

Todo

Revisit this

Implements marley::EventProcessor.

Definition at line 53 of file NucleusDecayer.cc.

{
  auto undecayed_residues = marley_hepmc3::get_particles_with_status(
    marley_hepmc3::NUHEPMC_UNDECAYED_RESIDUE_STATUS, event );
 
  MARLEY_LOG( DEBUG, "physics.deexcitation" ) << "NucleusDecayer: processing "
    << undecayed_residues.size() << " undecayed residue(s)";
 
  // Check the reaction process that created this event. The process types
  // distinguish between discrete and continuum reactions, which is helpful
  // below.
  int proc_id = event.attribute< HepMC3::IntAttribute >(
    "signal_process_id" )->value();
  auto proc_type = marley_hepmc3::from_nuhepmc_proc_id( proc_id );
  bool is_continuum_channel = false;
  if ( proc_type == marley::Reaction::ProcessType::NeutrinoCC_Continuum
    || proc_type == marley::Reaction::ProcessType::AntiNeutrinoCC_Continuum 
    || proc_type == marley::Reaction::ProcessType::NC_Continuum )
  {
    is_continuum_channel = true;
  }
 
  for ( auto residue : undecayed_residues ) {
 
    // Get the residue excitation energy from the event. These values represent
    // its state immediately following the initial two-two scattering reaction.
    double Ex = residue->attribute< HepMC3::DoubleAttribute >( "Ex" )->value();
    int twoJ = residue->attribute< HepMC3::IntAttribute >( "twoJ" )->value();
    int p_int = residue->attribute< HepMC3::IntAttribute >( "parity" )->value();
    marley::Parity P( p_int );
 
    // If the residue is in its ground state, then there's nothing for us to do.
    // Just continue the loop without comment.
    if ( Ex == 0. ) continue;
 
    MARLEY_LOG( DEBUG, "physics.deexcitation" ) << "De-exciting residue PDG "
      << residue->pid() << ": Ex = " << Ex << " MeV, 2J = " << twoJ
      << ", P = " << P;
 
    // The excitation energy should be nonnegative. Complain if it's not.
    if ( Ex < 0. ) throw marley::Error("Negative excitation energy Ex = "
      + std::to_string(Ex) + " MeV encountered in marley::NucleusDecayer::"
      "deexcite_residue()");
 
    // To prevent accidental double application of the de-excitation cascade,
    // check that the residue mass is consistent with the excitation energy
    // stored in the event record (and thus was never decayed).
    const auto& mt = marley::MassTable::Instance();
    int initial_residue_pdg = residue->pid();
    int qIon = marley_hepmc3::get_particle_charge( *residue );
 
    // Check that the residue PDG code makes sense. If it's not a nucleus,
    // warn the user and refuse to do the cascade.
    if ( !marley_utils::is_ion(initial_residue_pdg) ) {
      MARLEY_LOG( WARN, "physics.deexcitation" )
        << "Unrecognized nuclear PDG code "
        << initial_residue_pdg << " encountered in marley::NucleusDecayer::"
        << "deexcite_residue(). The de-excitation cascade will be skipped";
      continue;
    }
 
    double residue_mass = residue->generated_mass();
 
    // Ground-state residue mass
    double gs_residue_mass = mt.get_atomic_mass( initial_residue_pdg )
      - qIon*mt.get_particle_mass( marley_utils::ELECTRON );
 
    double expected_residue_mass = gs_residue_mass + Ex;
 
    if ( std::abs(residue_mass - expected_residue_mass) > EX_TOLERANCE ) {
 
      if ( std::abs(residue_mass - gs_residue_mass) <= EX_TOLERANCE ) {
        MARLEY_LOG( WARN, "physics.deexcitation" )
          << "Encountered ground-state nuclear remnant"
          << " in marley::NucleusDecay::deexcite_residue(). The de-excitation"
          << " cascade has already been applied.";
        continue;
      }
 
      // If we get here, then the residue is not in its ground state but also
      // not in the initial excited state given in the event record. Something
      // went wrong with a partial application of a de-excitation cascade.
      // Throw an error rather than trying to figure out how to do the right
      // thing.
      throw marley::Error("Partially de-excited nuclear remnant encountered"
        " in marley::NucleusDecay::deexcite_residue().");
    }
 
    // Decide whether we need to start the de-excitation cascade from a
    // discrete nuclear level or from the continuum. Do this by comparing the
    // excitation energy from the event record to the "unbound threshold" for
    // the residue. If we're above the unbound threshold, do a continuum decay.
    // Also start with a continuum decay if no discrete level data are
    // available for the residue.
    auto* ds = gen.get_structure_db().get_decay_scheme( initial_residue_pdg );
    double unbound_threshold = mt.unbound_threshold( initial_residue_pdg );
 
    // If Reaction::set_level_ptrs() changes, you'll want to change this too.
    // TODO: find a better way of keeping the two pieces of code in sync
    // TODO: numerical round-off can cause issues with the first test, so you
    // should revisit it again when interfacing MARLEY with other codes
    // that use it solely as a de-excitation model. For now, the third
    // option in the logical OR prevents issues with numerical round-off near
    // the unbound threshold.
    bool continuum = ( Ex > unbound_threshold )
      || ( !ds ) || ( is_continuum_channel );
 
    // Keep track of whether the cascade was started from the continuum
    // or not. If it was started from a discrete level, we'll double-check that
    // discrete level's excitation energy below.
    bool started_from_continuum = continuum;
 
    MARLEY_LOG( DEBUG, "physics.deexcitation" ) << "De-excitation path: "
      << ( continuum ? "continuum (Hauser-Feshbach)" : "discrete gamma cascade" );
 
    if ( continuum ) {
 
      // Particles used for storage of binary decay products during the
      // de-excitation cascade
      auto first = std::make_shared< HepMC3::GenParticle >();
      auto second = std::make_shared< HepMC3::GenParticle >();
 
      // The selected level is unbound, so handle its de-excitation using
      // the Hauser-Feshbach statistical model.
      while ( continuum && Ex > CONTINUUM_GS_CUTOFF ) {
 
        auto& sdb = gen.get_structure_db();
 
        marley::HauserFeshbachDecay hfd( residue, Ex, twoJ, P, sdb );
        MARLEY_LOG( DEBUG, "physics.deexcitation.hauser" ) << hfd;
 
        int q_second;
        const auto& exit_channel = hfd.do_decay( Ex, twoJ, P, first, second,
          q_second, gen );
 
        continuum = exit_channel.is_continuum();
 
        double width_tot = hfd.total_width();
        double width_ec = exit_channel.width();
 
        MARLEY_LOG( DEBUG, "physics.deexcitation.hauser" )
          << "Hauser-Feshbach decay to " << first->pid()
          << " and " << second->pid();
        MARLEY_LOG( DEBUG, "physics.deexcitation.hauser" )
          << second->pid() << " is at Ex = " << Ex << " MeV.";
 
        // Create a new binary decay vertex
        auto decay_vtx = std::make_shared< HepMC3::GenVertex >();
        decay_vtx->set_status( marley_hepmc3::NUHEPMC_HF_DECAY_VERTEX );
 
        decay_vtx->add_particle_in( residue );
        decay_vtx->add_particle_out( first );
        decay_vtx->add_particle_out( second );
 
        // Sample a decay time (MeV^{-1}) for emission of the chosen particle
        // and store this timing information in the new binary decay vertex
        marley_hepmc3::store_decay_time( width_ec, gen, decay_vtx, residue );
 
        // Add the decay vertex to the event record
        event.add_vertex( decay_vtx );
 
        // We can now set the charge of the daughter ion because it belongs
        // to the parent event (through the decay vertex)
        marley_hepmc3::set_particle_charge( *second, q_second );
 
        // We can also now set the attributes representing the daughter ion's
        // excitation energy, spin, and parity
        second->add_attribute( "Ex",
          std::make_shared< HepMC3::DoubleAttribute >(Ex) );
        second->add_attribute( "twoJ",
          std::make_shared< HepMC3::IntAttribute >(twoJ) );
        second->add_attribute( "parity",
          std::make_shared< HepMC3::IntAttribute >(static_cast<int>( P )) );
 
        // The daughter ion now takes the role of the residue for the next loop
        // iteration
        residue.swap( second );
 
        // Store some information about the total and partial widths of
        // the simulated compound nucleus decay in attributes attached to
        // the decay vertex
        decay_vtx->add_attribute( "TotalWidth",
          std::make_shared< HepMC3::DoubleAttribute >(width_tot) );
 
        decay_vtx->add_attribute( "ECWidth",
          std::make_shared< HepMC3::DoubleAttribute >(width_ec) );
 
        // In the case of a transition to the continuum, also store the partial
        // differential width for the chosen spin-parity of the daughter
        // nucleus
        if ( continuum ) {
          const auto& cec = dynamic_cast< const marley::ContinuumExitChannel& >(
            exit_channel );
 
          const auto* spw_ptr = cec.get_last_sampled_spw();
          double width_sp = spw_ptr->diff_width;
 
          decay_vtx->add_attribute( "SPWidth",
            std::make_shared< HepMC3::DoubleAttribute >(width_sp) );
 
          bool is_fragment_emission = exit_channel.emits_fragment();
 
          if ( is_fragment_emission ) {
            const auto* f_spw = static_cast< const marley
              ::FragmentContinuumExitChannel::FragmentSpinParityWidth* >(
              spw_ptr );
            decay_vtx->add_attribute( "two_j_frag",
              std::make_shared< HepMC3::IntAttribute >(f_spw->two_j_frag) );
            decay_vtx->add_attribute( "orb_l",
              std::make_shared< HepMC3::IntAttribute >(f_spw->orb_l) );
          }
 
          else {
            // Gamma-ray emission in the continuum
            const auto* g_spw = static_cast< const marley
              ::GammaContinuumExitChannel::GammaSpinParityWidth* >( spw_ptr );
 
            decay_vtx->add_attribute( "multipolarity",
              std::make_shared< HepMC3::IntAttribute >(g_spw->multipolarity) );
          }
        }
      }
    }
 
    if ( !continuum ) {
      // Either the selected initial level was bound (so it will only decay via
      // gamma emission) or the Hauser-Feshbach decay process has now accessed
      // a bound level in the residual nucleus. In either case, use gamma-ray
      // decay scheme data to sample the de-excitation gammas and add them to
      // this event's final particle list.
      marley::DecayScheme* dec_scheme = gen.get_structure_db()
        .get_decay_scheme( residue->pid() );
 
      // Start the gamma cascade from this discrete level
      marley::Level* lev = dec_scheme->get_pointer_to_closest_level( Ex );
 
      // If we get a null level pointer from the decay scheme, complain
      if ( !lev ) throw marley::Error( "Null nuclear level pointer encountered"
        " in marley::NucleusDecayer::deexcite_residue()" );
 
      // If we did not simulate any continuum decays before getting to this
      // point, then double-check that the excitation energy from the event
      // record and the initial level are consistent. If they're not, then
      // complain by throwing an error.
      if ( !started_from_continuum ) {
        double Ex_level = lev->energy();
        if ( std::abs(Ex - Ex_level) > EX_TOLERANCE ) {
          throw marley::Error( "Excitation energy mismatch encountered in"
            " marley::NucleusDecayer::deexcite_residue(). Event has Ex = "
            + std::to_string(Ex) + " MeV while the initial discrete level has "
            + std::to_string(Ex_level) + " MeV" );
        }
      }
 
      dec_scheme->do_cascade( *lev, event, gen, residue );
    }
 
  } // loop over undecayed residues
 
}

References marley::Reaction::AntiNeutrinoCC_Continuum, marley::ContinuumExitChannel::SpinParityWidth::diff_width, marley::DecayScheme::do_cascade(), marley::HauserFeshbachDecay::do_decay(), marley::Level::energy(), marley::StructureDatabase::get_decay_scheme(), marley::ContinuumExitChannel::get_last_sampled_spw(), marley::DecayScheme::get_pointer_to_closest_level(), marley::Generator::get_structure_db(), marley::MassTable::Instance(), marley::ExitChannel::is_continuum(), marley::Reaction::NC_Continuum, and marley::Reaction::NeutrinoCC_Continuum.

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

---

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

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