Merge pull request #14 from njoy/refactor/reaction

Introducing partial reaction identifiers for summation reactions

python/src/Reaction.python.cpp

@@ -32,22 +32,42 @@ void wrapReaction( python::module& module, python::module& ) {
   component
   .def(
 
-    python::init< ReactionID, ReactionType, TabulatedCrossSection,
+    python::init< ReactionID, TabulatedCrossSection,
                   std::vector< ReactionProduct >,
                   std::optional< double >, std::optional< double > >(),
-    python::arg( "id" ), python::arg( "type" ), python::arg( "xs" ),
+    python::arg( "id" ), python::arg( "xs" ),
     python::arg( "products" ) = std::vector< ReactionProduct >{},
     python::arg( "mass_q" ) = std::nullopt,
     python::arg( "reaction_q" ) = std::nullopt,
-    "Initialise the reaction\n\n"
+    "Initialise a primary reaction\n\n"
     "Arguments:\n"
     "    self         the reaction\n"
     "    id           the reaction identifier\n"
-    "    type         the reaction type\n"
     "    xs           the cross section of the reaction\n"
     "    products     the reaction products\n"
     "    mass_q       the mass difference Q value (optional)\n"
-    "    reaction_q   the reaction Q value (optional)\n"
+    "    reaction_q   the reaction Q value (optional)"
+  )
+  .def(
+
+    python::init< ReactionID,
+                  std::vector< ReactionID >,
+                  TabulatedCrossSection,
+                  std::vector< ReactionProduct >,
+                  std::optional< double >, std::optional< double > >(),
+    python::arg( "id" ), python::arg( "partials" ), python::arg( "xs" ),
+    python::arg( "products" ) = std::vector< ReactionProduct >{},
+    python::arg( "mass_q" ) = std::nullopt,
+    python::arg( "reaction_q" ) = std::nullopt,
+    "Initialise a summation reaction\n\n"
+    "Arguments:\n"
+    "    self        the reaction\n"
+    "    id          the reaction identifier\n"
+    "    partials    the identifiers of the partials of the reaction\n"
+    "    xs          the cross section of the reaction\n"
+    "    products     the reaction products\n"
+    "    mass_q       the mass difference Q value (optional)\n"
+    "    reaction_q   the reaction Q value (optional)"
   )
   .def_property_readonly(
 
@@ -61,6 +81,13 @@ void wrapReaction( python::module& module, python::module& ) {
     &Component::type,
     "The reaction type"
   )
+  .def_property_readonly(
+
+    "partial_reaction_identifiers",
+    &Component::partialReactionIdentifiers,
+    "The summation reaction identifiers (not defined if this is a primary\n"
+    "reaction)"
+  )
   .def_property_readonly(
 
     "mass_difference_qvalue",

python/test/Test_dryad_ProjectileTarget.py

@@ -147,12 +147,12 @@ def verify_linearised_chunk( self, chunk ) :
         # the data is given explicitly
         chunk = ProjectileTarget(
                     projectile = "n", target = "Fe56", type = InteractionType.Nuclear,
-                    reactions = [ Reaction( "n,Fe56->n,Fe56", ReactionType.Primary,
+                    reactions = [ Reaction( "n,Fe56->n,Fe56",
                                             TabulatedCrossSection( [ 1e-5, 20. ], [ 1000., 10. ],
                                                                    InterpolationType.LogLinear ),
                                             [],
                                             0, 0 ),
-                                  Reaction( "n,Fe56->n,Fe56_e1", ReactionType.Primary,
+                                  Reaction( "n,Fe56->n,Fe56_e1",
                                             TabulatedCrossSection( [ 1., 20. ], [ 0., 100. ],
                                                                    InterpolationType.LinearLinear ),
                                             [],

python/test/Test_dryad_Reaction.py

@@ -27,6 +27,9 @@ def verify_chunk( self, chunk ) :
             # reaction type
             self.assertEqual( ReactionType.Primary, chunk.type )
 
+            # partial identifiers
+            self.assertEqual( None, chunk.partial_reaction_identifiers )
+
             # q values
             self.assertAlmostEqual( 0, chunk.mass_difference_qvalue )
             self.assertAlmostEqual( -1, chunk.reaction_qvalue )
@@ -69,6 +72,11 @@ def verify_summation_chunk( self, chunk ) :
             # reaction type
             self.assertEqual( ReactionType.Summation, chunk.type )
 
+            # partial identifiers
+            self.assertEqual( 2, len( chunk.partial_reaction_identifiers ) )
+            self.assertEqual( 'n,Fe56->elastic', chunk.partial_reaction_identifiers[0] )
+            self.assertEqual( 'n,Fe56->2n,Fe56', chunk.partial_reaction_identifiers[1] )
+
             # q values
             self.assertEqual( None, chunk.mass_difference_qvalue )
             self.assertEqual( None, chunk.reaction_qvalue )
@@ -142,7 +150,7 @@ def verify_linearised_chunk( self, chunk ) :
             self.assertEqual( True, chunk.cross_section.is_linearised )
 
         # the data is given explicitly
-        chunk = Reaction( id = 'n,Fe56->n,Fe56_e1', type = ReactionType.Primary,
+        chunk = Reaction( id = 'n,Fe56->n,Fe56_e1',
                           mass_q = 0, reaction_q = -1,
                           xs = TabulatedCrossSection ( [ 1., 2., 2., 3., 4. ],
                                                        [ 4., 3., 4., 3., 2. ],
@@ -165,7 +173,8 @@ def verify_linearised_chunk( self, chunk ) :
         verify_linearised_chunk( self, copy )
 
         # the data is given explicitly for a summation reaction
-        chunk = Reaction( id = 'n,Fe56->total', type = ReactionType.Summation,
+        chunk = Reaction( id = 'n,Fe56->total',
+                          partials = [ 'n,Fe56->elastic', 'n,Fe56->2n,Fe56' ],
                           xs = TabulatedCrossSection ( [ 1., 2., 2., 3., 4. ],
                                                        [ 4., 3., 4., 3., 2. ],
                                                        [ 1, 4 ],

src/dryad/ProjectileTarget/test/ProjectileTarget.test.cpp

@@ -27,12 +27,12 @@ SCENARIO( "ProjectileTarget" ) {
 
       std::vector< Reaction > reactions = {
 
-        Reaction( id::ReactionID( "n,Fe56->n,Fe56" ), ReactionType::Primary,
+        Reaction( id::ReactionID( "n,Fe56->n,Fe56" ),
                   TabulatedCrossSection( { 1e-5, 20. }, { 1000., 10. },
                                            InterpolationType::LogLinear ),
                   {},
                   0, 0 ),
-        Reaction( id::ReactionID( "n,Fe56->n,Fe56_e1" ), ReactionType::Primary,
+        Reaction( id::ReactionID( "n,Fe56->n,Fe56_e1" ),
                   TabulatedCrossSection( { 1., 20. }, { 0., 100. },
                                            InterpolationType::LinearLinear ),
                   {},

src/dryad/Reaction.hpp

@@ -26,6 +26,7 @@ namespace dryad {
     id::ReactionID id_;
 
     ReactionType type_;
+    std::optional< std::vector< id::ReactionID > > partials_;
 
     std::optional< double > mass_difference_qvalue_;
     std::optional< double > reaction_qvalue_;
@@ -58,6 +59,16 @@ namespace dryad {
       return this->type_;
     }
 
+    /**
+     *  @brief Return the summation reaction identifiers (not defined if this is
+     *         a primary reaction)
+     */
+    const std::optional< std::vector< id::ReactionID > >&
+    partialReactionIdentifiers() const noexcept {
+
+      return this->partials_;
+    }
+
     /**
      *  @brief Return the mass difference Q value
      *
@@ -119,10 +130,18 @@ namespace dryad {
     Reaction linearise( ToleranceConvergence tolerance = {} ) const noexcept {
 
       TabulatedCrossSection xs = this->crossSection().linearise( tolerance );
-
-      return Reaction( this->identifier(), this->type(), std::move( xs ),
-                       this->products(), this->massDifferenceQValue(),
-                       this->reactionQValue() );
+      if ( this->type() == ReactionType::Primary ) {
+
+        return Reaction( this->identifier(), std::move( xs ),
+                         this->products(), this->massDifferenceQValue(),
+                         this->reactionQValue() );
+      }
+      else {
+
+        return Reaction( this->identifier(),
+                         this->partialReactionIdentifiers().value(),
+                         std::move( xs ) );
+      }
     }
 
     /**

src/dryad/Reaction/src/ctor.hpp

@@ -12,12 +12,14 @@
  */
 Reaction( id::ReactionID&& id,
           ReactionType&& type,
+          std::optional< std::vector< id::ReactionID > >&& partials,
           TabulatedCrossSection&& xs,
           std::vector< ReactionProduct >&& products,
           std::optional< double >&& mass_q,
           std::optional< double >&& reaction_q,
           bool linearised ) :
     id_( std::move( id ) ), type_( std::move( type ) ),
+    partials_( std::move( partials ) ),
     mass_difference_qvalue_( mass_q ),
     reaction_qvalue_( reaction_q ),
     xs_( std::move( xs ) ),
@@ -38,7 +40,7 @@ Reaction& operator=( const Reaction& ) = default;
 Reaction& operator=( Reaction&& ) = default;
 
 /**
- *  @brief Constructor
+ *  @brief Constructor for primary reactions
  *
  *  @param id           the reaction identifier
  *  @param xs           the cross section of the reaction
@@ -47,15 +49,37 @@ Reaction& operator=( Reaction&& ) = default;
  *  @param reaction_q   the reaction Q value
  */
 Reaction( id::ReactionID id,
-          ReactionType type,
           TabulatedCrossSection xs,
           std::vector< ReactionProduct > products = {},
           std::optional< double > mass_q = std::nullopt,
           std::optional< double > reaction_q = std::nullopt ) :
     Reaction( std::move( id ),
-              std::move( type ),
+              ReactionType::Primary,
+              std::nullopt,
               std::move( xs ),
               std::move( products ),
               std::move( mass_q ),
               std::move( reaction_q ),
               xs.isLinearised() ? true : false ) {}
+
+/**
+ *  @brief Constructor for summation reactions
+ *
+ *  @param id          the reaction identifier
+ *  @param xs          the cross section of the reaction
+ *  @param partials    the identifiers of the partials of the reaction
+ */
+Reaction( id::ReactionID id,
+          std::vector< id::ReactionID > partials,
+          TabulatedCrossSection xs,
+          std::vector< ReactionProduct > products = {},
+          std::optional< double > mass_q = std::nullopt,
+          std::optional< double > reaction_q = std::nullopt ) :
+    Reaction( std::move( id ),
+              ReactionType::Summation,
+              std::move( partials ),
+              std::move( xs ),
+              std::move( products ),
+              std::nullopt,
+              std::nullopt,
+              xs.isLinearised() ? true : false ) {}

src/dryad/Reaction/test/Reaction.test.cpp

@@ -22,7 +22,6 @@ SCENARIO( "Reaction" ) {
     WHEN( "the data is given explicitly" ) {
 
       id::ReactionID id( "n,Fe56->n,Fe56_e1" );
-      ReactionType type = ReactionType::Primary;
       TabulatedCrossSection xs( { 1., 2., 2., 3., 4. },
                                 { 4., 3., 4., 3., 2. },
                                 { 1, 4 },
@@ -32,7 +31,7 @@ SCENARIO( "Reaction" ) {
       double reaction_q = -1;
       std::vector< ReactionProduct > products = {};
 
-      Reaction chunk( std::move( id ), std::move( type ), std::move( xs ),
+      Reaction chunk( std::move( id ), std::move( xs ),
                       std::move( products ), mass_q, reaction_q );
 
       THEN( "a Reaction can be constructed and members can be tested" ) {
@@ -62,14 +61,14 @@ SCENARIO( "Reaction" ) {
     WHEN( "the data is given explicitly" ) {
 
       id::ReactionID id( "n,Fe56->total" );
-      ReactionType type = ReactionType::Summation;
+      std::vector< id::ReactionID > partials = { "n,Fe56->elastic", "n,Fe56->2n,Fe55" };
       TabulatedCrossSection xs( { 1., 2., 2., 3., 4. },
                                 { 4., 3., 4., 3., 2. },
                                 { 1, 4 },
                                 { InterpolationType::LinearLinear,
                                   InterpolationType::LinearLog } );
 
-      Reaction chunk( std::move( id ), std::move( type ), std::move( xs ) );
+      Reaction chunk( std::move( id ), std::move( partials ), std::move( xs ) );
 
       THEN( "a Reaction can be constructed and members can be tested" ) {
 
@@ -102,6 +101,9 @@ void verifyChunk( const Reaction& chunk ) {
   // reaction type
   CHECK( ReactionType::Primary == chunk.type() );
 
+  // partial identifiers
+  CHECK( std::nullopt == chunk.partialReactionIdentifiers() );
+
   // q values
   CHECK_THAT( 0, WithinRel( chunk.massDifferenceQValue().value() ) );
   CHECK_THAT( -1, WithinRel( chunk.reactionQValue().value() ) );
@@ -185,6 +187,13 @@ void verifySummationChunk( const Reaction& chunk ) {
   // reaction type
   CHECK( ReactionType::Summation == chunk.type() );
 
+  // partial identifiers
+  CHECK( std::nullopt != chunk.partialReactionIdentifiers() );
+  auto partials = chunk.partialReactionIdentifiers().value();
+  CHECK( 2 == partials.size() );
+  CHECK( id::ReactionID( "n,Fe56->elastic" ) == partials[0] );
+  CHECK( id::ReactionID( "n,Fe56->2n,Fe55" ) == partials[1] );
+
   // q values
   CHECK( std::nullopt == chunk.massDifferenceQValue() );
   CHECK( std::nullopt == chunk.reactionQValue() );

src/dryad/ReactionType.hpp

@@ -13,12 +13,20 @@ namespace dryad {
  *  @brief The reaction types
  *
  *  This enum is used to differentiate reaction types in the ProjectileTarget.
- *  We currently distinguish two types of reactions: primary and summation.
+ *  We currently distinguish four types of reactions: primary, summation,
+ *  production and incomplete
  */
 enum class ReactionType : short {
 
-  Primary = 1,   /**< The reaction is a primary independent reaction */
-  Summation = 2  /**< The reaction is a pure summation reaction */
+  /**
+   *  A primary independent reaction that contributes to the total cross section
+   */
+  Primary = 1,
+  /**
+   *  A summation reaction with or without reaction products that does not count
+   *  towards the total cross section
+   */
+  Summation = 2
 };
 
 } // dryad namespace

src/dryad/format/endf.hpp

@@ -6,6 +6,7 @@
 #include "dryad/format/endf/createTargetIdentifier.hpp"
 #include "dryad/format/endf/createProductIdentifier.hpp"
 #include "dryad/format/endf/createInteractionType.hpp"
+#include "dryad/format/endf/createReactionType.hpp"
 
 #include "dryad/format/endf/createTabulatedMultiplicity.hpp"
 #include "dryad/format/endf/createMultiplicity.hpp"