Computes event weights based on variations to the optical model potential. More...

#include <OMPWeightCalculator.hh>

Inheritance diagram for marley::OMPWeightCalculator:

Public Member Functions
	OMPWeightCalculator (const marley::JSON &config)

virtual double	weight (HepMC3::GenEvent &event, marley::Generator &gen) const override
	Compute the weight for the given event.

Public Member Functions inherited from marley::WeightCalculator
	WeightCalculator (const marley::JSON &config)

	WeightCalculator (const std::string &name)

const std::string &	name () const

Protected Attributes
std::shared_ptr< marley::StructureDatabase >	sdb_
	Owned StructureDatabase with a custom set of optical model parameters.

Protected Attributes inherited from marley::WeightCalculator
std::string	name_
	Name used to label the output event weight.

Detailed Description

Computes event weights based on variations to the optical model potential.

Definition at line 24 of file OMPWeightCalculator.hh.

Constructor & Destructor Documentation

◆ OMPWeightCalculator()

marley::OMPWeightCalculator::OMPWeightCalculator ( const marley::JSON & config )

Definition at line 40 of file OMPWeightCalculator.cc.

  : marley::WeightCalculator( config )
{
  if ( !config.is_object() ) {
    throw marley::Error( "Non-object JSON configuration passed to constructor"
      " of marley::OMPWeightCalculator" );
  }
 
  if ( !config.has_key("opt_mod") ) {
    throw marley::Error( "Missing \"opt_mod\" key in marley::OMPWeight"
      "Calculator JSON configuration" );
  }
 
  // Load the JSON configuration for the optical model parameters
  const auto& om_config = config.at( "opt_mod" );
  if ( !om_config.is_object() ) {
    throw marley::Error( "Non-object associated with \"opt_mod\" key in"
      " marley::OMPWeightCalculator JSON configuration" );
  }
 
  // Create a structure database with a custom set of optical model parameters
  sdb_ = std::make_shared< marley::StructureDatabase >();
  sdb_->load_optical_model_params( &om_config );
}

Member Function Documentation

◆ weight()

double marley::OMPWeightCalculator::weight	(	HepMC3::GenEvent &	event,
		marley::Generator &	gen ) const

overridevirtual

Compute the weight for the given event.

Note: If a derived class uses random numbers, it MUST own its own random number generator (seeded with a fixed value from its JSON config) rather than using the Generator's RNG. Random numbers MUST also only be used by derived classes during initialization, not on an event-by-event basis. These rules ensure that the calculated weights remain deterministic and replayable regardless of when or in what context the calculator is used. In particular, this allows the "resume" behavior of the "marley generate" command to be correct even after reweighting has been run on an existing sample.

Todo: Revisit this constraint if a clear use case is found for event-by-event random numbers sampled within a WeightCalculator derived class implementation.

Implements marley::WeightCalculator.

Definition at line 65 of file OMPWeightCalculator.cc.

{
  // Start with a result of unity
  double weight = 1.;
 
  // Get a vector of pointers to all decay vertices in the event that
  // were handled using the Hauser-Feshbach model
  auto hf_vtx_vec = marley_hepmc3::get_vertices_with_status(
    marley_hepmc3::NUHEPMC_HF_DECAY_VERTEX, ev );
 
  // Include each vertex that was found in the reweighting calculation
  for ( const auto& vtx : hf_vtx_vec ) {
 
    // Start with the default assumption that the decay was to a discrete
    // nuclear level
    bool decayed_to_continuum = false;
 
    const auto& in_vec = vtx->particles_in();
    const auto& out_vec = vtx->particles_out();
 
    // Double-check that the expected numbers of particles are present
    if ( in_vec.size() != 1u || out_vec.size() != 2u ) {
    throw marley::Error( "Hauser-Feshbach vertex encountered that"
      " is not a binary decay" );
    }
 
    // Retrieve the decay widths computed in the original simulation
    auto width_tot_attr = vtx->attribute< HepMC3::DoubleAttribute >(
      "TotalWidth" );
    double width_tot = width_tot_attr->value();
 
    auto width_ec_attr = vtx->attribute< HepMC3::DoubleAttribute >(
      "ECWidth" );
    double width_ec = width_ec_attr->value();
 
    // This attribute is only present for decays to the continuum
    auto width_sp_attr = vtx->attribute< HepMC3::DoubleAttribute >(
      "SPWidth" );
    double width_sp = 0.;
    if ( width_sp_attr ) {
      decayed_to_continuum = true;
      width_sp = width_sp_attr->value();
    }
 
    // Get access to the mother nucleus for the current binary decay
    const auto& mother = in_vec.front();
 
    // Retrieve the starting nuclear excitation energy, spin, and parity
    // TODO: add error handling here for missing attributes
    auto Exi_attr = mother->attribute< HepMC3::DoubleAttribute >( "Ex" );
    double Exi = Exi_attr->value();
 
    auto twoJi_attr = mother->attribute< HepMC3::IntAttribute >( "twoJ" );
    int twoJi = twoJi_attr->value();
 
    auto Pi_attr = mother->attribute< HepMC3::IntAttribute >( "parity" );
    marley::Parity Pi( Pi_attr->value() );
 
    // NucleusDecayer populates the outgoing particles in the decay
    // vertex in a specific order. The first is the emitted nuclear
    // fragment or gamma-ray, while the second is the daughter nucleus.
    const auto& emitted_particle = out_vec.front();
    const auto& daughter = out_vec.back();
 
    int emitted_pdg = emitted_particle->pid();
 
    // Retrieve the final nuclear excitation energy, spin, and parity
    // TODO: add error handling here for missing attributes
    auto Exf_attr = daughter->attribute< HepMC3::DoubleAttribute >( "Ex" );
    double Exf = Exf_attr->value();
 
    auto twoJf_attr = daughter->attribute< HepMC3::IntAttribute >( "twoJ" );
    int twoJf = twoJf_attr->value();
 
    auto Pf_attr = daughter->attribute< HepMC3::IntAttribute >( "parity" );
    marley::Parity Pf( Pf_attr->value() );
 
    // Construct a set of decay widths calculated with alternative settings
    marley::HauserFeshbachDecay hf_alt( mother, Exi, twoJi, Pi, *sdb_ );
 
    // Get the total decay width under the alternative calculation
    double width_tot_alt = hf_alt.total_width();
 
    // Find the ExitChannel in the alternative calculation corresponding
    // to the original decay that was sampled
    const auto& ec_vec = hf_alt.exit_channels();
    auto ec_iter = std::find_if( ec_vec.cbegin(), ec_vec.cend(),
    [ emitted_pdg, decayed_to_continuum, Exf ](
      const std::unique_ptr< marley::ExitChannel >& test_ec ) -> bool
    {
      if ( emitted_pdg != test_ec->emitted_particle_pdg() ) return false;
 
      if ( decayed_to_continuum ) {
      if ( !test_ec->is_continuum() ) return false;
      else return true;
      }
 
      // If we get to here, then we're dealing with a discrete transition
      if ( test_ec->is_continuum() ) return false;
 
      const auto* dec = dynamic_cast<
      const marley::DiscreteExitChannel* >( test_ec.get() );
      if ( !dec ) return false;
 
      // Check the excitation energy of the final discrete level as
      // a last confirmation that we've found the correct ExitChannel
      const auto& lev = dec->get_final_level();
      double Ex_level = lev.energy();
 
      // To deal with possible numerical precision issues, this check of
      // the final excitation energy allows for small differences
      if ( std::abs(Exf - Ex_level) > PRETTY_SMALL ) return false;
      return true;
    }
    );
 
    if ( ec_iter == ec_vec.cend() ) {
      MARLEY_LOG( WARN, "physics.opticalmodel" )
        << "Could not find ExitChannel during reweighting";
      weight = 0.;
      continue;
    }
 
    // Partial width for the chosen ExitChannel under the alternative
    // calculation
    double width_ec_alt = ( *ec_iter )->width();
 
    // Partial differential width for the chosen spin-parity under the
    // alternative calculation (applies only to decays to the continuum)
    double width_sp_alt = 0.;
    if ( decayed_to_continuum ) {
      // Evaluate the differential width at the sampled excitation energy.
      // This will populate the table of SpinParityWidth objects with the
      // correct values.
      const auto& cec = dynamic_cast<
        const marley::ContinuumExitChannel& >( *(*ec_iter) );
      cec.differential_width( Exf, true );
 
      // Retrieve extra information based on the kind of emitted particle
      int mpol = 0, two_j_frag = 0, orb_l = 0;
      bool emitted_gamma = false;
      auto mpol_attr = vtx->attribute< HepMC3::IntAttribute >(
        "multipolarity" );
      if ( mpol_attr ) {
        emitted_gamma = true;
        mpol = mpol_attr->value();
      }
      else {
        auto two_j_frag_attr = vtx->attribute< HepMC3::IntAttribute >(
          "two_j_frag" );
        auto orb_l_attr = vtx->attribute< HepMC3::IntAttribute >( "orb_l" );
 
        // TODO: add error handling for missing attributes here
 
        two_j_frag = two_j_frag_attr->value();
        orb_l = orb_l_attr->value();
      }
 
      // Find the SpinParityWidth object corresponding to the spin-parity
      // value that was actually sampled
      const auto& spw_vec = cec.get_spw_table();
      auto spw_iter = std::find_if( spw_vec.cbegin(), spw_vec.cend(),
        [ twoJf, Pf, emitted_gamma, mpol, two_j_frag, orb_l ](
          const std::unique_ptr< marley::ContinuumExitChannel
          ::SpinParityWidth >& spw ) -> bool
        {
          if ( Pf != spw->Pf ) return false;
          if ( twoJf != spw->twoJf ) return false;
          if ( emitted_gamma ) {
            const auto* g_spw = static_cast< const marley
              ::GammaContinuumExitChannel::GammaSpinParityWidth* >(
              spw.get() );
            if ( !g_spw ) return false;
            if ( mpol != g_spw->multipolarity ) return false;
          }
          else {
            // emitted fragment
            const auto* f_spw = static_cast< const marley
              ::FragmentContinuumExitChannel::FragmentSpinParityWidth* >(
              spw.get() );
            if ( !f_spw ) return false;
            if ( two_j_frag != f_spw->two_j_frag ) return false;
            if ( orb_l != f_spw->orb_l ) return false;
          }
          return true;
        }
      );
 
      if ( spw_iter == spw_vec.cend() ) {
        MARLEY_LOG( WARN, "physics.opticalmodel" )
          << "Could not find SpinParityWidth during reweighting";
        weight = 0.;
        continue;
      }
 
      // Store the partial differential width to this spin-parity state
      // under the alternative calculation
      width_sp_alt = spw_iter->get()->diff_width;
    }
 
    // We've accumulated all the width values we need to compute an
    // event weight. Do some sanity checks beforehand.
    std::string bad_width_name;
    if ( width_tot_alt <= 0. ) {
      bad_width_name = "Alternate total_width";
    }
    else if ( width_tot <= 0. ) {
      bad_width_name = "Original total width";
    }
    else if ( width_ec <= 0. ) {
      bad_width_name = "Original exit channel";
    }
    else if ( width_ec_alt <= 0. ) {
      bad_width_name = "Alternate exit channel";
    }
    else if ( decayed_to_continuum ) {
      if ( width_sp_alt <= 0. ) bad_width_name = "Alternate spin-parity";
      else if ( width_sp <= 0. ) bad_width_name = "Original spin-parity";
    }
 
    if ( !bad_width_name.empty() ) {
      MARLEY_LOG( WARN, "physics.opticalmodel" ) << bad_width_name << " is non-positive";
      weight = 0.;
      continue;
    }
 
    // Everything looks okay, so do the weight calculation using the
    // appropriate expression for a decay to the continuum or a discrete
    // nuclear level
    double  w = width_tot / width_tot_alt;
    if ( decayed_to_continuum ) {
      w *= width_sp_alt / width_sp;
    }
    else {
      w *= width_ec_alt / width_ec;
    }
 
    MARLEY_LOG( DEBUG, "physics.opticalmodel" ) << "OMP reweight: vertex"
      " weight w = " << w << " (tot_width_alt/tot_width = "
      << width_tot_alt << "/" << width_tot << ")";
 
    // Multiply the weight for the current decay vertex into the overall
    // event weight
    weight *= w;
 
  } // decay vertex loop
 
  MARLEY_LOG( DEBUG, "physics.opticalmodel" ) << "OMP reweight: final event"
    " weight = " << weight;
  return weight;
}

References marley::Level::energy(), marley::HauserFeshbachDecay::exit_channels(), marley::DiscreteExitChannel::get_final_level(), marley::ExitChannel::is_continuum(), sdb_, HepMC3::DoubleAttribute::value(), HepMC3::IntAttribute::value(), and weight().

Referenced by weight().

Member Data Documentation

◆ sdb_

std::shared_ptr< marley::StructureDatabase > marley::OMPWeightCalculator::sdb_

protected

Owned StructureDatabase with a custom set of optical model parameters.

Definition at line 38 of file OMPWeightCalculator.hh.