More updates

src/dryad/format/gnds/createMultiplicity.hpp

@@ -19,10 +19,14 @@ namespace gnds {
   /**
    *  @brief Create an integer or tabulated multiplicity from a parsed ENDF multiplicity
    */
-  std::variant< int, TabulatedMultiplicity > createMultiplicity( pugi::xml_node multiplicity ) {
+  std::variant< int, TabulatedMultiplicity > 
+  createMultiplicity( pugi::xml_node multiplicity, const std::string& style = "eval" ) {
 
+    // check that this is a valid multiplicity node
     throwExceptionOnWrongNode( multiplicity, "multiplicity" );
-    auto child = multiplicity.first_child();
+
+    // check the first child with the requested style and act accordingly
+    auto child = multiplicity.find_child_by_attribute( "label", style.c_str() );
     if ( strcmp( child.name(), "constant1d" ) == 0 ) {
 
       auto data = readConstant1dAsInteger( child );

src/dryad/format/gnds/createQValue.hpp

@@ -0,0 +1,48 @@
+#ifndef NJOY_DRYAD_FORMAT_GNDS_CREATEQVALUE
+#define NJOY_DRYAD_FORMAT_GNDS_CREATEQVALUE
+
+// system includes
+#include <variant>
+
+// other includes
+#include "pugixml.hpp"
+#include "tools/Log.hpp"
+#include "dryad/format/gnds/throwExceptionOnWrongNode.hpp"
+#include "dryad/format/gnds/createTabulatedMultiplicity.hpp"
+#include "dryad/format/gnds/readConstant1d.hpp"
+#include "dryad/format/gnds/convertEnergy.hpp"
+
+namespace njoy {
+namespace dryad {
+namespace format {
+namespace gnds {
+
+  /**
+   *  @brief Create a Q value from a GNDS q node
+   */
+  double createQValue( pugi::xml_node q, const std::string& style = "eval" ) {
+
+    // check that this is a valid q node
+    throwExceptionOnWrongNode( q, "Q" );
+
+    // check the first child with the requested style and act accordingly
+    auto child = q.find_child_by_attribute( "label", style.c_str() );
+    if ( strcmp( child.name(), "constant1d" ) == 0 ) {
+
+      auto data = readConstant1dAsDouble( child );
+      convertEnergy( data.first, data.second );
+      return data.first;
+    }
+    else {
+
+      Log::error( "The GNDS node named \'{}\' does not define supported q value data", q.name() );
+      throw std::exception();
+    }
+  }
+
+} // gnds namespace
+} // format namespace
+} // dryad namespace
+} // njoy namespace
+
+#endif

src/dryad/format/gnds/createReaction.hpp

@@ -8,7 +8,7 @@
 #include "pugixml.hpp"
 #include "tools/Log.hpp"
 #include "dryad/format/gnds/convertEnergy.hpp"
-#include "dryad/format/gnds/readConstant1d.hpp"
+#include "dryad/format/gnds/createQValue.hpp"
 #include "dryad/format/gnds/createReactionProducts.hpp"
 #include "dryad/format/gnds/createTabulatedCrossSection.hpp"
 #include "dryad/Reaction.hpp"
@@ -22,7 +22,8 @@ namespace gnds {
    *  @brief Create a Reaction from an unparsed ENDF material
    */
   Reaction createReaction( const id::ParticleID& projectile, const id::ParticleID& target,
-                           pugi::xml_node suite, pugi::xml_node reaction ) {
+                           pugi::xml_node suite, pugi::xml_node reaction, 
+                           const std::string& style = "eval" ) {
 
     if ( strcmp( reaction.name(), "reaction" ) == 0 ) {
 
@@ -31,14 +32,12 @@ namespace gnds {
 
       // cross section
       auto section = reaction.child( "crossSection" );
-      TabulatedCrossSection xs = createTabulatedCrossSection( section );
+      TabulatedCrossSection xs = createTabulatedCrossSection( section, style );
 
       // Q values
       auto output = reaction.child( "outputChannel" );
-      auto q = readConstant1d( output.child( "Q" ).child( "constant1d" ) );
-      convertEnergy( q.first, q.second );
       std::optional< double > mass_q = std::nullopt;
-      std::optional< double > reaction_q = q.first;
+      std::optional< double > reaction_q = createQValue( output.child( "Q" ), style );
 
       // reaction products
       std::vector< ReactionProduct > products = 

src/dryad/format/gnds/createReactionProduct.hpp

@@ -10,6 +10,7 @@
 #include "dryad/format/gnds/throwExceptionOnWrongNode.hpp"
 #include "dryad/format/gnds/createMultiplicity.hpp"
 #include "dryad/format/gnds/createTwoBodyDistributionData.hpp"
+#include "dryad/format/gnds/createUncorrelatedDistributionData.hpp"
 #include "dryad/ReactionProduct.hpp"
 
 namespace njoy {
@@ -22,31 +23,50 @@ namespace gnds {
    */
   ReactionProduct
   createReactionProduct( const id::ParticleID& projectile, const id::ParticleID& target,
-                         pugi::xml_node suite, pugi::xml_node product ) {
+                         pugi::xml_node suite, pugi::xml_node product,
+                         const std::string& style = "eval" ) {
 
+    // check that this is a valid product node
     throwExceptionOnWrongNode( product, "product" );
+
+    // get the secondary particle identifier and adjust as required
     std::string pid( product.attribute( "pid" ).as_string() );
     if ( pid == "photon" ) {
 
       pid = "g";
     }
     id::ParticleID id( pid );
-    auto multiplicity = createMultiplicity( product.child( "multiplicity" ) );
 
+    // create the multiplicity
+    auto multiplicity = createMultiplicity( product.child( "multiplicity" ), style );
+
+    // get distribution
     auto distribution = product.child( "distribution" );
     if ( distribution ) {
 
-      auto first = distribution.first_child();
-      if ( strcmp( first.name(), "angularTwoBody" ) == 0 ) {
+      // get the first node of the requested style and act accordingly
+      auto node = distribution.find_child_by_attribute( "label", style.c_str() );
+      if ( strcmp( node.name(), "angularTwoBody" ) == 0 ) {
 
-        auto recoil = first.child( "recoil" );
+        // ignore recoil distributions for now
+        auto recoil = node.child( "recoil" );
         if ( ! recoil ) {
 
-          return ReactionProduct( id, multiplicity, createTwoBodyDistributionData( first ) );
+          return ReactionProduct( id, multiplicity, createTwoBodyDistributionData( node ) );
         }
       }
+      if ( strcmp( node.name(), "uncorrelated" ) == 0 ) {
+
+        return ReactionProduct( id, multiplicity, createUncorrelatedDistributionData( node ) );
+      }
+      else {
+
+        // for now, return a basic reaction product
+        return ReactionProduct( id, multiplicity );
+      }
     }
-
+
+    // return a basic reaction product
     return ReactionProduct( id, multiplicity );
   }
 

src/dryad/format/gnds/createReactionProducts.hpp

@@ -20,16 +20,21 @@ namespace gnds {
    */
   std::vector< ReactionProduct >
   createReactionProducts( const id::ParticleID& projectile, const id::ParticleID& target,
-                          pugi::xml_node suite, pugi::xml_node data ) {
+                          pugi::xml_node suite, pugi::xml_node products,
+                         const std::string& style = "eval" ) {
 
-    std::vector< ReactionProduct > products;
-    for ( pugi::xml_node product = data.child( "product" ); product; 
+    // check that this is a valid products node
+    throwExceptionOnWrongNode( products, "products" );
+
+    // loop over product children
+    std::vector< ReactionProduct > data;
+    for ( pugi::xml_node product = products.child( "product" ); product; 
           product = product.next_sibling( "product" ) ) {
 
-      products.emplace_back( createReactionProduct( projectile, target, suite, product ) );
+      data.emplace_back( createReactionProduct( projectile, target, suite, product, style ) );
     }
 
-    return products;
+    return data;
   }
 
 } // gnds namespace

src/dryad/format/gnds/createReactions.hpp

@@ -20,37 +20,47 @@ namespace gnds {
    */
   std::vector< Reaction >
   createReactions( const id::ParticleID& projectile, const id::ParticleID& target,
-                   pugi::xml_node suite ) {
+                   pugi::xml_node suite, const std::string& style = "eval" ) {
+
+    // check that this is a valid reaction node
+    throwExceptionOnWrongNode( suite, "reactionSuite" );
 
     std::vector< Reaction > reactions;
 
+    // get the children that contain the reaction data
     pugi::xml_node primaries = suite.child( "reactions" );
     pugi::xml_node sums = suite.child( "sums" ).child( "crossSectionSums" );
     pugi::xml_node incomplete = suite.child( "incompleteReactions" );
 
+    // treat primaries
     if ( primaries ) {
 
+      // loop over reaction nodes
       for ( pugi::xml_node reaction = primaries.child( "reaction" ); 
             reaction; reaction = reaction.next_sibling( "reaction" ) ) {
 
         Log::info( "Reading data for MT{}", reaction.attribute( "ENDF_MT" ).as_string() );
-        reactions.emplace_back( createReaction( projectile, target, suite, reaction ) );
+        reactions.emplace_back( createReaction( projectile, target, suite, reaction, style ) );
       }
     }
 
+    // treat summation
     if ( sums ) {
 
+      // loop over crossSectionSum nodes
       for ( pugi::xml_node reaction = sums.child( "crossSectionSum" ); 
             reaction; reaction = reaction.next_sibling( "crossSectionSum" ) ) {
 
         Log::info( "Reading data for MT{}", reaction.attribute( "ENDF_MT" ).as_string() );
-        reactions.emplace_back( createReaction( projectile, target, suite, reaction ) );
+        reactions.emplace_back( createReaction( projectile, target, suite, reaction, style ) );
       }
     }
 
+    // treat incomplete reactions
     if ( incomplete ) {
 
-      for ( pugi::xml_node reaction = incomplete.child( "reaction" ); 
+      // loop over reaction nodes
+     for ( pugi::xml_node reaction = incomplete.child( "reaction" ); 
             reaction; reaction = reaction.next_sibling( "reaction" ) ) {
 
         Log::info( "Reading data for MT{}", reaction.attribute( "ENDF_MT" ).as_string() );

src/dryad/format/gnds/createTabulatedCrossSection.hpp

@@ -21,7 +21,9 @@ namespace gnds {
   /**
    *  @brief Create a TabulatedCrossSection from a parsed ENDF section
    */
-  TabulatedCrossSection createTabulatedCrossSection( const pugi::xml_node& xs ) {
+  TabulatedCrossSection 
+  createTabulatedCrossSection( const pugi::xml_node& xs, 
+                               const std::string& style = "eval" ) {
 
     Log::info( "Reading cross section data" );
 
@@ -30,11 +32,14 @@ namespace gnds {
     std::vector< std::size_t > boundaries;
     std::vector< InterpolationType > interpolants;
 
-    auto xys1d = xs.child( "XYs1d" );
-    if ( xys1d ) {
+    // check that this is a valid q node
+    throwExceptionOnWrongNode( xs, "crossSection" );
+
+    auto node = xs.find_child_by_attribute( "label", style.c_str() );
+    if ( strcmp( node.name(), "XYs1d" ) == 0 ) {
 
       // read the cross section data
-      auto data = readXYs1D( xys1d );
+      auto data = readXYs1D( node );
 
       // get the interpolation type
       auto interpolant = createInterpolationType( std::get< 6 >( data ) );
@@ -49,52 +54,50 @@ namespace gnds {
       boundaries.emplace_back( energies.size() - 1 );
       interpolants.emplace_back( interpolant );
     }
-    else {
+    else if ( strcmp( node.name(), "regions1d" ) == 0 ) {
 
-      auto regions1d = xs.child( "regions1d" );
-      if ( regions1d ) {
+      // get the units
+      auto units = readAxes( node.child( "axes" ) );
 
-        auto units = readAxes( regions1d.child( "axes" ) );
+      // loop over the children of function1ds
+      pugi::xml_node function1ds = node.child( "function1ds" );
+      for ( pugi::xml_node xys1d = function1ds.child( "XYs1d" ); 
+            xys1d; xys1d = xys1d.next_sibling(  "XYs1d"  ) ) {
 
-        pugi::xml_node function1ds = regions1d.child( "function1ds" );
-        for ( pugi::xml_node xys1d = function1ds.child( "XYs1d" ); 
-              xys1d; xys1d = xys1d.next_sibling(  "XYs1d"  ) ) {
+        // read the current interpolation region
+        auto data = readXYs1D( xys1d, units );
 
-          // read the current interpolation region
-          auto data = readXYs1D( xys1d, units );
+        // get the interpolation type
+        auto interpolant = createInterpolationType( std::get< 6 >( data ) );
 
-          // get the interpolation type
-          auto interpolant = createInterpolationType( std::get< 6 >( data ) );
+        // convert units - if necessary
+        convertEnergies( std::get< 2 >( data ), std::get< 3 >( data ) );
+        convertCrossSections( std::get< 4 >( data ), std::get< 5 >( data ) );
 
-          // convert units - if necessary
-          convertEnergies( std::get< 2 >( data ), std::get< 3 >( data ) );
-          convertCrossSections( std::get< 4 >( data ), std::get< 5 >( data ) );
+        // check for duplicate points at interpolation region boundaries
+        std::size_t offset = 0;
+        if ( energies.size() > 0 ) {
 
-          // check for duplicate points at interpolation region boundaries
-          std::size_t offset = 0;
-          if ( energies.size() > 0 ) {
+          if ( energies.back() == std::get< 2 >( data ).front() &&
+               values.back() == std::get< 4 >( data ).front() ) {
 
-            if ( energies.back() == std::get< 2 >( data ).front() &&
-                 values.back() == std::get< 4 >( data ).front() ) {
-
-              offset = 1;
-            }
+            offset = 1;
           }
-
-          // grow the data accordingly
-          energies.insert( energies.end(), std::get< 2 >( data ).begin() + offset, std::get< 2 >( data ).end() );
-          values.insert( values.end(), std::get< 4 >( data ).begin() + offset, std::get< 4 >( data ).end() );
-          boundaries.emplace_back( energies.size() - 1 );
-          interpolants.emplace_back( interpolant );
         }
-      }
-      else {
 
-        Log::error( "Expected either an XYs1d node or regions1d node with XYs1d nodes"
-                    "for cross section data" );
-        throw std::exception();
+        // grow the data accordingly
+        energies.insert( energies.end(), std::get< 2 >( data ).begin() + offset, std::get< 2 >( data ).end() );
+        values.insert( values.end(), std::get< 4 >( data ).begin() + offset, std::get< 4 >( data ).end() );
+        boundaries.emplace_back( energies.size() - 1 );
+        interpolants.emplace_back( interpolant );
       }
     }
+    else {
+
+      Log::error( "Expected either an XYs1d node or regions1d node with XYs1d nodes"
+                  "for cross section data" );
+      throw std::exception();
+    }
 
     return TabulatedCrossSection(
              std::move( energies ), std::move( values ),