marley::StrengthVariationWeightCalculator Class Reference

WeightCalculator that varies nuclear matrix element strengths according to their experimental uncertainties. Three modes are available: "multisim" (dimidiated Gaussian random draws), "shift" (deterministic systematic shifts of all matrix elements), and "unisim" (deterministic systematic shifts of a single matrix element at a time). More...

#include <StrengthVariationWeightCalculator.hh>

Inheritance diagram for marley::StrengthVariationWeightCalculator:

Public Types
enum class	VariationMode { multisim , shift , unisim }
	Supported variation modes. More...

Public Member Functions
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

Static Public Member Functions
static std::vector< std::shared_ptr< StrengthVariationWeightCalculator > >	create_instances (const marley::JSON &config, marley::Generator &gen)
	Static factory: validates JSON configuration and creates all variation instances. Called by the Weighter for each strength_variation entry in the weights config array.

Protected Member Functions
void	ensure_initialized (marley::Generator &gen) const
	Lazy initialization of per-instance varied strengths.

Protected Attributes
const DiscreteNuclearReaction *	dnr_ = nullptr
	Pointer to the matched DiscreteNuclearReaction.
bool	initialized_ = false
	Whether lazy initialization has been completed.
size_t	me_idx_ = 0
	Index of the single matrix element to vary (unisim mode)
VariationMode	mode_
	Variation mode for this instance.
Reaction::ProcessType	process_type_ = Reaction::ProcessType::Unknown
	Process type of the matched Reaction.
std::string	resolved_reaction_file_
	Resolved path of the reaction input file of interest.
std::shared_ptr< std::mt19937_64 >	rng_
	Shared RNG (seeded once, shared across N instances)
double	sigma_factor_ = 0.
	Signed sigma factor for systematic shifts (shift and unisim modes)
int	target_pdg_ = 0
	Target nucleus PDG code for the matched Reaction.
std::vector< double >	varied_
	Pre-generated varied strengths (one per matrix element in the matched Reaction)
Protected Attributes inherited from marley::WeightCalculator
std::string	name_
	Name used to label the output event weight.

Detailed Description

WeightCalculator that varies nuclear matrix element strengths according to their experimental uncertainties. Three modes are available: "multisim" (dimidiated Gaussian random draws), "shift" (deterministic systematic shifts of all matrix elements), and "unisim" (deterministic systematic shifts of a single matrix element at a time).

Definition at line 45 of file StrengthVariationWeightCalculator.hh.

Member Enumeration Documentation

VariationMode

enum class marley::StrengthVariationWeightCalculator::VariationMode

strong

Supported variation modes.

Definition at line 50 of file StrengthVariationWeightCalculator.hh.

{ multisim, shift, unisim };

Member Function Documentation

create_instances()

std::vector< std::shared_ptr< marley::StrengthVariationWeightCalculator > > marley::StrengthVariationWeightCalculator::create_instances ( const marley::JSON & config, marley::Generator & gen )

static

Static factory: validates JSON configuration and creates all variation instances. Called by the Weighter for each strength_variation entry in the weights config array.

Parameters

config	JSON configuration for this calculator
gen	Generator whose loaded reactions are used by the unisim mode to determine the number of matrix elements

Returns

Vector of shared pointers to the created instances

Definition at line 127 of file StrengthVariationWeightCalculator.cc.

{
  // Resolve the reaction file
  if ( !config.has_key( "reaction_file" ) ) {
    throw marley::Error( "Missing \"reaction_file\" key in a"
      " strength_variation weight calculator JSON configuration" );
  }
  std::string reaction_file = config.at( "reaction_file" ).to_string();
  std::string resolved_reaction_file = marley::FileManager::Instance()
    .find_file( reaction_file );
  if ( resolved_reaction_file.empty() ) {
    throw marley::Error( "Could not find reaction data file \""
      + reaction_file + "\" requested by a strength_variation"
      " weight calculator" );
  }
 
  // Read the optional variation mode (default "multisim")
  std::string mode_str = "multisim";
  if ( config.has_key( "mode" ) ) {
    mode_str = config.at( "mode" ).to_string();
  }
 
  // Extract the base name
  if ( !config.has_key( "name" ) ) {
    throw marley::Error( "Missing \"name\" key in a"
      " strength_variation weight calculator JSON configuration" );
  }
  std::string base_name = config.at( "name" ).to_string();
 
  std::vector< std::shared_ptr<
    StrengthVariationWeightCalculator > > instances;
 
  if ( mode_str == "multisim" ) {
 
    // Reject sigma_factor key in multisim mode
    if ( config.has_key( "sigma_factor" ) ) {
      throw marley::Error( "The \"sigma_factor\" key is not allowed"
        " in multisim mode for strength_variation weight calculators" );
    }
 
    // Read and validate the number of variations
    if ( !config.has_key( "num_variations" ) ) {
      throw marley::Error( "Missing \"num_variations\" key in a"
        " strength_variation weight calculator JSON configuration" );
    }
    const auto& nv = config.at( "num_variations" );
    if ( !nv.is_integer() ) {
      throw marley::Error( "The \"num_variations\" value must be a"
        " positive integer" );
    }
    long num_instances = nv.to_long();
    if ( num_instances <= 0 ) {
      throw marley::Error( "The \"num_variations\" value must be a"
        " positive integer" );
    }
 
    // Read the optional seed (default 0)
    long seed = 0;
    if ( config.has_key( "seed" ) ) {
      seed = config.at( "seed" ).to_long();
    }
 
    // Create one shared RNG for all instances
    auto rng = std::make_shared< std::mt19937_64 >(
      static_cast< std::mt19937_64::result_type >( seed ) );
 
    // Create num_instances weight calculators, each sharing the RNG
    for ( long idx = 0; idx < num_instances; ++idx ) {
      instances.push_back( std::shared_ptr<
        StrengthVariationWeightCalculator >(
          new StrengthVariationWeightCalculator(
            base_name + '_' + std::to_string( idx ),
            rng, resolved_reaction_file ) ) );
    }
  }
  else if ( mode_str == "shift" ) {
 
    // Reject keys inappropriate for shift
    if ( config.has_key( "num_variations" ) ) {
      throw marley::Error( "The \"num_variations\" key is not allowed"
        " in shift mode for strength_variation weight"
        " calculators" );
    }
    if ( config.has_key( "seed" ) ) {
      throw marley::Error( "The \"seed\" key is not allowed"
        " in shift mode for strength_variation weight"
        " calculators" );
    }
 
    // Read and validate sigma_factor
    if ( !config.has_key( "sigma_factor" ) ) {
      throw marley::Error( "Missing \"sigma_factor\" key in a"
        " strength_variation weight calculator JSON configuration"
        " with shift mode" );
    }
    const auto& sf = config.at( "sigma_factor" );
 
    std::vector< double > factors;
    if ( sf.is_array() ) {
      for ( const auto& elem : sf.array_range() ) {
        factors.push_back( elem.to_double_or_throw() );
      }
      if ( factors.empty() ) {
        throw marley::Error( "The \"sigma_factor\" array must have"
          " at least one element" );
      }
    }
    else {
      factors.push_back( sf.to_double_or_throw() );
    }
 
    // Reject duplicate sigma_factor values (exact equality)
    {
      std::set< double > seen;
      for ( double f : factors ) {
        if ( !seen.insert( f ).second ) {
          throw marley::Error( "Duplicate sigma_factor value \""
            + std::to_string( f ) + "\" in strength_variation"
            " weight calculator configuration" );
        }
      }
    }
 
    // Generate shortest roundtrip labels for all factors
    auto labels = shortest_roundtrip_labels( factors );
 
    // Create two instances per factor (+k and -k)
    for ( size_t i = 0; i < factors.size(); ++i ) {
      double k = factors[ i ];
 
      // +k "up" instance
      instances.push_back( std::shared_ptr<
        StrengthVariationWeightCalculator >(
          new StrengthVariationWeightCalculator(
            base_name + "-up@" + labels[ i ],
            +k, resolved_reaction_file ) ) );
 
      // -k "down" instance
      instances.push_back( std::shared_ptr<
        StrengthVariationWeightCalculator >(
          new StrengthVariationWeightCalculator(
            base_name + "-down@" + labels[ i ],
            -k, resolved_reaction_file ) ) );
    }
  }
  else if ( mode_str == "unisim" ) {
 
    // Reject keys inappropriate for unisim
    if ( config.has_key( "num_variations" ) ) {
      throw marley::Error( "The \"num_variations\" key is not allowed"
        " in unisim mode for strength_variation weight calculators" );
    }
    if ( config.has_key( "seed" ) ) {
      throw marley::Error( "The \"seed\" key is not allowed"
        " in unisim mode for strength_variation weight calculators" );
    }
 
    // Read and validate sigma_factor
    if ( !config.has_key( "sigma_factor" ) ) {
      throw marley::Error( "Missing \"sigma_factor\" key in a"
        " strength_variation weight calculator JSON configuration"
        " with unisim mode" );
    }
    const auto& sf = config.at( "sigma_factor" );
 
    std::vector< double > factors;
    if ( sf.is_array() ) {
      for ( const auto& elem : sf.array_range() ) {
        factors.push_back( elem.to_double_or_throw() );
      }
      if ( factors.empty() ) {
        throw marley::Error( "The \"sigma_factor\" array must have"
          " at least one element" );
      }
    }
    else {
      factors.push_back( sf.to_double_or_throw() );
    }
 
    // Reject duplicate sigma_factor values (exact equality)
    {
      std::set< double > seen;
      for ( double f : factors ) {
        if ( !seen.insert( f ).second ) {
          throw marley::Error( "Duplicate sigma_factor value \""
            + std::to_string( f ) + "\" in strength_variation"
            " weight calculator configuration" );
        }
      }
    }
 
    // Find the matching DiscreteNuclearReaction from the Generator
    const DiscreteNuclearReaction* dnr = nullptr;
    for ( const auto& rptr : gen.get_reactions() ) {
      if ( rptr->source_file() != resolved_reaction_file ) continue;
      auto pt = rptr->process_type();
      if ( pt != Reaction::ProcessType::NeutrinoCC_Discrete
        && pt != Reaction::ProcessType::AntiNeutrinoCC_Discrete
        && pt != Reaction::ProcessType::NC_Discrete )
      {
        continue;
      }
      dnr = dynamic_cast< const DiscreteNuclearReaction* >(
        rptr.get() );
      if ( !dnr ) throw marley::Error( "Reaction data file \""
        + resolved_reaction_file + "\" was loaded as a discrete nuclear"
        " reaction type but the corresponding Reaction object is not a"
        " DiscreteNuclearReaction. This should not happen and likely"
        " indicates a bug in MARLEY." );
      break;
    }
    if ( !dnr ) throw marley::Error( "Could not find a discrete nuclear"
      " reaction with the source file \"" + resolved_reaction_file
      + "\" in the Generator." );
 
    size_t M = dnr->matrix_elements().size();
 
    // Generate shortest roundtrip labels for all factors
    auto labels = shortest_roundtrip_labels( factors );
 
    // Create two instances per factor per matrix element
    for ( size_t i = 0; i < factors.size(); ++i ) {
      double k = factors[ i ];
      for ( size_t m = 0; m < M; ++m ) {
        // +k "up" instance for this matrix element
        instances.push_back( std::shared_ptr<
          StrengthVariationWeightCalculator >(
            new StrengthVariationWeightCalculator(
              base_name + "-me" + std::to_string( m )
                + "_up@" + labels[ i ],
              +k, m, resolved_reaction_file ) ) );
 
        // -k "down" instance for this matrix element
        instances.push_back( std::shared_ptr<
          StrengthVariationWeightCalculator >(
            new StrengthVariationWeightCalculator(
              base_name + "-me" + std::to_string( m )
                + "_down@" + labels[ i ],
              -k, m, resolved_reaction_file ) ) );
      }
    }
  }
  else {
    throw marley::Error( "Unrecognized variation mode \""
      + mode_str + "\" for strength_variation weight calculator."
      " Allowed values are \"multisim\", \"shift\","
      " and \"unisim\"" );
  }
 
  return instances;
}

References marley::Reaction::AntiNeutrinoCC_Discrete, marley::FileManager::find_file(), marley::Generator::get_reactions(), marley::FileManager::Instance(), marley::DiscreteNuclearReaction::matrix_elements(), marley::Reaction::NC_Discrete, and marley::Reaction::NeutrinoCC_Discrete.

ensure_initialized()

void marley::StrengthVariationWeightCalculator::ensure_initialized ( marley::Generator & gen ) const

protected

Lazy initialization of per-instance varied strengths.

Definition at line 380 of file StrengthVariationWeightCalculator.cc.

{
  if ( initialized_ ) return;
 
  // Scan the Generator's reactions to find the one matching our
  // resolved reaction file path and a discrete nuclear process type
  for ( const auto& reaction_ptr : gen.get_reactions() ) {
 
    if ( reaction_ptr->source_file() != resolved_reaction_file_ ) continue;
 
    // Only nuclear reactions populating discrete nuclear levels may be
    // handled by this weight calculator
    auto pt = reaction_ptr->process_type();
    if ( pt != Reaction::ProcessType::NeutrinoCC_Discrete
      && pt != Reaction::ProcessType::AntiNeutrinoCC_Discrete
      && pt != Reaction::ProcessType::NC_Discrete )
    {
      continue;
    }
 
    auto* dnr = dynamic_cast< const DiscreteNuclearReaction* >(
      reaction_ptr.get() );
    if ( !dnr ) throw marley::Error( "Reaction data file \""
      + resolved_reaction_file_ + "\" was loaded as a discrete nuclear"
      " reaction type but the corresponding Reaction object is not a"
      " DiscreteNuclearReaction. This should not happen and likely"
      " indicates a bug in MARLEY." );
 
    dnr_ = dnr;
    process_type_ = pt;
    target_pdg_ = dnr->pdg_b();
 
    // Pre-generate varied strengths depending on the variation mode
    const auto& matrix_els = dnr->matrix_elements();
    varied_.reserve( matrix_els.size() );
 
    if ( mode_ == VariationMode::multisim ) {
      std::normal_distribution< double > normal_dist;
      for ( const auto& me : matrix_els ) {
        double nom = me.strength();
        double err_low = me.strength_err_low();
        double err_high = me.strength_err_high();
 
        double varied;
        if ( err_low == 0. && err_high == 0. ) {
          varied = nom;
        }
        else {
          double u = normal_dist( *rng_ );
          double sigma = ( u >= 0. ) ? err_high : err_low;
          varied = nom + u * sigma;
          if ( varied < 0. ) varied = 0.;
        }
        varied_.push_back( varied );
      }
    }
    else if ( mode_ == VariationMode::shift ) {
      for ( const auto& me : matrix_els ) {
        double nom = me.strength();
        double err_low = me.strength_err_low();
        double err_high = me.strength_err_high();
 
        double varied;
        if ( err_low == 0. && err_high == 0. ) {
          varied = nom;
        }
        else {
          double sigma = ( sigma_factor_ >= 0. ) ? err_high : err_low;
          varied = nom + sigma_factor_ * sigma;
          if ( varied < 0. ) varied = 0.;
        }
        varied_.push_back( varied );
      }
    }
    else { // unisim
      for ( const auto& me : matrix_els ) {
        varied_.push_back( me.strength() );
      }
      const auto& target_me = matrix_els[ me_idx_ ];
      double nom = target_me.strength();
      double err_low = target_me.strength_err_low();
      double err_high = target_me.strength_err_high();
      if ( err_low != 0. || err_high != 0. ) {
        double sigma = ( sigma_factor_ >= 0. ) ? err_high : err_low;
        double varied = nom + sigma_factor_ * sigma;
        if ( varied < 0. ) varied = 0.;
        varied_[ me_idx_ ] = varied;
      }
    }
 
    initialized_ = true;
    return;
  }
 
  throw marley::Error( "Could not find a discrete nuclear reaction"
    " with the source file \"" + resolved_reaction_file_
    + "\" in the Generator." );
}

References marley::Reaction::AntiNeutrinoCC_Discrete, dnr_, marley::Generator::get_reactions(), initialized_, me_idx_, mode_, marley::Reaction::NC_Discrete, marley::Reaction::NeutrinoCC_Discrete, process_type_, resolved_reaction_file_, rng_, sigma_factor_, target_pdg_, and varied_.

Referenced by weight().

weight()

double marley::StrengthVariationWeightCalculator::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 480 of file StrengthVariationWeightCalculator.cc.

{
  ensure_initialized( gen );
 
  // Check that the event's process type matches this calculator's
  auto sp_attr = event.attribute< HepMC3::IntAttribute >(
    "signal_process_id" );
  if ( !sp_attr ) return 1.;
  auto event_pt = marley_hepmc3::from_nuhepmc_proc_id(
    sp_attr->value() );
  if ( event_pt != process_type_ ) return 1.;
 
  // Verify that the event's target nucleus PDG code matches
  auto target = marley_hepmc3::get_target( event );
  if ( !target || target->pdg_id() != target_pdg_ ) return 1.;
 
  // Read the matrix element index from the event
  auto mi_attr = event.attribute< HepMC3::IntAttribute >(
    "me_index" );
  if ( !mi_attr ) return 1.;
  size_t mi = static_cast< size_t >( mi_attr->value() );
  if ( mi >= varied_.size() ) return 1.;
 
  // Compute the weight as the ratio of the varied strength to the
  // nominal strength
  double nom = dnr_->matrix_elements().at( mi ).strength();
  if ( nom == 0. ) return 1.;
  return varied_.at( mi ) / nom;
}

References dnr_, ensure_initialized(), process_type_, target_pdg_, and varied_.

Member Data Documentation

dnr_

const DiscreteNuclearReaction* marley::StrengthVariationWeightCalculator::dnr_ = nullptr

mutableprotected

Pointer to the matched DiscreteNuclearReaction.

Definition at line 112 of file StrengthVariationWeightCalculator.hh.

Referenced by ensure_initialized(), and weight().

initialized_

bool marley::StrengthVariationWeightCalculator::initialized_ = false

mutableprotected

Whether lazy initialization has been completed.

Definition at line 109 of file StrengthVariationWeightCalculator.hh.

Referenced by ensure_initialized().

me_idx_

size_t marley::StrengthVariationWeightCalculator::me_idx_ = 0

protected

Index of the single matrix element to vary (unisim mode)

Definition at line 106 of file StrengthVariationWeightCalculator.hh.

Referenced by ensure_initialized().

mode_

VariationMode marley::StrengthVariationWeightCalculator::mode_

protected

Variation mode for this instance.

Definition at line 99 of file StrengthVariationWeightCalculator.hh.

Referenced by ensure_initialized().

process_type_

Reaction::ProcessType marley::StrengthVariationWeightCalculator::process_type_ = Reaction::ProcessType::Unknown

mutableprotected

Process type of the matched Reaction.

Definition at line 115 of file StrengthVariationWeightCalculator.hh.

Referenced by ensure_initialized(), and weight().

resolved_reaction_file_

std::string marley::StrengthVariationWeightCalculator::resolved_reaction_file_

protected

Resolved path of the reaction input file of interest.

Definition at line 96 of file StrengthVariationWeightCalculator.hh.

Referenced by ensure_initialized().

rng_

std::shared_ptr< std::mt19937_64 > marley::StrengthVariationWeightCalculator::rng_

protected

Shared RNG (seeded once, shared across N instances)

Definition at line 93 of file StrengthVariationWeightCalculator.hh.

Referenced by ensure_initialized().

sigma_factor_

double marley::StrengthVariationWeightCalculator::sigma_factor_ = 0.

protected

Signed sigma factor for systematic shifts (shift and unisim modes)

Definition at line 103 of file StrengthVariationWeightCalculator.hh.

Referenced by ensure_initialized().

target_pdg_

int marley::StrengthVariationWeightCalculator::target_pdg_ = 0

mutableprotected

Target nucleus PDG code for the matched Reaction.

Definition at line 119 of file StrengthVariationWeightCalculator.hh.

Referenced by ensure_initialized(), and weight().

varied_

std::vector< double > marley::StrengthVariationWeightCalculator::varied_

mutableprotected

Pre-generated varied strengths (one per matrix element in the matched Reaction)

Definition at line 123 of file StrengthVariationWeightCalculator.hh.

Referenced by ensure_initialized(), and weight().

---

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

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