diff --git a/docs/source/usersguide/tallies.rst b/docs/source/usersguide/tallies.rst index 75ca142ba..69e8506d9 100644 --- a/docs/source/usersguide/tallies.rst +++ b/docs/source/usersguide/tallies.rst @@ -108,9 +108,11 @@ The following tables show all valid scores: +----------------------+---------------------------------------------------+ |Score | Description | +======================+===================================================+ - |absorption |Total absorption rate. This accounts for all | - | |reactions which do not produce secondary neutrons | - | |as well as fission. | + |absorption |Total absorption rate. For incident neutrons, this | + | |accounts for all reactions that do not produce | + | |secondary neutrons as well as fission. For incident| + | |photons, this includes photoelectric and pair | + | |production. | +----------------------+---------------------------------------------------+ |elastic |Elastic scattering reaction rate. | +----------------------+---------------------------------------------------+ @@ -194,6 +196,14 @@ The following tables show all valid scores: +----------------------+---------------------------------------------------+ |(n,da) |(n,d\ :math:`\alpha`\ ) reaction rate. | +----------------------+---------------------------------------------------+ + |coherent-scatter |Coherent (Rayleigh) scattering reaction rate. | + +----------------------+---------------------------------------------------+ + |incoherent-scatter |Incoherent (Compton) scattering reaction rate. | + +----------------------+---------------------------------------------------+ + |photoelectric |Photoelectric absorption reaction rate. | + +----------------------+---------------------------------------------------+ + |pair-production |Pair production reaction rate. | + +----------------------+---------------------------------------------------+ |*Arbitrary integer* |An arbitrary integer is interpreted to mean the | | |reaction rate for a reaction with a given ENDF MT | | |number. | diff --git a/src/reaction.cpp b/src/reaction.cpp index fbe2c34de..b3559b10c 100644 --- a/src/reaction.cpp +++ b/src/reaction.cpp @@ -159,11 +159,11 @@ const std::unordered_map REACTION_NAME_MAP { {N_XA, "(n,Xa)"}, {HEATING, "heating"}, {DAMAGE_ENERGY, "damage-energy"}, - {COHERENT, "coherent scatter"}, - {INCOHERENT, "incoherent scatter"}, - {PAIR_PROD_ELEC, "pair production, electron"}, - {PAIR_PROD, "pair production"}, - {PAIR_PROD_NUC, "pair production, nuclear"}, + {COHERENT, "coherent-scatter"}, + {INCOHERENT, "incoherent-scatter"}, + {PAIR_PROD_ELEC, "pair-production-electron"}, + {PAIR_PROD, "pair-production"}, + {PAIR_PROD_NUC, "pair-production-nuclear"}, {PHOTOELECTRIC, "photoelectric"}, {N_PC, "(n,pc)"}, {N_DC, "(n,dc)"}, diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 924975358..c025755d1 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -216,6 +216,14 @@ score_str_to_int(std::string score_str) return N_XA; if (score_str == "damage-energy") return DAMAGE_ENERGY; + if (score_str == "coherent-scatter") + return COHERENT; + if (score_str == "incoherent-scatter") + return INCOHERENT; + if (score_str == "pair-production") + return PAIR_PROD; + if (score_str == "photoelectric") + return PHOTOELECTRIC; // So far we have not identified this score string. Check to see if it is a // deprecated score. diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 5c971d1cf..3ca005cfa 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -530,6 +530,8 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, // Get the pre-collision energy of the particle. auto E = p.E_last_; + using Type = Particle::Type; + for (auto i = 0; i < tally.scores_.size(); ++i) { auto score_bin = tally.scores_[i]; auto score_index = start_index + i; @@ -549,8 +551,7 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, score = p.wgt_last_; } - if (p.type_ == Particle::Type::neutron || - p.type_ == Particle::Type::photon) { + if (p.type_ == Type::neutron || p.type_ == Type::photon) { score *= flux / p.macro_xs_.total; } else { score = 0.; @@ -575,9 +576,9 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, } else { if (i_nuclide >= 0) { - if (p.type_ == Particle::Type::neutron) { + if (p.type_ == Type::neutron) { score = p.neutron_xs_[i_nuclide].total * atom_density * flux; - } else if (p.type_ == Particle::Type::photon) { + } else if (p.type_ == Type::photon) { score = p.photon_xs_[i_nuclide].total * atom_density * flux; } } else { @@ -588,7 +589,7 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_INVERSE_VELOCITY: - if (p.type_ != Particle::Type::neutron) continue; + if (p.type_ != Type::neutron) continue; if (tally.estimator_ == TallyEstimator::ANALOG) { // All events score to an inverse velocity bin. We actually use a @@ -611,7 +612,7 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_SCATTER: - if (p.type_ != Particle::Type::neutron) continue; + if (p.type_ != Type::neutron && p.type_ != Type::photon) continue; if (tally.estimator_ == TallyEstimator::ANALOG) { // Skip any event where the particle didn't scatter @@ -621,17 +622,26 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, score = p.wgt_last_ * flux; } else { if (i_nuclide >= 0) { - score = (p.neutron_xs_[i_nuclide].total - - p.neutron_xs_[i_nuclide].absorption) * atom_density * flux; + if (p.type_ == Type::neutron) { + const auto& micro = p.neutron_xs_[i_nuclide]; + score = (micro.total - micro.absorption) * atom_density * flux; + } else { + const auto& micro = p.photon_xs_[i_nuclide]; + score = (micro.coherent + micro.incoherent) * atom_density * flux; + } } else { - score = (p.macro_xs_.total - p.macro_xs_.absorption) * flux; + if (p.type_ == Type::neutron) { + score = (p.macro_xs_.total - p.macro_xs_.absorption) * flux; + } else { + score = (p.macro_xs_.coherent + p.macro_xs_.incoherent) * flux; + } } } break; case SCORE_NU_SCATTER: - if (p.type_ != Particle::Type::neutron) continue; + if (p.type_ != Type::neutron) continue; // Only analog estimators are available. // Skip any event where the particle didn't scatter @@ -654,7 +664,7 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case SCORE_ABSORPTION: - if (p.type_ != Particle::Type::neutron) continue; + if (p.type_ != Type::neutron && p.type_ != Type::photon) continue; if (tally.estimator_ == TallyEstimator::ANALOG) { if (settings::survival_biasing) { @@ -670,10 +680,18 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, } } else { if (i_nuclide >= 0) { - score = p.neutron_xs_[i_nuclide].absorption * atom_density - * flux; + if (p.type_ == Type::neutron) { + score = p.neutron_xs_[i_nuclide].absorption * atom_density * flux; + } else { + const auto& xs = p.photon_xs_[i_nuclide]; + score = (xs.total - xs.coherent - xs.incoherent) * atom_density * flux; + } } else { - score = p.macro_xs_.absorption * flux; + if (p.type_ == Type::neutron) { + score = p.macro_xs_.absorption * flux; + } else { + score = (p.macro_xs_.photoelectric + p.macro_xs_.pair_production) * flux; + } } } break; @@ -1212,7 +1230,7 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case ELASTIC: - if (p.type_ != Particle::Type::neutron) continue; + if (p.type_ != Type::neutron) continue; if (tally.estimator_ == TallyEstimator::ANALOG) { // Check if event MT matches @@ -1253,7 +1271,7 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, case N_GAMMA: case N_P: case N_A: - if (p.type_ != Particle::Type::neutron) continue; + if (p.type_ != Type::neutron) continue; if (tally.estimator_ == TallyEstimator::ANALOG) { // Check if the event MT matches @@ -1288,8 +1306,46 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, break; + case COHERENT: + case INCOHERENT: + case PHOTOELECTRIC: + case PAIR_PROD: + if (p.type_ != Type::photon) continue; + + if (tally.estimator_ == TallyEstimator::ANALOG) { + if (score_bin == PHOTOELECTRIC) { + // Photoelectric events are assigned an MT value corresponding to the + // shell cross section. Also, photons below the energy cutoff are + // assumed to have been absorbed via photoelectric absorption + if ((p.event_mt_ < 534 || p.event_mt_ > 572) && + p.event_mt_ != REACTION_NONE) continue; + } else { + if (p.event_mt_ != score_bin) continue; + } + score = p.wgt_last_ * flux; + } else { + if (i_nuclide >= 0) { + const auto& micro = p.photon_xs_[i_nuclide]; + double xs = + (score_bin == COHERENT) ? micro.coherent : + (score_bin == INCOHERENT) ? micro.incoherent : + (score_bin == PHOTOELECTRIC) ? micro.photoelectric : + micro.pair_production; + score = xs * atom_density * flux; + } else { + double xs = + (score_bin == COHERENT) ? p.macro_xs_.coherent : + (score_bin == INCOHERENT) ? p.macro_xs_.incoherent : + (score_bin == PHOTOELECTRIC) ? p.macro_xs_.photoelectric : + p.macro_xs_.pair_production; + score = xs * flux; + } + } + break; + + case HEATING: - if (p.type_ == Particle::Type::neutron) { + if (p.type_ == Type::neutron) { score = score_neutron_heating(p, tally, flux, HEATING, i_nuclide, atom_density); } else { @@ -1312,7 +1368,7 @@ score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, // The default block is really only meant for redundant neutron reactions // (e.g. 444, 901) - if (p.type_ != Particle::Type::neutron) continue; + if (p.type_ != Type::neutron) continue; if (tally.estimator_ == TallyEstimator::ANALOG) {