diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 224a9ef14..20329e0c7 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -23,7 +23,7 @@ class FilterBinIter public: //! Construct an iterator over bins that match a given particle's state. - FilterBinIter(const Tally& tally, Particle* p); + FilterBinIter(const Tally& tally, Particle& p); //! Construct an iterator over all filter bin combinations. // @@ -41,7 +41,7 @@ public: int index_ {1}; double weight_ {1.}; - + std::vector& filter_matches_; private: @@ -63,21 +63,21 @@ private: //! since collisions do not occur in voids. // //! \param p The particle being tracked -void score_collision_tally(Particle* p); +void score_collision_tally(Particle& p); //! Score tallies based on a simple count of events (for continuous energy). // //! Analog tallies are triggered at every collision, not every event. // //! \param p The particle being tracked -void score_analog_tally_ce(Particle* p); +void score_analog_tally_ce(Particle& p); //! Score tallies based on a simple count of events (for multigroup). // //! Analog tallies are triggered at every collision, not every event. // //! \param p The particle being tracked -void score_analog_tally_mg(Particle* p); +void score_analog_tally_mg(Particle& p); //! Score tallies using a tracklength estimate of the flux. // @@ -87,13 +87,13 @@ void score_analog_tally_mg(Particle* p); // //! \param p The particle being tracked //! \param distance The distance in [cm] traveled by the particle -void score_tracklength_tally(Particle* p, double distance); +void score_tracklength_tally(Particle& p, double distance); //! Score surface or mesh-surface tallies for particle currents. // //! \param p The particle being tracked //! \param tallies A vector of tallies to score to -void score_surface_tally(Particle* p, const std::vector& tallies); +void score_surface_tally(Particle& p, const std::vector& tallies); } // namespace openmc diff --git a/src/particle.cpp b/src/particle.cpp index 96cf3fdcf..f7c7f3c0e 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -223,7 +223,7 @@ Particle::event_advance() // Score track-length tallies if (!model::active_tracklength_tallies.empty()) { - score_tracklength_tally(this, distance); + score_tracklength_tally(*this, distance); } // Score track-length estimate of k-eff @@ -264,7 +264,7 @@ Particle::event_cross_surface() } // Score cell to cell partial currents if (!model::active_surface_tallies.empty()) { - score_surface_tally(this, model::active_surface_tallies); + score_surface_tally(*this, model::active_surface_tallies); } } @@ -283,7 +283,7 @@ Particle::event_collide() // pre-collision direction to figure out what mesh surfaces were crossed if (!model::active_meshsurf_tallies.empty()) - score_surface_tally(this, model::active_meshsurf_tallies); + score_surface_tally(*this, model::active_meshsurf_tallies); // Clear surface component surface_ = 0; @@ -297,12 +297,12 @@ Particle::event_collide() // Score collision estimator tallies -- this is done after a collision // has occurred rather than before because we need information on the // outgoing energy for any tallies with an outgoing energy filter - if (!model::active_collision_tallies.empty()) score_collision_tally(this); + if (!model::active_collision_tallies.empty()) score_collision_tally(*this); if (!model::active_analog_tallies.empty()) { if (settings::run_CE) { - score_analog_tally_ce(this); + score_analog_tally_ce(*this); } else { - score_analog_tally_mg(this); + score_analog_tally_mg(*this); } } @@ -429,7 +429,7 @@ Particle::cross_surface() // physically moving the particle forward slightly this->r() += TINY_BIT * this->u(); - score_surface_tally(this, model::active_meshsurf_tallies); + score_surface_tally(*this, model::active_meshsurf_tallies); } // Score to global leakage tally @@ -462,14 +462,14 @@ Particle::cross_surface() // with a mesh boundary if (!model::active_surface_tallies.empty()) { - score_surface_tally(this, model::active_surface_tallies); + score_surface_tally(*this, model::active_surface_tallies); } if (!model::active_meshsurf_tallies.empty()) { Position r {this->r()}; this->r() -= TINY_BIT * this->u(); - score_surface_tally(this, model::active_meshsurf_tallies); + score_surface_tally(*this, model::active_meshsurf_tallies); this->r() = r; } @@ -524,7 +524,7 @@ Particle::cross_surface() if (!model::active_meshsurf_tallies.empty()) { Position r {this->r()}; this->r() -= TINY_BIT * this->u(); - score_surface_tally(this, model::active_meshsurf_tallies); + score_surface_tally(*this, model::active_meshsurf_tallies); this->r() = r; } diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 54bf5cec6..78ff6dd1f 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -26,8 +26,8 @@ namespace openmc { // FilterBinIter implementation //============================================================================== -FilterBinIter::FilterBinIter(const Tally& tally, Particle* p) - : filter_matches_{p->filter_matches_}, tally_{tally} +FilterBinIter::FilterBinIter(const Tally& tally, Particle& p) + : filter_matches_{p.filter_matches_}, tally_{tally} { // Find all valid bins in each relevant filter if they have not already been // found for this event. @@ -36,7 +36,7 @@ FilterBinIter::FilterBinIter(const Tally& tally, Particle* p) if (!match.bins_present_) { match.bins_.clear(); match.weights_.clear(); - model::tally_filters[i_filt]->get_all_bins(*p, tally_.estimator_, match); + model::tally_filters[i_filt]->get_all_bins(p, tally_.estimator_, match); match.bins_present_ = true; } @@ -177,15 +177,15 @@ score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index, //! Helper function to retrieve fission q value from a nuclide -double get_nuc_fission_q(const Nuclide& nuc, const Particle* p, int score_bin) +double get_nuc_fission_q(const Nuclide& nuc, const Particle& p, int score_bin) { if (score_bin == SCORE_FISS_Q_PROMPT) { if (nuc.fission_q_prompt_) { - return (*nuc.fission_q_prompt_)(p->E_last_); + return (*nuc.fission_q_prompt_)(p.E_last_); } } else if (score_bin == SCORE_FISS_Q_RECOV) { if (nuc.fission_q_recov_) { - return (*nuc.fission_q_recov_)(p->E_last_); + return (*nuc.fission_q_recov_)(p.E_last_); } } return 0.0; @@ -196,47 +196,47 @@ double get_nuc_fission_q(const Nuclide& nuc, const Particle* p, int score_bin) //! Pulled out to support both the fission_q scores and energy deposition //! score -double score_fission_q(const Particle* p, int score_bin, const Tally& tally, +double score_fission_q(const Particle& p, int score_bin, const Tally& tally, double flux, int i_nuclide, double atom_density) { if (tally.estimator_ == TallyEstimator::ANALOG) { - const Nuclide& nuc {*data::nuclides[p->event_nuclide_]}; + const Nuclide& nuc {*data::nuclides[p.event_nuclide_]}; if (settings::survival_biasing) { // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value - if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { - return p->wgt_absorb_ * get_nuc_fission_q(nuc, p, score_bin) - * p->neutron_xs_[p->event_nuclide_].fission * flux - / p->neutron_xs_[p->event_nuclide_].absorption; + if (p.neutron_xs_[p.event_nuclide_].absorption > 0) { + return p.wgt_absorb_ * get_nuc_fission_q(nuc, p, score_bin) + * p.neutron_xs_[p.event_nuclide_].fission * flux + / p.neutron_xs_[p.event_nuclide_].absorption; } } else { // Skip any non-absorption events - if (p->event_ == TallyEvent::SCATTER) return 0.0; + if (p.event_ == TallyEvent::SCATTER) return 0.0; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { - return p->wgt_last_ * get_nuc_fission_q(nuc, p, score_bin) - * p->neutron_xs_[p->event_nuclide_].fission * flux - / p->neutron_xs_[p->event_nuclide_].absorption; + if (p.neutron_xs_[p.event_nuclide_].absorption > 0) { + return p.wgt_last_ * get_nuc_fission_q(nuc, p, score_bin) + * p.neutron_xs_[p.event_nuclide_].fission * flux + / p.neutron_xs_[p.event_nuclide_].absorption; } } } else { if (i_nuclide >= 0) { const Nuclide& nuc {*data::nuclides[i_nuclide]}; return get_nuc_fission_q(nuc, p, score_bin) * atom_density * flux - * p->neutron_xs_[i_nuclide].fission; + * p.neutron_xs_[i_nuclide].fission; } else { - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; + if (p.material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p.material_]}; double score {0.0}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); const Nuclide& nuc {*data::nuclides[j_nuclide]}; score += get_nuc_fission_q(nuc, p, score_bin) * atom_density - * p->neutron_xs_[j_nuclide].fission; + * p.neutron_xs_[j_nuclide].fission; } return score * flux; } @@ -247,28 +247,28 @@ double score_fission_q(const Particle* p, int score_bin, const Tally& tally, //! Helper function to obtain the kerma coefficient for a given nuclide -double get_nuclide_neutron_heating(const Particle* p, const Nuclide& nuc, +double get_nuclide_neutron_heating(const Particle& p, const Nuclide& nuc, int rxn_index, int i_nuclide) { size_t mt = nuc.reaction_index_[rxn_index]; if (mt == C_NONE) return 0.0; - auto i_temp = p->neutron_xs_[i_nuclide].index_temp; + auto i_temp = p.neutron_xs_[i_nuclide].index_temp; if (i_temp < 0) return 0.0; // Can be true due to multipole const auto& rxn {*nuc.reactions_[mt]}; const auto& xs {rxn.xs_[i_temp]}; - auto i_grid = p->neutron_xs_[i_nuclide].index_grid; + auto i_grid = p.neutron_xs_[i_nuclide].index_grid; if (i_grid < xs.threshold) return 0.0; - auto f = p->neutron_xs_[i_nuclide].interp_factor; + auto f = p.neutron_xs_[i_nuclide].interp_factor; return (1.0 - f) * xs.value[i_grid-xs.threshold] + f * xs.value[i_grid-xs.threshold+1]; } //! Helper function to obtain neutron heating [eV] -double score_neutron_heating(const Particle* p, const Tally& tally, double flux, +double score_neutron_heating(const Particle& p, const Tally& tally, double flux, int rxn_bin, int i_nuclide, double atom_density) { double score; @@ -278,14 +278,14 @@ double score_neutron_heating(const Particle* p, const Tally& tally, double flux, const Nuclide& nuc {*data::nuclides[i_nuclide]}; heating_xs = get_nuclide_neutron_heating(p, nuc, rxn_bin, i_nuclide); if (tally.estimator_ == TallyEstimator::ANALOG) { - heating_xs /= p->neutron_xs_[i_nuclide].total; + heating_xs /= p.neutron_xs_[i_nuclide].total; } else { heating_xs *= atom_density; } } else { - if (p->material_ != MATERIAL_VOID) { + if (p.material_ != MATERIAL_VOID) { heating_xs = 0.0; - const Material& material {*model::materials[p->material_]}; + const Material& material {*model::materials[p.material_]}; for (auto i = 0; i< material.nuclide_.size(); ++i) { int j_nuclide = material.nuclide_[i]; double atom_density {material.atom_density_(i)}; @@ -293,7 +293,7 @@ double score_neutron_heating(const Particle* p, const Tally& tally, double flux, heating_xs += atom_density * get_nuclide_neutron_heating(p, nuc, rxn_bin, j_nuclide); } if (tally.estimator_ == TallyEstimator::ANALOG) { - heating_xs /= p->macro_xs_.total; + heating_xs /= p.macro_xs_.total; } } } @@ -304,9 +304,9 @@ double score_neutron_heating(const Particle* p, const Tally& tally, double flux, // reaction-wise heating cross section if (settings::survival_biasing) { // Account for the fact that some weight has been absorbed - score *= p->wgt_last_ + p->wgt_absorb_; + score *= p.wgt_last_ + p.wgt_absorb_; } else { - score *= p->wgt_last_; + score *= p.wgt_last_; } } return score; @@ -318,12 +318,12 @@ double score_neutron_heating(const Particle* p, const Tally& tally, double flux, //! neutrons produced with different energies. void -score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) +score_fission_eout(Particle& p, int i_tally, int i_score, int score_bin) { auto& tally {*model::tallies[i_tally]}; auto i_eout_filt = tally.filters()[tally.energyout_filter_]; - auto i_bin = p->filter_matches_[i_eout_filt].i_bin_; - auto bin_energyout = p->filter_matches_[i_eout_filt].bins_[i_bin]; + auto i_bin = p.filter_matches_[i_eout_filt].i_bin_; + auto bin_energyout = p.filter_matches_[i_eout_filt].bins_[i_bin]; const EnergyoutFilter& eo_filt {*dynamic_cast(model::tally_filters[i_eout_filt].get())}; @@ -334,8 +334,8 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) // rate. Otherwise, the sum of all nu-fission rates would be ~1.0. // loop over number of particles banked - for (auto i = 0; i < p->n_bank_; ++i) { - const auto& bank = p->nu_bank_[i]; + for (auto i = 0; i < p.n_bank_; ++i) { + const auto& bank = p.nu_bank_[i]; // get the delayed group auto g = bank.delayed_group; @@ -347,7 +347,7 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) // i_nuclide and atom_density arguments do not matter since this is an // analog estimator. if (tally.deriv_ != C_NONE) - apply_derivative_to_score(*p, i_tally, 0, 0., SCORE_NU_FISSION, score); + apply_derivative_to_score(p, i_tally, 0, 0., SCORE_NU_FISSION, score); if (!settings::run_CE && eo_filt.matches_transport_groups()) { @@ -359,7 +359,7 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) g_out = eo_filt.n_bins() - g_out - 1; // change outgoing energy bin - p->filter_matches_[i_eout_filt].bins_[i_bin] = g_out; + p.filter_matches_[i_eout_filt].bins_[i_bin] = g_out; } else { @@ -376,7 +376,7 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) } else { auto i_match = lower_bound_index(eo_filt.bins().begin(), eo_filt.bins().end(), E_out); - p->filter_matches_[i_eout_filt].bins_[i_bin] = i_match; + p.filter_matches_[i_eout_filt].bins_[i_bin] = i_match; } } @@ -390,7 +390,7 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) double filter_weight = 1.0; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); - auto& match {p->filter_matches_[i_filt]}; + auto& match {p.filter_matches_[i_filt]}; auto i_bin = match.i_bin_; filter_index += match.bins_[i_bin] * tally.strides(j); filter_weight *= match.weights_[i_bin]; @@ -419,13 +419,13 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) double filter_weight = 1.; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); - auto& match {p->filter_matches_[i_filt]}; + auto& match {p.filter_matches_[i_filt]}; auto i_bin = match.i_bin_; filter_weight *= match.weights_[i_bin]; } score_fission_delayed_dg(i_tally, d_bin, score*filter_weight, - i_score, p->filter_matches_); + i_score, p.filter_matches_); } } @@ -437,7 +437,7 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) double filter_weight = 1.; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); - auto& match {p->filter_matches_[i_filt]}; + auto& match {p.filter_matches_[i_filt]}; auto i_bin = match.i_bin_; filter_index += match.bins_[i_bin] * tally.strides(j); filter_weight *= match.weights_[i_bin]; @@ -451,7 +451,7 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) } // Reset outgoing energy bin and score index - p->filter_matches_[i_eout_filt].bins_[i_bin] = bin_energyout; + p.filter_matches_[i_eout_filt].bins_[i_bin] = bin_energyout; } //! Update tally results for continuous-energy tallies with any estimator. @@ -461,13 +461,13 @@ score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) //! is not used for analog tallies. void -score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, +score_general_ce(Particle& p, int i_tally, int start_index, int filter_index, double filter_weight, int i_nuclide, double atom_density, double flux) { Tally& tally {*model::tallies[i_tally]}; // Get the pre-collision energy of the particle. - auto E = p->E_last_; + auto E = p.E_last_; for (auto i = 0; i < tally.scores_.size(); ++i) { auto score_bin = tally.scores_[i]; @@ -483,14 +483,14 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->wgt_last_ + p->wgt_absorb_; + score = p.wgt_last_ + p.wgt_absorb_; } else { - score = p->wgt_last_; + score = p.wgt_last_; } - if (p->type_ == Particle::Type::neutron || - p->type_ == Particle::Type::photon) { - score *= flux / p->macro_xs_.total; + if (p.type_ == Particle::Type::neutron || + p.type_ == Particle::Type::photon) { + score *= flux / p.macro_xs_.total; } else { score = 0.; } @@ -507,20 +507,20 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = (p->wgt_last_ + p->wgt_absorb_) * flux; + score = (p.wgt_last_ + p.wgt_absorb_) * flux; } else { - score = p->wgt_last_ * flux; + score = p.wgt_last_ * flux; } } else { if (i_nuclide >= 0) { - if (p->type_ == Particle::Type::neutron) { - score = p->neutron_xs_[i_nuclide].total * atom_density * flux; - } else if (p->type_ == Particle::Type::photon) { - score = p->photon_xs_[i_nuclide].total * atom_density * flux; + if (p.type_ == Particle::Type::neutron) { + score = p.neutron_xs_[i_nuclide].total * atom_density * flux; + } else if (p.type_ == Particle::Type::photon) { + score = p.photon_xs_[i_nuclide].total * atom_density * flux; } } else { - score = p->macro_xs_.total * flux; + score = p.macro_xs_.total * flux; } } break; @@ -534,11 +534,11 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->wgt_last_ + p->wgt_absorb_; + score = p.wgt_last_ + p.wgt_absorb_; } else { - score = p->wgt_last_; + score = p.wgt_last_; } - score *= flux / p->macro_xs_.total; + score *= flux / p.macro_xs_.total; } else { score = flux; } @@ -550,17 +550,17 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, case SCORE_SCATTER: if (tally.estimator_ == TallyEstimator::ANALOG) { // Skip any event where the particle didn't scatter - if (p->event_ != TallyEvent::SCATTER) continue; + if (p.event_ != TallyEvent::SCATTER) continue; // Since only scattering events make it here, again we can use the // weight entering the collision as the estimator for the reaction rate - score = p->wgt_last_ * flux; + 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; + score = (p.neutron_xs_[i_nuclide].total + - p.neutron_xs_[i_nuclide].absorption) * atom_density * flux; } else { - score = (p->macro_xs_.total - - p->macro_xs_.absorption) * flux; + score = (p.macro_xs_.total + - p.macro_xs_.absorption) * flux; } } break; @@ -569,21 +569,21 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, case SCORE_NU_SCATTER: // Only analog estimators are available. // Skip any event where the particle didn't scatter - if (p->event_ != TallyEvent::SCATTER) continue; + if (p.event_ != TallyEvent::SCATTER) continue; // For scattering production, we need to use the pre-collision weight // times the yield as the estimate for the number of neutrons exiting a // reaction with neutrons in the exit channel - if (p->event_mt_ == ELASTIC || p->event_mt_ == N_LEVEL || - (p->event_mt_ >= N_N1 && p->event_mt_ <= N_NC)) { + if (p.event_mt_ == ELASTIC || p.event_mt_ == N_LEVEL || + (p.event_mt_ >= N_N1 && p.event_mt_ <= N_NC)) { // Don't waste time on very common reactions we know have // multiplicities of one. - score = p->wgt_last_ * flux; + score = p.wgt_last_ * flux; } else { // Get yield and apply to score auto m = - data::nuclides[p->event_nuclide_]->reaction_index_[p->event_mt_]; - const auto& rxn {*data::nuclides[p->event_nuclide_]->reactions_[m]}; - score = p->wgt_last_ * flux * (*rxn.products_[0].yield_)(E); + data::nuclides[p.event_nuclide_]->reaction_index_[p.event_mt_]; + const auto& rxn {*data::nuclides[p.event_nuclide_]->reactions_[m]}; + score = p.wgt_last_ * flux * (*rxn.products_[0].yield_)(E); } break; @@ -593,63 +593,63 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // No absorption events actually occur if survival biasing is on -- // just use weight absorbed in survival biasing - score = p->wgt_absorb_ * flux; + score = p.wgt_absorb_ * flux; } else { // Skip any event where the particle wasn't absorbed - if (p->event_ == TallyEvent::SCATTER) continue; + if (p.event_ == TallyEvent::SCATTER) continue; // All fission and absorption events will contribute here, so we // can just use the particle's weight entering the collision - score = p->wgt_last_ * flux; + score = p.wgt_last_ * flux; } } else { if (i_nuclide >= 0) { - score = p->neutron_xs_[i_nuclide].absorption * atom_density + score = p.neutron_xs_[i_nuclide].absorption * atom_density * flux; } else { - score = p->macro_xs_.absorption * flux; + score = p.macro_xs_.absorption * flux; } } break; case SCORE_FISSION: - if (p->macro_xs_.absorption == 0) continue; + if (p.macro_xs_.absorption == 0) continue; if (tally.estimator_ == TallyEstimator::ANALOG) { if (settings::survival_biasing) { // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission - if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { - score = p->wgt_absorb_ - * p->neutron_xs_[p->event_nuclide_].fission - / p->neutron_xs_[p->event_nuclide_].absorption * flux; + if (p.neutron_xs_[p.event_nuclide_].absorption > 0) { + score = p.wgt_absorb_ + * p.neutron_xs_[p.event_nuclide_].fission + / p.neutron_xs_[p.event_nuclide_].absorption * flux; } else { score = 0.; } } else { // Skip any non-absorption events - if (p->event_ == TallyEvent::SCATTER) continue; + if (p.event_ == TallyEvent::SCATTER) continue; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - score = p->wgt_last_ - * p->neutron_xs_[p->event_nuclide_].fission - / p->neutron_xs_[p->event_nuclide_].absorption * flux; + score = p.wgt_last_ + * p.neutron_xs_[p.event_nuclide_].fission + / p.neutron_xs_[p.event_nuclide_].absorption * flux; } } else { if (i_nuclide >= 0) { - score = p->neutron_xs_[i_nuclide].fission * atom_density * flux; + score = p.neutron_xs_[i_nuclide].fission * atom_density * flux; } else { - score = p->macro_xs_.fission * flux; + score = p.macro_xs_.fission * flux; } } break; case SCORE_NU_FISSION: - if (p->macro_xs_.absorption == 0) continue; + if (p.macro_xs_.absorption == 0) continue; if (tally.estimator_ == TallyEstimator::ANALOG) { - if (settings::survival_biasing || p->fission_) { + if (settings::survival_biasing || p.fission_) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -661,38 +661,38 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // nu-fission - if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { - score = p->wgt_absorb_ - * p->neutron_xs_[p->event_nuclide_].nu_fission - / p->neutron_xs_[p->event_nuclide_].absorption * flux; + if (p.neutron_xs_[p.event_nuclide_].absorption > 0) { + score = p.wgt_absorb_ + * p.neutron_xs_[p.event_nuclide_].nu_fission + / p.neutron_xs_[p.event_nuclide_].absorption * flux; } else { score = 0.; } } else { // Skip any non-fission events - if (!p->fission_) continue; + if (!p.fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since // this was weighted by 1/keff, we multiply by keff to get the proper // score. - score = simulation::keff * p->wgt_bank_ * flux; + score = simulation::keff * p.wgt_bank_ * flux; } } else { if (i_nuclide >= 0) { - score = p->neutron_xs_[i_nuclide].nu_fission * atom_density + score = p.neutron_xs_[i_nuclide].nu_fission * atom_density * flux; } else { - score = p->macro_xs_.nu_fission * flux; + score = p.macro_xs_.nu_fission * flux; } } break; case SCORE_PROMPT_NU_FISSION: - if (p->macro_xs_.absorption == 0) continue; + if (p.macro_xs_.absorption == 0) continue; if (tally.estimator_ == TallyEstimator::ANALOG) { - if (settings::survival_biasing || p->fission_) { + if (settings::survival_biasing || p.fission_) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -704,43 +704,43 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // prompt-nu-fission - if (p->neutron_xs_[p->event_nuclide_].absorption > 0) { - score = p->wgt_absorb_ - * p->neutron_xs_[p->event_nuclide_].fission - * data::nuclides[p->event_nuclide_] + if (p.neutron_xs_[p.event_nuclide_].absorption > 0) { + score = p.wgt_absorb_ + * p.neutron_xs_[p.event_nuclide_].fission + * data::nuclides[p.event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::prompt) - / p->neutron_xs_[p->event_nuclide_].absorption * flux; + / p.neutron_xs_[p.event_nuclide_].absorption * flux; } else { score = 0.; } } else { // Skip any non-fission events - if (!p->fission_) continue; + if (!p.fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since // this was weighted by 1/keff, we multiply by keff to get the proper // score. - auto n_delayed = std::accumulate(p->n_delayed_bank_, - p->n_delayed_bank_+MAX_DELAYED_GROUPS, 0); - auto prompt_frac = 1. - n_delayed / static_cast(p->n_bank_); - score = simulation::keff * p->wgt_bank_ * prompt_frac * flux; + auto n_delayed = std::accumulate(p.n_delayed_bank_, + p.n_delayed_bank_+MAX_DELAYED_GROUPS, 0); + auto prompt_frac = 1. - n_delayed / static_cast(p.n_bank_); + score = simulation::keff * p.wgt_bank_ * prompt_frac * flux; } } else { if (i_nuclide >= 0) { - score = p->neutron_xs_[i_nuclide].fission + score = p.neutron_xs_[i_nuclide].fission * data::nuclides[i_nuclide] ->nu(E, ReactionProduct::EmissionMode::prompt) * atom_density * flux; } else { score = 0.; // Add up contributions from each nuclide in the material. - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; + if (p.material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p.material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - score += p->neutron_xs_[j_nuclide].fission + score += p.neutron_xs_[j_nuclide].fission * data::nuclides[j_nuclide] ->nu(E, ReactionProduct::EmissionMode::prompt) * atom_density * flux; @@ -752,9 +752,9 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, case SCORE_DELAYED_NU_FISSION: - if (p->macro_xs_.absorption == 0) continue; + if (p.macro_xs_.absorption == 0) continue; if (tally.estimator_ == TallyEstimator::ANALOG) { - if (settings::survival_biasing || p->fission_) { + if (settings::survival_biasing || p.fission_) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -766,8 +766,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - if (p->neutron_xs_[p->event_nuclide_].absorption > 0 - && data::nuclides[p->event_nuclide_]->fissionable_) { + if (p.neutron_xs_[p.event_nuclide_].absorption > 0 + && data::nuclides[p.event_nuclide_]->fissionable_) { if (tally.delayedgroup_filter_ != C_NONE) { auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; const DelayedGroupFilter& filt @@ -776,29 +776,29 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin]; - auto yield = data::nuclides[p->event_nuclide_] + auto yield = data::nuclides[p.event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::delayed, d); - score = p->wgt_absorb_ * yield - * p->neutron_xs_[p->event_nuclide_].fission - / p->neutron_xs_[p->event_nuclide_].absorption * flux; + score = p.wgt_absorb_ * yield + * p.neutron_xs_[p.event_nuclide_].fission + / p.neutron_xs_[p.event_nuclide_].absorption * flux; score_fission_delayed_dg(i_tally, d_bin, score, - score_index, p->filter_matches_); + score_index, p.filter_matches_); } continue; } else { // If the delayed group filter is not present, compute the score // by multiplying the absorbed weight by the fraction of the // delayed-nu-fission xs to the absorption xs - score = p->wgt_absorb_ - * p->neutron_xs_[p->event_nuclide_].fission - * data::nuclides[p->event_nuclide_] + score = p.wgt_absorb_ + * p.neutron_xs_[p.event_nuclide_].fission + * data::nuclides[p.event_nuclide_] ->nu(E, ReactionProduct::EmissionMode::delayed) - / p->neutron_xs_[p->event_nuclide_].absorption *flux; + / p.neutron_xs_[p.event_nuclide_].absorption *flux; } } } else { // Skip any non-fission events - if (!p->fission_) continue; + if (!p.fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since @@ -814,16 +814,16 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin]; - score = simulation::keff * p->wgt_bank_ / p->n_bank_ - * p->n_delayed_bank_[d-1] * flux; - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); + score = simulation::keff * p.wgt_bank_ / p.n_bank_ + * p.n_delayed_bank_[d-1] * flux; + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); } continue; } else { // Add the contribution from all delayed groups - auto n_delayed = std::accumulate(p->n_delayed_bank_, - p->n_delayed_bank_+MAX_DELAYED_GROUPS, 0); - score = simulation::keff * p->wgt_bank_ / p->n_bank_ * n_delayed + auto n_delayed = std::accumulate(p.n_delayed_bank_, + p.n_delayed_bank_+MAX_DELAYED_GROUPS, 0); + score = simulation::keff * p.wgt_bank_ / p.n_bank_ * n_delayed * flux; } } @@ -839,15 +839,15 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, auto d = filt.groups()[d_bin]; auto yield = data::nuclides[i_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed, d); - score = p->neutron_xs_[i_nuclide].fission * yield + score = p.neutron_xs_[i_nuclide].fission * yield * atom_density * flux; - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); } continue; } else { // If the delayed group filter is not present, compute the score // by multiplying the delayed-nu-fission macro xs by the flux - score = p->neutron_xs_[i_nuclide].fission + score = p.neutron_xs_[i_nuclide].fission * data::nuclides[i_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed) * atom_density * flux; @@ -858,8 +858,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, const DelayedGroupFilter& filt {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; + if (p.material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p.material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -868,22 +868,22 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, auto d = filt.groups()[d_bin]; auto yield = data::nuclides[j_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed, d); - score = p->neutron_xs_[j_nuclide].fission * yield + score = p.neutron_xs_[j_nuclide].fission * yield * atom_density * flux; score_fission_delayed_dg(i_tally, d_bin, score, - score_index, p->filter_matches_); + score_index, p.filter_matches_); } } } continue; } else { score = 0.; - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; + if (p.material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p.material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - score += p->neutron_xs_[j_nuclide].fission + score += p.neutron_xs_[j_nuclide].fission * data::nuclides[j_nuclide] ->nu(E, ReactionProduct::EmissionMode::delayed) * atom_density * flux; @@ -896,14 +896,14 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, case SCORE_DECAY_RATE: - if (p->macro_xs_.absorption == 0) continue; + if (p.macro_xs_.absorption == 0) continue; if (tally.estimator_ == TallyEstimator::ANALOG) { if (settings::survival_biasing) { // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // delayed-nu-fission - const auto& nuc {*data::nuclides[p->event_nuclide_]}; - if (p->neutron_xs_[p->event_nuclide_].absorption > 0 + const auto& nuc {*data::nuclides[p.event_nuclide_]}; + if (p.neutron_xs_[p.event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; if (tally.delayedgroup_filter_ != C_NONE) { @@ -917,12 +917,12 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; - score = p->wgt_absorb_ * yield - * p->neutron_xs_[p->event_nuclide_].fission - / p->neutron_xs_[p->event_nuclide_].absorption + score = p.wgt_absorb_ * yield + * p.neutron_xs_[p.event_nuclide_].fission + / p.neutron_xs_[p.event_nuclide_].absorption * rate * flux; score_fission_delayed_dg(i_tally, d_bin, score, - score_index, p->filter_matches_); + score_index, p.filter_matches_); } continue; } else { @@ -940,15 +940,15 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; - score += rate * p->wgt_absorb_ - * p->neutron_xs_[p->event_nuclide_].fission * yield - / p->neutron_xs_[p->event_nuclide_].absorption * flux; + score += rate * p.wgt_absorb_ + * p.neutron_xs_[p.event_nuclide_].fission * yield + / p.neutron_xs_[p.event_nuclide_].absorption * flux; } } } } else { // Skip any non-fission events - if (!p->fission_) continue; + if (!p.fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since @@ -957,11 +957,11 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, // ones are delayed. If a delayed neutron is encountered, add its // contribution to the fission bank to the score. score = 0.; - for (auto i = 0; i < p->n_bank_; ++i) { - const auto& bank = p->nu_bank_[i]; + for (auto i = 0; i < p.n_bank_; ++i) { + const auto& bank = p.nu_bank_[i]; auto g = bank.delayed_group; if (g != 0) { - const auto& nuc {*data::nuclides[p->event_nuclide_]}; + const auto& nuc {*data::nuclides[p.event_nuclide_]}; const auto& rxn {*nuc.fission_rx_[0]}; auto rate = rxn.products_[g].decay_rate_; score += simulation::keff * bank.wgt * rate * flux; @@ -975,7 +975,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, auto d = filt.groups()[d_bin]; if (d == g) score_fission_delayed_dg(i_tally, d_bin, score, - score_index, p->filter_matches_); + score_index, p.filter_matches_); } score = 0.; } @@ -998,9 +998,9 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; - score = p->neutron_xs_[i_nuclide].fission * yield * flux + score = p.neutron_xs_[i_nuclide].fission * yield * flux * atom_density * rate; - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); } continue; } else { @@ -1014,7 +1014,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; - score += p->neutron_xs_[i_nuclide].fission * flux + score += p.neutron_xs_[i_nuclide].fission * flux * yield * atom_density * rate; } } @@ -1024,8 +1024,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, const DelayedGroupFilter& filt {*dynamic_cast( model::tally_filters[i_dg_filt].get())}; - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; + if (p.material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p.material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1038,10 +1038,10 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d); auto rate = rxn.products_[d].decay_rate_; - score = p->neutron_xs_[j_nuclide].fission * yield + score = p.neutron_xs_[j_nuclide].fission * yield * flux * atom_density * rate; score_fission_delayed_dg(i_tally, d_bin, score, - score_index, p->filter_matches_); + score_index, p.filter_matches_); } } } @@ -1049,8 +1049,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, continue; } else { score = 0.; - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; + if (p.material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p.material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1066,7 +1066,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, auto yield = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1); auto rate = rxn.products_[d+1].decay_rate_; - score += p->neutron_xs_[j_nuclide].fission + score += p.neutron_xs_[j_nuclide].fission * yield * atom_density * flux * rate; } } @@ -1079,7 +1079,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, case SCORE_KAPPA_FISSION: - if (p->macro_xs_.absorption == 0.) continue; + if (p.macro_xs_.absorption == 0.) continue; score = 0.; // Kappa-fission values are determined from the Q-value listed for the // fission cross section. @@ -1088,27 +1088,27 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value - const auto& nuc {*data::nuclides[p->event_nuclide_]}; - if (p->neutron_xs_[p->event_nuclide_].absorption > 0 + const auto& nuc {*data::nuclides[p.event_nuclide_]}; + if (p.neutron_xs_[p.event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score = p->wgt_absorb_ * rxn.q_value_ - * p->neutron_xs_[p->event_nuclide_].fission - / p->neutron_xs_[p->event_nuclide_].absorption * flux; + score = p.wgt_absorb_ * rxn.q_value_ + * p.neutron_xs_[p.event_nuclide_].fission + / p.neutron_xs_[p.event_nuclide_].absorption * flux; } } else { // Skip any non-absorption events - if (p->event_ == TallyEvent::SCATTER) continue; + if (p.event_ == TallyEvent::SCATTER) continue; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - const auto& nuc {*data::nuclides[p->event_nuclide_]}; - if (p->neutron_xs_[p->event_nuclide_].absorption > 0 + const auto& nuc {*data::nuclides[p.event_nuclide_]}; + if (p.neutron_xs_[p.event_nuclide_].absorption > 0 && nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score = p->wgt_last_ * rxn.q_value_ - * p->neutron_xs_[p->event_nuclide_].fission - / p->neutron_xs_[p->event_nuclide_].absorption * flux; + score = p.wgt_last_ * rxn.q_value_ + * p.neutron_xs_[p.event_nuclide_].fission + / p.neutron_xs_[p.event_nuclide_].absorption * flux; } } } else { @@ -1116,19 +1116,19 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, const auto& nuc {*data::nuclides[i_nuclide]}; if (nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score = rxn.q_value_ * p->neutron_xs_[i_nuclide].fission + score = rxn.q_value_ * p.neutron_xs_[i_nuclide].fission * atom_density * flux; } } else { - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; + if (p.material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p.material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); const auto& nuc {*data::nuclides[j_nuclide]}; if (nuc.fissionable_) { const auto& rxn {*nuc.fission_rx_[0]}; - score += rxn.q_value_ * p->neutron_xs_[j_nuclide].fission + score += rxn.q_value_ * p.neutron_xs_[j_nuclide].fission * atom_density * flux; } } @@ -1148,23 +1148,23 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, case ELASTIC: if (tally.estimator_ == TallyEstimator::ANALOG) { // Check if event MT matches - if (p->event_mt_ != ELASTIC) continue; - score = p->wgt_last_ * flux; + if (p.event_mt_ != ELASTIC) continue; + score = p.wgt_last_ * flux; } else { if (i_nuclide >= 0) { - if (p->neutron_xs_[i_nuclide].elastic == CACHE_INVALID) - data::nuclides[i_nuclide]->calculate_elastic_xs(*p); - score = p->neutron_xs_[i_nuclide].elastic * atom_density * flux; + if (p.neutron_xs_[i_nuclide].elastic == CACHE_INVALID) + data::nuclides[i_nuclide]->calculate_elastic_xs(p); + score = p.neutron_xs_[i_nuclide].elastic * atom_density * flux; } else { score = 0.; - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; + if (p.material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p.material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - if (p->neutron_xs_[j_nuclide].elastic == CACHE_INVALID) - data::nuclides[j_nuclide]->calculate_elastic_xs(*p); - score += p->neutron_xs_[j_nuclide].elastic * atom_density + if (p.neutron_xs_[j_nuclide].elastic == CACHE_INVALID) + data::nuclides[j_nuclide]->calculate_elastic_xs(p); + score += p.neutron_xs_[j_nuclide].elastic * atom_density * flux; } } @@ -1174,7 +1174,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, case SCORE_FISS_Q_PROMPT: case SCORE_FISS_Q_RECOV: - if (p->macro_xs_.absorption == 0.) continue; + if (p.macro_xs_.absorption == 0.) continue; score = score_fission_q(p, score_bin, tally, flux, i_nuclide, atom_density); break; @@ -1187,8 +1187,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, case N_A: if (tally.estimator_ == TallyEstimator::ANALOG) { // Check if the event MT matches - if (p->event_mt_ != score_bin) continue; - score = p->wgt_last_ * flux; + if (p.event_mt_ != score_bin) continue; + score = p.wgt_last_ * flux; } else { int m; switch (score_bin) { @@ -1200,16 +1200,16 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, case N_4N: m = 5; break; } if (i_nuclide >= 0) { - score = p->neutron_xs_[i_nuclide].reaction[m] * atom_density + score = p.neutron_xs_[i_nuclide].reaction[m] * atom_density * flux; } else { score = 0.; - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; + if (p.material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p.material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); - score += p->neutron_xs_[j_nuclide].reaction[m] + score += p.neutron_xs_[j_nuclide].reaction[m] * atom_density * flux; } } @@ -1219,22 +1219,22 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, case HEATING: - if (p->type_ == Particle::Type::neutron) { + if (p.type_ == Particle::Type::neutron) { score = score_neutron_heating(p, tally, flux, HEATING, i_nuclide, atom_density); } else { // The energy deposited is the difference between the pre-collision and // post-collision energy... - score = E - p->E_; + score = E - p.E_; // ...less the energy of any secondary particles since they will be // transported individually later - const auto& bank = p->secondary_bank_; - for (auto it = bank.end() - p->n_bank_second_; it < bank.end(); ++it) { + const auto& bank = p.secondary_bank_; + for (auto it = bank.end() - p.n_bank_second_; it < bank.end(); ++it) { score -= it->E; } - score *= p->wgt_last_; + score *= p.wgt_last_; } break; @@ -1242,14 +1242,14 @@ 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_ != Particle::Type::neutron) continue; if (tally.estimator_ == TallyEstimator::ANALOG) { // Any other score is assumed to be a MT number. Thus, we just need // to check if it matches the MT number of the event - if (p->event_mt_ != score_bin) continue; - score = p->wgt_last_*flux; + if (p.event_mt_ != score_bin) continue; + score = p.wgt_last_*flux; } else { // Any other cross section has to be calculated on-the-fly if (score_bin < 2) fatal_error("Invalid score type on tally " @@ -1260,10 +1260,10 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, auto m = nuc.reaction_index_[score_bin]; if (m == C_NONE) continue; const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = p->neutron_xs_[i_nuclide].index_temp; + auto i_temp = p.neutron_xs_[i_nuclide].index_temp; if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = p->neutron_xs_[i_nuclide].index_grid; - auto f = p->neutron_xs_[i_nuclide].interp_factor; + auto i_grid = p.neutron_xs_[i_nuclide].index_grid; + auto f = p.neutron_xs_[i_nuclide].interp_factor; const auto& xs {rxn.xs_[i_temp]}; if (i_grid >= xs.threshold) { score = ((1.0 - f) * xs.value[i_grid-xs.threshold] @@ -1271,8 +1271,8 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, } } } else { - if (p->material_ != MATERIAL_VOID) { - const Material& material {*model::materials[p->material_]}; + if (p.material_ != MATERIAL_VOID) { + const Material& material {*model::materials[p.material_]}; for (auto i = 0; i < material.nuclide_.size(); ++i) { auto j_nuclide = material.nuclide_[i]; auto atom_density = material.atom_density_(i); @@ -1280,10 +1280,10 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, auto m = nuc.reaction_index_[score_bin]; if (m == C_NONE) continue; const auto& rxn {*nuc.reactions_[m]}; - auto i_temp = p->neutron_xs_[j_nuclide].index_temp; + auto i_temp = p.neutron_xs_[j_nuclide].index_temp; if (i_temp >= 0) { // Can be false due to multipole - auto i_grid = p->neutron_xs_[j_nuclide].index_grid; - auto f = p->neutron_xs_[j_nuclide].interp_factor; + auto i_grid = p.neutron_xs_[j_nuclide].index_grid; + auto f = p.neutron_xs_[j_nuclide].interp_factor; const auto& xs {rxn.xs_[i_temp]}; if (i_grid >= xs.threshold) { score += ((1.0 - f) * xs.value[i_grid-xs.threshold] @@ -1299,7 +1299,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, // Add derivative information on score for differential tallies. if (tally.deriv_ != C_NONE) - apply_derivative_to_score(*p, i_tally, i_nuclide, atom_density, score_bin, + apply_derivative_to_score(p, i_tally, i_nuclide, atom_density, score_bin, score); // Update tally results @@ -1314,7 +1314,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index, //! argument is really just used for filter weights. void -score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, +score_general_mg(Particle& p, int i_tally, int start_index, int filter_index, double filter_weight, int i_nuclide, double atom_density, double flux) { auto& tally {*model::tallies[i_tally]}; @@ -1329,40 +1329,40 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, // Then we either are alive and had a scatter (and so g changed), // or are dead and g did not change - if (p->alive_) { - p_u = p->u_last_; - p_g = p->g_last_; + if (p.alive_) { + p_u = p.u_last_; + p_g = p.g_last_; } else { - p_u = p->u_local(); - p_g = p->g_; + p_u = p.u_local(); + p_g = p.g_; } - } else if (p->event_ == TallyEvent::SCATTER) { + } else if (p.event_ == TallyEvent::SCATTER) { // Then the energy group has been changed by the scattering routine // meaning gin is now in p % last_g - p_u = p->u_last_; - p_g = p->g_last_; + p_u = p.u_last_; + p_g = p.g_last_; } else { // No scatter, no change in g. - p_u = p->u_local(); - p_g = p->g_; + p_u = p.u_local(); + p_g = p.g_; } } else { // No actual collision so g has not changed. - p_u = p->u_local(); - p_g = p->g_; + p_u = p.u_local(); + p_g = p.g_; } // For shorthand, assign pointers to the material and nuclide xs set auto& nuc_xs = data::mg.nuclides_[i_nuclide]; - auto& macro_xs = data::mg.macro_xs_[p->material_]; + auto& macro_xs = data::mg.macro_xs_[p.material_]; // Find the temperature and angle indices of interest macro_xs.set_angle_index(p_u); if (i_nuclide >= 0) { - nuc_xs.set_temperature_index(p->sqrtkT_); + nuc_xs.set_temperature_index(p.sqrtkT_); nuc_xs.set_angle_index(p_u); } @@ -1383,11 +1383,11 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->wgt_last_ + p->wgt_absorb_; + score = p.wgt_last_ + p.wgt_absorb_; } else { - score = p->wgt_last_; + score = p.wgt_last_; } - score *= flux / p->macro_xs_.total; + score *= flux / p.macro_xs_.total; } else { score = flux; } @@ -1401,9 +1401,9 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->wgt_last_ + p->wgt_absorb_; + score = p.wgt_last_ + p.wgt_absorb_; } else { - score = p->wgt_last_; + score = p.wgt_last_; } //TODO: should flux be multiplied in above instead of below? if (i_nuclide >= 0) { @@ -1414,7 +1414,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, if (i_nuclide >= 0) { score = atom_density * flux * nuc_xs.get_xs(MgxsType::TOTAL, p_g); } else { - score = p->macro_xs_.total * flux; + score = p.macro_xs_.total * flux; } } break; @@ -1429,9 +1429,9 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // We need to account for the fact that some weight was already // absorbed - score = p->wgt_last_ + p->wgt_absorb_; + score = p.wgt_last_ + p.wgt_absorb_; } else { - score = p->wgt_last_; + score = p.wgt_last_; } if (i_nuclide >= 0) { score *= flux * nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g) @@ -1453,23 +1453,23 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, case SCORE_SCATTER: if (tally.estimator_ == TallyEstimator::ANALOG) { // Skip any event where the particle didn't scatter - if (p->event_ != TallyEvent::SCATTER) continue; + if (p.event_ != TallyEvent::SCATTER) continue; // Since only scattering events make it here, again we can use the // weight entering the collision as the estimator for the reaction rate - score = p->wgt_last_ * flux; + score = p.wgt_last_ * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs( - MgxsType::SCATTER_FMU_MULT, p->g_last_, &p->g_, &p->mu_, nullptr) + MgxsType::SCATTER_FMU_MULT, p.g_last_, &p.g_, &p.mu_, nullptr) / macro_xs.get_xs( - MgxsType::SCATTER_FMU_MULT, p->g_last_, &p->g_, &p->mu_, nullptr); + MgxsType::SCATTER_FMU_MULT, p.g_last_, &p.g_, &p.mu_, nullptr); } } else { if (i_nuclide >= 0) { score = atom_density * flux * nuc_xs.get_xs( - MgxsType::SCATTER_MULT, p_g, nullptr, &p->mu_, nullptr); + MgxsType::SCATTER_MULT, p_g, nullptr, &p.mu_, nullptr); } else { score = flux * macro_xs.get_xs( - MgxsType::SCATTER_MULT, p_g, nullptr, &p->mu_, nullptr); + MgxsType::SCATTER_MULT, p_g, nullptr, &p.mu_, nullptr); } } break; @@ -1478,20 +1478,20 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, case SCORE_NU_SCATTER: if (tally.estimator_ == TallyEstimator::ANALOG) { // Skip any event where the particle didn't scatter - if (p->event_ != TallyEvent::SCATTER) continue; + if (p.event_ != TallyEvent::SCATTER) continue; // For scattering production, we need to use the pre-collision weight // times the multiplicity as the estimate for the number of neutrons // exiting a reaction with neutrons in the exit channel - score = p->wgt_ * flux; + score = p.wgt_ * flux; // Since we transport based on material data, the angle selected // was not selected from the f(mu) for the nuclide. Therefore // adjust the score by the actual probability for that nuclide. if (i_nuclide >= 0) { score *= atom_density - * nuc_xs.get_xs(MgxsType::SCATTER_FMU, p->g_last_, &p->g_, - &p->mu_, nullptr) - / macro_xs.get_xs(MgxsType::SCATTER_FMU, p->g_last_, &p->g_, - &p->mu_, nullptr); + * nuc_xs.get_xs(MgxsType::SCATTER_FMU, p.g_last_, &p.g_, + &p.mu_, nullptr) + / macro_xs.get_xs(MgxsType::SCATTER_FMU, p.g_last_, &p.g_, + &p.mu_, nullptr); } } else { if (i_nuclide >= 0) { @@ -1508,13 +1508,13 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, if (settings::survival_biasing) { // No absorption events actually occur if survival biasing is on -- // just use weight absorbed in survival biasing - score = p->wgt_absorb_ * flux; + score = p.wgt_absorb_ * flux; } else { // Skip any event where the particle wasn't absorbed - if (p->event_ == TallyEvent::SCATTER) continue; + if (p.event_ == TallyEvent::SCATTER) continue; // All fission and absorption events will contribute here, so we // can just use the particle's weight entering the collision - score = p->wgt_last_ * flux; + score = p.wgt_last_ * flux; } if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::ABSORPTION, p_g) @@ -1525,7 +1525,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, score = atom_density * flux * nuc_xs.get_xs(MgxsType::ABSORPTION, p_g); } else { - score = p->macro_xs_.absorption * flux; + score = p.macro_xs_.absorption * flux; } } break; @@ -1537,14 +1537,14 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission - score = p->wgt_absorb_ * flux; + score = p.wgt_absorb_ * flux; } else { // Skip any non-absorption events - if (p->event_ == TallyEvent::SCATTER) continue; + if (p.event_ == TallyEvent::SCATTER) continue; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - score = p->wgt_last_ * flux; + score = p.wgt_last_ * flux; } if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) @@ -1565,7 +1565,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, case SCORE_NU_FISSION: if (tally.estimator_ == TallyEstimator::ANALOG) { - if (settings::survival_biasing || p->fission_) { + if (settings::survival_biasing || p.fission_) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -1577,7 +1577,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // nu-fission - score = p->wgt_absorb_ * flux; + score = p.wgt_absorb_ * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::NU_FISSION, p_g) / macro_xs.get_xs(MgxsType::ABSORPTION, p_g); @@ -1587,13 +1587,13 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, } } else { // Skip any non-fission events - if (!p->fission_) continue; + if (!p.fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since // this was weighted by 1/keff, we multiply by keff to get the proper // score. - score = simulation::keff * p->wgt_bank_ * flux; + score = simulation::keff * p.wgt_bank_ * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) / macro_xs.get_xs(MgxsType::FISSION, p_g); @@ -1612,7 +1612,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, case SCORE_PROMPT_NU_FISSION: if (tally.estimator_ == TallyEstimator::ANALOG) { - if (settings::survival_biasing || p->fission_) { + if (settings::survival_biasing || p.fission_) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -1624,7 +1624,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // prompt-nu-fission - score = p->wgt_absorb_ * flux; + score = p.wgt_absorb_ * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g) @@ -1635,16 +1635,16 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, } } else { // Skip any non-fission events - if (!p->fission_) continue; + if (!p.fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since // this was weighted by 1/keff, we multiply by keff to get the proper // score. - auto n_delayed = std::accumulate(p->n_delayed_bank_, - p->n_delayed_bank_+MAX_DELAYED_GROUPS, 0); - auto prompt_frac = 1. - n_delayed / static_cast(p->n_bank_); - score = simulation::keff * p->wgt_bank_ * prompt_frac * flux; + auto n_delayed = std::accumulate(p.n_delayed_bank_, + p.n_delayed_bank_+MAX_DELAYED_GROUPS, 0); + auto prompt_frac = 1. - n_delayed / static_cast(p.n_bank_); + score = simulation::keff * p.wgt_bank_ * prompt_frac * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) / macro_xs.get_xs(MgxsType::FISSION, p_g); @@ -1663,7 +1663,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, case SCORE_DELAYED_NU_FISSION: if (tally.estimator_ == TallyEstimator::ANALOG) { - if (settings::survival_biasing || p->fission_) { + if (settings::survival_biasing || p.fission_) { if (tally.energyout_filter_ != C_NONE) { // Fission has multiple outgoing neutrons so this helper function // is used to handle scoring the multiple filter bins. @@ -1685,7 +1685,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin] - 1; - score = p->wgt_absorb_ * flux; + score = p.wgt_absorb_ * flux; if (i_nuclide >= 0) { score *= nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) @@ -1695,14 +1695,14 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, nullptr, nullptr, &d) / abs_xs; } - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); } continue; } else { // If the delayed group filter is not present, compute the score // by multiplying the absorbed weight by the fraction of the // delayed-nu-fission xs to the absorption xs - score = p->wgt_absorb_ * flux; + score = p.wgt_absorb_ * flux; if (i_nuclide >= 0) { score *= nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g) / abs_xs; @@ -1714,7 +1714,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, } } else { // Skip any non-fission events - if (!p->fission_) continue; + if (!p.fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since @@ -1730,20 +1730,20 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin]; - score = simulation::keff * p->wgt_bank_ / p->n_bank_ - * p->n_delayed_bank_[d-1] * flux; + score = simulation::keff * p.wgt_bank_ / p.n_bank_ + * p.n_delayed_bank_[d-1] * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) / macro_xs.get_xs(MgxsType::FISSION, p_g); } - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); } continue; } else { // Add the contribution from all delayed groups - auto n_delayed = std::accumulate(p->n_delayed_bank_, - p->n_delayed_bank_+MAX_DELAYED_GROUPS, 0); - score = simulation::keff * p->wgt_bank_ / p->n_bank_ * n_delayed + auto n_delayed = std::accumulate(p.n_delayed_bank_, + p.n_delayed_bank_+MAX_DELAYED_GROUPS, 0); + score = simulation::keff * p.wgt_bank_ / p.n_bank_ * n_delayed * flux; if (i_nuclide >= 0) { score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g) @@ -1767,7 +1767,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, score = flux * macro_xs.get_xs( MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); } continue; } else { @@ -1798,7 +1798,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, // Tally each delayed group bin individually for (auto d_bin = 0; d_bin < filt.n_bins(); ++d_bin) { auto d = filt.groups()[d_bin] - 1; - score = p->wgt_absorb_ * flux; + score = p.wgt_absorb_ * flux; if (i_nuclide >= 0) { score *= nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) @@ -1810,7 +1810,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) / abs_xs; } - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); } continue; } else { @@ -1821,13 +1821,13 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, score = 0.; for (auto d = 0; d < data::mg.num_delayed_groups_; ++d) { if (i_nuclide >= 0) { - score += p->wgt_absorb_ * flux + score += p.wgt_absorb_ * flux * nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) / abs_xs; } else { - score += p->wgt_absorb_ * flux + score += p.wgt_absorb_ * flux * macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d) * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, @@ -1838,7 +1838,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, } } else { // Skip any non-fission events - if (!p->fission_) continue; + if (!p.fission_) continue; // If there is no outgoing energy filter, than we only need to score // to one bin. For the score to be 'analog', we need to score the // number of particles that were banked in the fission bank. Since @@ -1847,8 +1847,8 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, // ones are delayed. If a delayed neutron is encountered, add its // contribution to the fission bank to the score. score = 0.; - for (auto i = 0; i < p->n_bank_; ++i) { - const auto& bank = p->nu_bank_[i]; + for (auto i = 0; i < p.n_bank_; ++i) { + const auto& bank = p.nu_bank_[i]; auto d = bank.delayed_group - 1; if (d != -1) { if (i_nuclide >= 0) { @@ -1870,7 +1870,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, auto dg = filt.groups()[d_bin]; if (dg == d + 1) score_fission_delayed_dg(i_tally, d_bin, score, - score_index, p->filter_matches_); + score_index, p.filter_matches_); } score = 0.; } @@ -1900,7 +1900,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } - score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p.filter_matches_); } continue; } else { @@ -1927,14 +1927,14 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, // No fission events occur if survival biasing is on -- need to // calculate fraction of absorptions that would have resulted in // fission scaled by the Q-value - score = p->wgt_absorb_ * flux; + score = p.wgt_absorb_ * flux; } else { // Skip any non-absorption events - if (p->event_ == TallyEvent::SCATTER) continue; + if (p.event_ == TallyEvent::SCATTER) continue; // All fission events will contribute, so again we can use particle's // weight entering the collision as the estimate for the fission // reaction rate - score = p->wgt_last_ * flux; + score = p.wgt_last_ * flux; } if (i_nuclide >= 0) { score *= atom_density @@ -1974,11 +1974,11 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, //! Tally rates for when the user requests a tally on all nuclides. void -score_all_nuclides(Particle* p, int i_tally, double flux, +score_all_nuclides(Particle& p, int i_tally, double flux, int filter_index, double filter_weight) { const Tally& tally {*model::tallies[i_tally]}; - const Material& material {*model::materials[p->material_]}; + const Material& material {*model::materials[p.material_]}; // Score all individual nuclide reaction rates. for (auto i = 0; i < material.nuclide_.size(); ++i) { @@ -2009,13 +2009,13 @@ score_all_nuclides(Particle* p, int i_tally, double flux, } } -void score_analog_tally_ce(Particle* p) +void score_analog_tally_ce(Particle& p) { // Since electrons/positrons are not transported, we assign a flux of zero. // Note that the heating score does NOT use the flux and will be non-zero for // electrons/positrons. double flux = - (p->type_ == Particle::Type::neutron || p->type_ == Particle::Type::photon) ? + (p.type_ == Particle::Type::neutron || p.type_ == Particle::Type::photon) ? 1.0 : 0.0; for (auto i_tally : model::active_analog_tallies) { @@ -2025,7 +2025,7 @@ void score_analog_tally_ce(Particle* p) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true, &p->filter_matches_); + auto end = FilterBinIter(tally, true, &p.filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2041,7 +2041,7 @@ void score_analog_tally_ce(Particle* p) // Tally this event in the present nuclide bin if that bin represents // the event nuclide or the total material. Note that the atomic // density argument for score_general is not used for analog tallies. - if (i_nuclide == p->event_nuclide_ || i_nuclide == -1) + if (i_nuclide == p.event_nuclide_ || i_nuclide == -1) score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index, filter_weight, i_nuclide, -1.0, flux); } @@ -2050,7 +2050,7 @@ void score_analog_tally_ce(Particle* p) // In the case that the user has requested to tally all nuclides, we // can take advantage of the fact that we know exactly how nuclide // bins correspond to nuclide indices. First, tally the nuclide. - auto i = p->event_nuclide_; + auto i = p.event_nuclide_; score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index, filter_weight, -1, -1.0, flux); @@ -2069,11 +2069,11 @@ void score_analog_tally_ce(Particle* p) } // Reset all the filter matches for the next tally event. - for (auto& match : p->filter_matches_) + for (auto& match : p.filter_matches_) match.bins_present_ = false; } -void score_analog_tally_mg(Particle* p) +void score_analog_tally_mg(Particle& p) { for (auto i_tally : model::active_analog_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2082,7 +2082,7 @@ void score_analog_tally_mg(Particle* p) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true, &p->filter_matches_); + auto end = FilterBinIter(tally, true, &p.filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2096,9 +2096,9 @@ void score_analog_tally_mg(Particle* p) double atom_density = 0.; if (i_nuclide >= 0) { - auto j = model::materials[p->material_]->mat_nuclide_index_[i_nuclide]; + auto j = model::materials[p.material_]->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; - atom_density = model::materials[p->material_]->atom_density_(j); + atom_density = model::materials[p.material_]->atom_density_(j); } score_general_mg(p, i_tally, i*tally.scores_.size(), filter_index, @@ -2114,15 +2114,15 @@ void score_analog_tally_mg(Particle* p) } // Reset all the filter matches for the next tally event. - for (auto& match : p->filter_matches_) + for (auto& match : p.filter_matches_) match.bins_present_ = false; } void -score_tracklength_tally(Particle* p, double distance) +score_tracklength_tally(Particle& p, double distance) { // Determine the tracklength estimate of the flux - double flux = p->wgt_ * distance; + double flux = p.wgt_ * distance; for (auto i_tally : model::active_tracklength_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2131,7 +2131,7 @@ score_tracklength_tally(Particle* p, double distance) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true, &p->filter_matches_); + auto end = FilterBinIter(tally, true, &p.filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2141,7 +2141,7 @@ score_tracklength_tally(Particle* p, double distance) // Loop over nuclide bins. if (tally.all_nuclides_) { - if (p->material_ != MATERIAL_VOID) + if (p.material_ != MATERIAL_VOID) score_all_nuclides(p, i_tally, flux*filter_weight, filter_index, filter_weight); @@ -2151,10 +2151,10 @@ score_tracklength_tally(Particle* p, double distance) double atom_density = 0.; if (i_nuclide >= 0) { - if (p->material_ != MATERIAL_VOID) { - auto j = model::materials[p->material_]->mat_nuclide_index_[i_nuclide]; + if (p.material_ != MATERIAL_VOID) { + auto j = model::materials[p.material_]->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; - atom_density = model::materials[p->material_]->atom_density_(j); + atom_density = model::materials[p.material_]->atom_density_(j); } } @@ -2179,19 +2179,19 @@ score_tracklength_tally(Particle* p, double distance) } // Reset all the filter matches for the next tally event. - for (auto& match : p->filter_matches_) + for (auto& match : p.filter_matches_) match.bins_present_ = false; } -void score_collision_tally(Particle* p) +void score_collision_tally(Particle& p) { // Determine the collision estimate of the flux double flux = 0.0; - if (p->type_ == Particle::Type::neutron || p->type_ == Particle::Type::photon) { + if (p.type_ == Particle::Type::neutron || p.type_ == Particle::Type::photon) { if (!settings::survival_biasing) { - flux = p->wgt_last_ / p->macro_xs_.total; + flux = p.wgt_last_ / p.macro_xs_.total; } else { - flux = (p->wgt_last_ + p->wgt_absorb_) / p->macro_xs_.total; + flux = (p.wgt_last_ + p.wgt_absorb_) / p.macro_xs_.total; } } @@ -2202,7 +2202,7 @@ void score_collision_tally(Particle* p) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true, &p->filter_matches_); + auto end = FilterBinIter(tally, true, &p.filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2221,9 +2221,9 @@ void score_collision_tally(Particle* p) double atom_density = 0.; if (i_nuclide >= 0) { - auto j = model::materials[p->material_]->mat_nuclide_index_[i_nuclide]; + auto j = model::materials[p.material_]->mat_nuclide_index_[i_nuclide]; if (j == C_NONE) continue; - atom_density = model::materials[p->material_]->atom_density_(j); + atom_density = model::materials[p.material_]->atom_density_(j); } //TODO: consider replacing this "if" with pointers or templates @@ -2246,15 +2246,15 @@ void score_collision_tally(Particle* p) } // Reset all the filter matches for the next tally event. - for (auto& match : p->filter_matches_) + for (auto& match : p.filter_matches_) match.bins_present_ = false; } void -score_surface_tally(Particle* p, const std::vector& tallies) +score_surface_tally(Particle& p, const std::vector& tallies) { // No collision, so no weight change when survival biasing - double flux = p->wgt_; + double flux = p.wgt_; for (auto i_tally : tallies) { auto& tally {*model::tallies[i_tally]}; @@ -2263,7 +2263,7 @@ score_surface_tally(Particle* p, const std::vector& tallies) // no valid combinations, use a continue statement to ensure we skip the // assume_separate break below. auto filter_iter = FilterBinIter(tally, p); - auto end = FilterBinIter(tally, true, &p->filter_matches_); + auto end = FilterBinIter(tally, true, &p.filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2290,7 +2290,7 @@ score_surface_tally(Particle* p, const std::vector& tallies) } // Reset all the filter matches for the next tally event. - for (auto& match : p->filter_matches_) + for (auto& match : p.filter_matches_) match.bins_present_ = false; }