diff --git a/python/src/Reaction.python.cpp b/python/src/Reaction.python.cpp index 31d74df..c31f9c4 100644 --- a/python/src/Reaction.python.cpp +++ b/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", diff --git a/python/test/Test_dryad_ProjectileTarget.py b/python/test/Test_dryad_ProjectileTarget.py index 56c4ab9..6540a3c 100644 --- a/python/test/Test_dryad_ProjectileTarget.py +++ b/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 ), [], diff --git a/python/test/Test_dryad_Reaction.py b/python/test/Test_dryad_Reaction.py index c54dc80..610ef8c 100644 --- a/python/test/Test_dryad_Reaction.py +++ b/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 ], diff --git a/src/dryad/ProjectileTarget/test/ProjectileTarget.test.cpp b/src/dryad/ProjectileTarget/test/ProjectileTarget.test.cpp index 7be1793..5422c67 100644 --- a/src/dryad/ProjectileTarget/test/ProjectileTarget.test.cpp +++ b/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 ), {}, diff --git a/src/dryad/Reaction.hpp b/src/dryad/Reaction.hpp index 9397e5e..ff19763 100644 --- a/src/dryad/Reaction.hpp +++ b/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 ) ); + } } /** diff --git a/src/dryad/Reaction/src/ctor.hpp b/src/dryad/Reaction/src/ctor.hpp index 3ed4f02..c23080c 100644 --- a/src/dryad/Reaction/src/ctor.hpp +++ b/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 ) {} diff --git a/src/dryad/Reaction/test/Reaction.test.cpp b/src/dryad/Reaction/test/Reaction.test.cpp index 8861c1b..d65b615 100644 --- a/src/dryad/Reaction/test/Reaction.test.cpp +++ b/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() ); diff --git a/src/dryad/ReactionType.hpp b/src/dryad/ReactionType.hpp index 1141de4..729981e 100644 --- a/src/dryad/ReactionType.hpp +++ b/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 diff --git a/src/dryad/format/endf.hpp b/src/dryad/format/endf.hpp index dbe36c0..d62cee1 100644 --- a/src/dryad/format/endf.hpp +++ b/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" diff --git a/src/dryad/format/endf/ReactionIdentifiers.hpp b/src/dryad/format/endf/ReactionIdentifiers.hpp index a8f48a1..9bb7c05 100644 --- a/src/dryad/format/endf/ReactionIdentifiers.hpp +++ b/src/dryad/format/endf/ReactionIdentifiers.hpp @@ -2,6 +2,7 @@ #define NJOY_DRYAD_FORMAT_ENDF_REACTIONIDENTIFIERS // system includes +#include #include // other includes @@ -24,11 +25,11 @@ namespace endf { }; inline static const std::unordered_set< int > summation_or_primary_ = { - 16, 18, 103, 104, 105, 106, 107 + 16, 18, 103, 104, 105, 106, 107, 526 }; inline static const std::unordered_set< int > summation_ = { - 1, 3, 4, 27, 101, 501, 516, 522, 526 + 1, 3, 4, 27, 101, 501, 516, 522 }; inline static const std::unordered_set< int > primary_ = { @@ -78,13 +79,61 @@ namespace endf { 875, 876, 877, 878, 879, 880, 881, 882, 883, 884, 885, 886, 887, 888, 889, 890, 891 }; + inline static const std::map< int, std::vector< int > > partials_ = { + + //! @todo complete this list + { 1, { } }, + { 3, { } }, + { 4, { } }, + { 16, { } }, + { 18, { 19, 20, 21, 38 } }, + { 27, { } }, + { 101, { } }, + { 103, { 600, 601, 602, 603, 604, 605, 606, 607, 608, 609, + 610, 611, 612, 613, 614, 615, 616, 617, 618, 619, + 620, 621, 622, 623, 624, 625, 626, 627, 628, 629, + 630, 631, 632, 633, 634, 635, 636, 637, 638, 639, + 640, 641, 642, 643, 644, 645, 646, 647, 648, 649 } }, + { 104, { 650, 651, 652, 653, 654, 655, 656, 657, 658, 659, + 660, 661, 662, 663, 664, 665, 666, 667, 668, 669, + 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, + 680, 681, 682, 683, 684, 685, 686, 687, 688, 689, + 690, 691, 692, 693, 694, 695, 696, 697, 698, 699 } }, + { 105, { 700, 701, 702, 703, 704, 705, 706, 707, 708, 709, + 710, 711, 712, 713, 714, 715, 716, 717, 718, 719, + 720, 721, 722, 723, 724, 725, 726, 727, 728, 729, + 730, 731, 732, 733, 734, 735, 736, 737, 738, 739, + 740, 741, 742, 743, 744, 745, 746, 747, 748, 749 } }, + { 106, { 750, 751, 752, 753, 754, 755, 756, 757, 758, 759, + 760, 761, 762, 763, 764, 765, 766, 767, 768, 769, + 770, 771, 772, 773, 774, 775, 776, 777, 778, 779, + 780, 781, 782, 783, 784, 785, 786, 787, 788, 789, + 790, 791, 792, 793, 794, 795, 796, 797, 798, 799 } }, + { 107, { 800, 801, 802, 803, 804, 805, 806, 807, 808, 809, + 810, 811, 812, 813, 814, 815, 816, 817, 818, 819, + 820, 821, 822, 823, 824, 825, 826, 827, 828, 829, + 830, 831, 832, 833, 834, 835, 836, 837, 838, 839, + 840, 841, 842, 843, 844, 845, 846, 847, 848, 849 } }, + { 501, { 502, 504, 515, 517, 525, 527, 528, 534, 535, + 536, 537, 538, 539, 540, 541, 542, 543, 544, 545, + 546, 547, 548, 549, 550, 551, 552, 553, 554, 555, + 556, 557, 558, 559, 560, 561, 562, 563, 564, 565, + 566, 567, 568, 569, 570, 571, 572 } }, + { 516, { 515, 517 } }, + { 522, { 534, 535, 536, 537, 538, 539, 540, 541, 542, 543, + 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, + 554, 555, 556, 557, 558, 559, 560, 561, 562, 563, + 564, 565, 566, 567, 568, 569, 570, 571, 572 } }, + { 526, { 525 } } + }; public: /** - * @brief Return whther or not the MT number is valid + * @brief Return whether or not the MT number is valid * - * @param[in] mt the MT number + * @param[in] material the ENDF material + * @param[in] mt the MT number */ static bool isValid( const ENDFtk::tree::Material&, int mt ) { @@ -105,9 +154,10 @@ namespace endf { } /** - * @brief Return whther or not the MT number is for a derived reaction + * @brief Return whether or not the MT number is for a derived reaction * - * @param[in] mt the MT number + * @param[in] material the ENDF material + * @param[in] mt the MT number */ static bool isDerived( const ENDFtk::tree::Material&, int mt ) { @@ -115,21 +165,23 @@ namespace endf { } /** - * @brief Return whther or not the MT number is for a summation reaction + * @brief Return whether or not the MT number is for a summation reaction * - * @param[in] mt the MT number + * @param[in] material the ENDF material + * @param[in] mt the MT number */ static bool isSummation( const ENDFtk::tree::Material& material, int mt ) { switch ( mt ) { - case 16 : return material.hasSection( 3, 875 ) ? true : false; - case 18 : return material.hasSection( 3, 19 ) ? true : false; - case 103 : return material.hasSection( 3, 600 ) ? true : false; - case 104 : return material.hasSection( 3, 650 ) ? true : false; - case 105 : return material.hasSection( 3, 700 ) ? true : false; - case 106 : return material.hasSection( 3, 750 ) ? true : false; - case 107 : return material.hasSection( 3, 800 ) ? true : false; + case 16 : return material.hasSection( 3, 875 ) ? true : false; + case 18 : return material.hasSection( 3, 19 ) ? true : false; + case 103 : return material.hasSection( 3, 600 ) ? true : false; + case 104 : return material.hasSection( 3, 650 ) ? true : false; + case 105 : return material.hasSection( 3, 700 ) ? true : false; + case 106 : return material.hasSection( 3, 750 ) ? true : false; + case 107 : return material.hasSection( 3, 800 ) ? true : false; + case 526 : return material.hasSection( 23, 525 ) ? true : false; default : { return summation_.find( mt ) != summation_.end(); @@ -138,27 +190,59 @@ namespace endf { } /** - * @brief Return whther or not the MT number is for a primary reaction + * @brief Return whether or not the MT number is for a primary reaction * - * @param[in] mt the MT number + * @param[in] material the ENDF material + * @param[in] mt the MT number */ static bool isPrimary( const ENDFtk::tree::Material& material, int mt ) { switch ( mt ) { - case 16 : return material.hasSection( 3, 875 ) ? false : true; - case 18 : return material.hasSection( 3, 19 ) ? false : true; - case 103 : return material.hasSection( 3, 600 ) ? false : true; - case 104 : return material.hasSection( 3, 650 ) ? false : true; - case 105 : return material.hasSection( 3, 700 ) ? false : true; - case 106 : return material.hasSection( 3, 750 ) ? false : true; - case 107 : return material.hasSection( 3, 800 ) ? false : true; + case 16 : return material.hasSection( 3, 875 ) ? false : true; + case 18 : return material.hasSection( 3, 19 ) ? false : true; + case 103 : return material.hasSection( 3, 600 ) ? false : true; + case 104 : return material.hasSection( 3, 650 ) ? false : true; + case 105 : return material.hasSection( 3, 700 ) ? false : true; + case 106 : return material.hasSection( 3, 750 ) ? false : true; + case 107 : return material.hasSection( 3, 800 ) ? false : true; + case 526 : return material.hasSection( 23, 525 ) ? false : true; default : { return primary_.find( mt ) != primary_.end(); } } } + + /** + * @brief Return the partial mt numbers for a summation mt number + * + * @param[in] material the ENDF material + * @param[in] mt the MT number + */ + static std::vector< int > partials( const ENDFtk::tree::Material& material, + int mf, int mt ) { + + std::vector< int > partials; + if ( isSummation( material, mt ) ) { + + // create a vector of partial mt numbers that are present + std::copy_if( partials_.at( mt ).begin(), partials_.at( mt ).end(), + std::back_inserter( partials ), + [&material, mf] ( auto&& value ) { return material.hasSection( mf, value ); } ); + + // add deficiency mt number for total elastic in electro-atomic data + auto iter = std::lower_bound( partials.begin(), partials.end(), 525 ); + if ( iter != partials.end() ) { + + if ( *iter == 525 ) { + + partials.insert( iter + 1, -526 ); + } + } + } + return partials; + } }; } // endf namespace diff --git a/src/dryad/format/endf/createReaction.hpp b/src/dryad/format/endf/createReaction.hpp index b2a5354..af9b6ef 100644 --- a/src/dryad/format/endf/createReaction.hpp +++ b/src/dryad/format/endf/createReaction.hpp @@ -26,62 +26,94 @@ namespace endf { if ( material.hasSection( 3, mt ) ) { - auto section = material.section( 3, mt ).parse< 3 >(); - // metadata and miscellaneous information id::ReactionID id = std::to_string( mt ); ReactionType type = createReactionType( material, mt ); - std::optional< double > mass_q = std::nullopt; - std::optional< double > reaction_q = std::nullopt; - if ( type == ReactionType::Primary ) { - - mass_q = section.massDifferenceQValue(); - reaction_q = section.reactionQValue(); - } // cross section + auto section = material.section( 3, mt ).parse< 3 >(); TabulatedCrossSection xs = createTabulatedCrossSection( section ); + // Q values + std::optional< double > mass_q = std::nullopt; + std::optional< double > reaction_q = std::nullopt; + // reaction products - std::vector< ReactionProduct > products = {}; + std::vector< ReactionProduct > products = createReactionProducts( projectile, target, material, mt ); + if ( type == ReactionType::Primary ) { - products = createReactionProducts( projectile, target, material, mt ); + // Q values + if ( mt == 18 ) { + + //! @todo handle fission Q value defined on MT18 when no partials are defined + } + else { + + mass_q = section.massDifferenceQValue(); + reaction_q = section.massDifferenceQValue(); + } + + // return the reaction data + return Reaction( std::move( id ), std::move( xs ), + std::move( products ), std::move( mass_q ), + std::move( reaction_q ) ); } + else { + + if ( mt == 18 ) { + + //! @todo handle fission Q value defined on MT18 if partials are defined + } + + // partials + std::vector< int > endf_partials = ReactionIdentifiers::partials( material, 3, mt ); + std::vector< id::ReactionID > partials( endf_partials.size() ); + std::transform( endf_partials.begin(), endf_partials.end(), partials.begin(), + [] ( auto&& value ) { return std::to_string( value ); } ); - // return the reaction data - return Reaction( std::move( id ), std::move( type ), std::move( xs ), - std::move( products ), std::move( mass_q ), - std::move( reaction_q ) ); + // return the reaction data + return Reaction( std::move( id ), std::move( partials ), std::move( xs ) ); + } } else if ( material.hasSection( 23, mt ) ) { - auto section = material.section( 23, mt ).parse< 23 >(); - // metadata and miscellaneous information id::ReactionID id = std::to_string( mt ); ReactionType type = createReactionType( material, mt ); - std::optional< double > mass_q = std::nullopt; - std::optional< double > reaction_q = std::nullopt; - if ( type == ReactionType::Primary ) { - - reaction_q = -section.subshellBindingEnergy(); - } // cross section + auto section = material.section( 23, mt ).parse< 23 >(); TabulatedCrossSection xs = createTabulatedCrossSection( section ); + // q values + std::optional< double > mass_q = std::nullopt; + std::optional< double > reaction_q = std::nullopt; + // reaction products - std::vector< ReactionProduct > products = {}; + std::vector< ReactionProduct > products = createReactionProducts( projectile, target, material, mt ); + if ( type == ReactionType::Primary ) { - products = createReactionProducts( projectile, target, material, mt ); + // q values + reaction_q = -section.subshellBindingEnergy(); + + // return the reaction data + return Reaction( std::move( id ), std::move( xs ), + std::move( products ), std::move( mass_q ), + std::move( reaction_q ) ); } + else { - // return the reaction data - return Reaction( std::move( id ), std::move( type ), std::move( xs ), - std::move( products ), std::move( mass_q ), - std::move( reaction_q ) ); + // partials + std::vector< int > endf_partials = ReactionIdentifiers::partials( material, 23, mt ); + std::vector< id::ReactionID > partials( endf_partials.size() ); + std::transform( endf_partials.begin(), endf_partials.end(), partials.begin(), + [] ( auto&& value ) { return std::to_string( value ); } ); + + // return the reaction data + return Reaction( std::move( id ), std::move( partials ), std::move( xs ) ); + } } else { diff --git a/src/dryad/format/endf/createReactions.hpp b/src/dryad/format/endf/createReactions.hpp index 3d3efb4..4448f88 100644 --- a/src/dryad/format/endf/createReactions.hpp +++ b/src/dryad/format/endf/createReactions.hpp @@ -32,6 +32,8 @@ namespace endf { if ( material.hasFile( 3 ) || material.hasFile( 23 ) ) { int source = material.hasFile( 3 ) ? 3 : 23; + + // loop over the available reactions and create the reaction objects reactions.reserve( material.file( source ).sectionNumbers().size() ); for ( auto mt : material.file( source ).sectionNumbers() ) { @@ -45,6 +47,17 @@ namespace endf { Log::warning( "Skipping data for derived MT{}", mt ); } } + + // calculate deficit reaction for elastic scattering in electro-atomic data + if ( material.hasSection( 23, 526 ) && ReactionIdentifiers::isSummation( material, 526 ) ) { + + auto total = material.section( 23, 526 ).parse< 23 >(); + TabulatedCrossSection deficit = createTabulatedCrossSection( total ).linearise(); + auto partial = material.section( 23, 525 ).parse< 23 >(); + deficit -= createTabulatedCrossSection( partial ).linearise(); + + reactions.push_back( Reaction( "-526", deficit, {}, std::nullopt, 0. ) ); + } } return reactions; } diff --git a/src/dryad/format/endf/test/createProjectileTarget.test.cpp b/src/dryad/format/endf/test/createProjectileTarget.test.cpp index 594fd9f..8b7fac1 100644 --- a/src/dryad/format/endf/test/createProjectileTarget.test.cpp +++ b/src/dryad/format/endf/test/createProjectileTarget.test.cpp @@ -42,6 +42,8 @@ SCENARIO( "createProjectileTarget" ) { CHECK( true == H1.hasReaction( id::ReactionID( "102" ) ) ); CHECK( false == H1.hasReaction( id::ReactionID( "some unknown reaction" ) ) ); + CHECK( 3 == H1.reactions().size() ); + auto total = H1.reactions()[0]; verifyNeutronTotalReaction( total ); @@ -83,15 +85,18 @@ SCENARIO( "createProjectileTarget" ) { CHECK( std::nullopt == H0.resonances() ); - CHECK( true == H0.hasReaction( id::ReactionID( "501" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "522" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "525" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "526" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "527" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "528" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "534" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "501" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "522" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "525" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "526" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "527" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "528" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "534" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "-526" ) ) ); CHECK( false == H0.hasReaction( id::ReactionID( "some unknown reaction" ) ) ); + CHECK( 8 == H0.reactions().size() ); + auto total = H0.reactions()[0]; verifyElectronTotalReaction( total ); @@ -113,6 +118,9 @@ SCENARIO( "createProjectileTarget" ) { auto subionisation = H0.reactions()[6]; verifyElectronSubshellIonisationReaction( subionisation ); + auto deficit = H0.reactions()[7]; + verifyElectronElasticDeficitReaction( deficit ); + total = H0.reaction( id::ReactionID( "501" ) ); verifyElectronTotalReaction( total ); @@ -133,6 +141,9 @@ SCENARIO( "createProjectileTarget" ) { subionisation = H0.reaction( id::ReactionID( "534" ) ); verifyElectronSubshellIonisationReaction( subionisation ); + + deficit = H0.reaction( id::ReactionID( "-526" ) ); + verifyElectronElasticDeficitReaction( deficit ); } // THEN } // WHEN } // GIVEN @@ -157,16 +168,18 @@ SCENARIO( "createProjectileTarget" ) { CHECK( std::nullopt == H0.resonances() ); - CHECK( true == H0.hasReaction( id::ReactionID( "501" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "502" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "504" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "515" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "516" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "517" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "522" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "534" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "501" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "502" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "504" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "515" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "516" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "517" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "522" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "534" ) ) ); CHECK( false == H0.hasReaction( id::ReactionID( "some unknown reaction" ) ) ); + CHECK( 8 == H0.reactions().size() ); + auto total = H0.reactions()[0]; verifyPhotonTotalReaction( total ); diff --git a/src/dryad/format/endf/test/createProjectileTargetFromFile.test.cpp b/src/dryad/format/endf/test/createProjectileTargetFromFile.test.cpp index e851efc..65207a9 100644 --- a/src/dryad/format/endf/test/createProjectileTargetFromFile.test.cpp +++ b/src/dryad/format/endf/test/createProjectileTargetFromFile.test.cpp @@ -39,6 +39,8 @@ SCENARIO( "projectileTarget" ) { CHECK( true == H1.hasReaction( id::ReactionID( "102" ) ) ); CHECK( false == H1.hasReaction( id::ReactionID( "some unknown reaction" ) ) ); + CHECK( 3 == H1.reactions().size() ); + auto total = H1.reactions()[0]; verifyNeutronTotalReaction( total ); @@ -77,15 +79,18 @@ SCENARIO( "projectileTarget" ) { CHECK( std::nullopt == H0.resonances() ); - CHECK( true == H0.hasReaction( id::ReactionID( "501" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "522" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "525" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "526" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "527" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "528" ) ) ); - CHECK( true == H0.hasReaction( id::ReactionID( "534" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "501" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "522" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "525" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "526" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "527" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "528" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "534" ) ) ); + CHECK( true == H0.hasReaction( id::ReactionID( "-526" ) ) ); CHECK( false == H0.hasReaction( id::ReactionID( "some unknown reaction" ) ) ); + CHECK( 8 == H0.reactions().size() ); + auto total = H0.reactions()[0]; verifyElectronTotalReaction( total ); @@ -107,6 +112,9 @@ SCENARIO( "projectileTarget" ) { auto subionisation = H0.reactions()[6]; verifyElectronSubshellIonisationReaction( subionisation ); + auto deficit = H0.reactions()[7]; + verifyElectronElasticDeficitReaction( deficit ); + total = H0.reaction( id::ReactionID( "501" ) ); verifyElectronTotalReaction( total ); @@ -127,6 +135,9 @@ SCENARIO( "projectileTarget" ) { subionisation = H0.reaction( id::ReactionID( "534" ) ); verifyElectronSubshellIonisationReaction( subionisation ); + + deficit = H0.reaction( id::ReactionID( "-526" ) ); + verifyElectronElasticDeficitReaction( deficit ); } // THEN } // WHEN } // GIVEN @@ -158,6 +169,8 @@ SCENARIO( "projectileTarget" ) { CHECK( true == H0.hasReaction( id::ReactionID( "534" ) ) ); CHECK( false == H0.hasReaction( id::ReactionID( "some unknown reaction" ) ) ); + CHECK( 8 == H0.reactions().size() ); + auto total = H0.reactions()[0]; verifyPhotonTotalReaction( total ); diff --git a/src/dryad/format/endf/test/createReactions.test.cpp b/src/dryad/format/endf/test/createReactions.test.cpp index 7635a8e..e3e282a 100644 --- a/src/dryad/format/endf/test/createReactions.test.cpp +++ b/src/dryad/format/endf/test/createReactions.test.cpp @@ -30,6 +30,8 @@ SCENARIO( "createReactions" ) { id::ParticleID target( "H1" ); std::vector< Reaction > reactions = format::endf::createReactions( projectile, target, material ); + CHECK( 3 == reactions.size() ); + auto total = reactions[0]; verifyNeutronTotalReaction( total ); @@ -51,10 +53,12 @@ SCENARIO( "createReactions" ) { THEN( "all reactions can be created" ) { - id::ParticleID projectile( "n" ); - id::ParticleID target( "H1" ); + id::ParticleID projectile( "e-" ); + id::ParticleID target( "H0" ); std::vector< Reaction > reactions = format::endf::createReactions( projectile, target, material ); + CHECK( 8 == reactions.size() ); + auto total = reactions[0]; verifyElectronTotalReaction( total ); @@ -75,6 +79,51 @@ SCENARIO( "createReactions" ) { auto subionisation = reactions[6]; verifyElectronSubshellIonisationReaction( subionisation ); + + auto deficit = reactions[7]; + verifyElectronElasticDeficitReaction( deficit ); + } // THEN + } // WHEN + } // GIVEN + + GIVEN( "ENDF materials - photo-atomic" ) { + + auto tape = njoy::ENDFtk::tree::fromFile( "photoat-001_H_000.endf" ); + auto material = tape.materials().front(); + + WHEN( "a single ENDF material is given" ) { + + THEN( "all reactions can be created" ) { + + id::ParticleID projectile( "G" ); + id::ParticleID target( "H0" ); + std::vector< Reaction > reactions = format::endf::createReactions( projectile, target, material ); + + CHECK( 8 == reactions.size() ); + + auto total = reactions[0]; + verifyPhotonTotalReaction( total ); + + auto coherent = reactions[1]; + verifyPhotonCoherentReaction( coherent ); + + auto incoherent = reactions[2]; + verifyPhotonIncoherentReaction( incoherent ); + + auto epairproduction = reactions[3]; + verifyPhotonElectronFieldPairProductionReaction( epairproduction ); + + auto tpairproduction = reactions[4]; + verifyPhotonTotalPairProductionReaction( tpairproduction ); + + auto npairproduction = reactions[5]; + verifyPhotonNuclearFieldPairProductionReaction( npairproduction ); + + auto tionisation = reactions[6]; + verifyPhotonTotalIonisationReaction( tionisation ); + + auto ionisation = reactions[7]; + verifyPhotonIonisationReaction( ionisation ); } // THEN } // WHEN } // GIVEN diff --git a/src/dryad/format/endf/test/test_verification_functions.hpp b/src/dryad/format/endf/test/test_verification_functions.hpp index ad60b0d..9b55772 100644 --- a/src/dryad/format/endf/test/test_verification_functions.hpp +++ b/src/dryad/format/endf/test/test_verification_functions.hpp @@ -5,6 +5,12 @@ void verifyNeutronTotalReaction( const Reaction& total ) { CHECK( false == total.hasProducts() ); CHECK( false == total.isLinearised() ); + CHECK( std::nullopt != total.partialReactionIdentifiers() ); + auto partials = total.partialReactionIdentifiers().value(); +// CHECK( 2 == partials.size() ); +// CHECK( id::ReactionID( "2" ) == partials[0] ); +// CHECK( id::ReactionID( "102" ) == partials[1] ); + CHECK( std::nullopt == total.massDifferenceQValue() ); CHECK( std::nullopt == total.reactionQValue() ); @@ -102,6 +108,15 @@ void verifyElectronTotalReaction( const Reaction& total ) { CHECK( false == total.hasProducts() ); CHECK( true == total.isLinearised() ); + CHECK( std::nullopt != total.partialReactionIdentifiers() ); + auto partials = total.partialReactionIdentifiers().value(); + CHECK( 5 == partials.size() ); + CHECK( id::ReactionID( "525" ) == partials[0] ); + CHECK( id::ReactionID( "-526" ) == partials[1] ); + CHECK( id::ReactionID( "527" ) == partials[2] ); + CHECK( id::ReactionID( "528" ) == partials[3] ); + CHECK( id::ReactionID( "534" ) == partials[4] ); + CHECK( std::nullopt == total.massDifferenceQValue() ); CHECK( std::nullopt == total.reactionQValue() ); @@ -129,6 +144,11 @@ void verifyElectronTotalIonisationReaction( const Reaction& ionisation ) { CHECK( false == ionisation.hasProducts() ); CHECK( true == ionisation.isLinearised() ); + CHECK( std::nullopt != ionisation.partialReactionIdentifiers() ); + auto partials = ionisation.partialReactionIdentifiers().value(); + CHECK( 1 == partials.size() ); + CHECK( id::ReactionID( "534" ) == partials[0] ); + CHECK( std::nullopt == ionisation.massDifferenceQValue() ); CHECK( std::nullopt == ionisation.reactionQValue() ); @@ -247,6 +267,12 @@ void verifyElectronTotalElasticReaction( const Reaction& telastic ) { CHECK( false == telastic.hasProducts() ); CHECK( true == telastic.isLinearised() ); + CHECK( std::nullopt != telastic.partialReactionIdentifiers() ); + auto partials = telastic.partialReactionIdentifiers().value(); + CHECK( 2 == partials.size() ); + CHECK( id::ReactionID( "525" ) == partials[0] ); + CHECK( id::ReactionID( "-526" ) == partials[1] ); + CHECK( std::nullopt == telastic.massDifferenceQValue() ); CHECK( std::nullopt == telastic.reactionQValue() ); @@ -542,6 +568,15 @@ void verifyPhotonTotalReaction( const Reaction& total ) { CHECK( false == total.hasProducts() ); CHECK( true == total.isLinearised() ); + CHECK( std::nullopt != total.partialReactionIdentifiers() ); + auto partials = total.partialReactionIdentifiers().value(); + CHECK( 5 == partials.size() ); + CHECK( id::ReactionID( "502" ) == partials[0] ); + CHECK( id::ReactionID( "504" ) == partials[1] ); + CHECK( id::ReactionID( "515" ) == partials[2] ); + CHECK( id::ReactionID( "517" ) == partials[3] ); + CHECK( id::ReactionID( "534" ) == partials[4] ); + CHECK( std::nullopt == total.massDifferenceQValue() ); CHECK( std::nullopt == total.reactionQValue() ); @@ -797,6 +832,12 @@ void verifyPhotonTotalPairProductionReaction( const Reaction& tpairproduction ) CHECK( false == tpairproduction.hasProducts() ); CHECK( true == tpairproduction.isLinearised() ); + CHECK( std::nullopt != tpairproduction.partialReactionIdentifiers() ); + auto partials = tpairproduction.partialReactionIdentifiers().value(); + CHECK( 2 == partials.size() ); + CHECK( id::ReactionID( "515" ) == partials[0] ); + CHECK( id::ReactionID( "517" ) == partials[1] ); + CHECK( std::nullopt == tpairproduction.massDifferenceQValue() ); CHECK( std::nullopt == tpairproduction.reactionQValue() ); @@ -852,6 +893,11 @@ void verifyPhotonTotalIonisationReaction( const Reaction& tionisation ) { CHECK( false == tionisation.hasProducts() ); CHECK( true == tionisation.isLinearised() ); + CHECK( std::nullopt != tionisation.partialReactionIdentifiers() ); + auto partials = tionisation.partialReactionIdentifiers().value(); + CHECK( 1 == partials.size() ); + CHECK( id::ReactionID( "534" ) == partials[0] ); + CHECK( std::nullopt == tionisation.massDifferenceQValue() ); CHECK( std::nullopt == tionisation.reactionQValue() ); @@ -871,3 +917,20 @@ void verifyPhotonTotalIonisationReaction( const Reaction& tionisation ) { CHECK( 0 == tionisation.products().size() ); } + +void verifyElectronElasticDeficitReaction( const Reaction& deficit ) { + + CHECK( id::ReactionID( "-526" ) == deficit.identifier() ); + CHECK( ReactionType::Primary == deficit.type() ); + CHECK( false == deficit.hasProducts() ); + CHECK( true == deficit.isLinearised() ); + + CHECK( std::nullopt == deficit.partialReactionIdentifiers() ); + + CHECK( std::nullopt == deficit.massDifferenceQValue() ); + CHECK( std::nullopt != deficit.reactionQValue() ); + + CHECK( true == deficit.crossSection().isLinearised() ); + + CHECK( 0 == deficit.products().size() ); +}