Merge branch 'develop' into rng-stream-changes

src/physics.cpp

@@ -38,46 +38,46 @@ namespace openmc {
 void collision(Particle& p)
 {
   // Add to collision counter for particle
-  ++(p.n_collision_);
+  ++(p.n_collision());
 
   // Sample reaction for the material the particle is in
-  switch (p.type_) {
-  case Particle::Type::neutron:
+  switch (p.type()) {
+  case ParticleType::neutron:
     sample_neutron_reaction(p);
     break;
-  case Particle::Type::photon:
+  case ParticleType::photon:
     sample_photon_reaction(p);
     break;
-  case Particle::Type::electron:
+  case ParticleType::electron:
     sample_electron_reaction(p);
     break;
-  case Particle::Type::positron:
+  case ParticleType::positron:
     sample_positron_reaction(p);
     break;
   }
 
   // Kill particle if energy falls below cutoff
-  int type = static_cast<int>(p.type_);
-  if (p.E_ < settings::energy_cutoff[type]) {
-    p.alive_ = false;
-    p.wgt_ = 0.0;
+  int type = static_cast<int>(p.type());
+  if (p.E() < settings::energy_cutoff[type]) {
+    p.alive() = false;
+    p.wgt() = 0.0;
   }
 
   // Display information about collision
-  if (settings::verbosity >= 10 || p.trace_) {
+  if (settings::verbosity >= 10 || p.trace()) {
     std::string msg;
-    if (p.event_ == TallyEvent::KILL) {
-      msg = fmt::format("    Killed. Energy = {} eV.", p.E_);
-    } else if (p.type_ == Particle::Type::neutron) {
+    if (p.event() == TallyEvent::KILL) {
+      msg = fmt::format("    Killed. Energy = {} eV.", p.E());
+    } else if (p.type() == ParticleType::neutron) {
       msg = fmt::format("    {} with {}. Energy = {} eV.",
-        reaction_name(p.event_mt_), data::nuclides[p.event_nuclide_]->name_,
-        p.E_);
-    } else if (p.type_ == Particle::Type::photon) {
+        reaction_name(p.event_mt()), data::nuclides[p.event_nuclide()]->name_,
+        p.E());
+    } else if (p.type() == ParticleType::photon) {
       msg = fmt::format("    {} with {}. Energy = {} eV.",
-        reaction_name(p.event_mt_),
-        to_element(data::nuclides[p.event_nuclide_]->name_), p.E_);
+        reaction_name(p.event_mt()),
+        to_element(data::nuclides[p.event_nuclide()]->name_), p.E());
     } else {
-      msg = fmt::format("    Disappeared. Energy = {} eV.", p.E_);
+      msg = fmt::format("    Disappeared. Energy = {} eV.", p.E());
     }
     write_message(msg, 1);
   }
@@ -89,7 +89,7 @@ void sample_neutron_reaction(Particle& p)
   int i_nuclide = sample_nuclide(p);
 
   // Save which nuclide particle had collision with
-  p.event_nuclide_ = i_nuclide;
+  p.event_nuclide() = i_nuclide;
 
   // Create fission bank sites. Note that while a fission reaction is sampled,
   // it never actually "happens", i.e. the weight of the particle does not
@@ -108,7 +108,7 @@ void sample_neutron_reaction(Particle& p)
 
       // Make sure particle population doesn't grow out of control for
       // subcritical multiplication problems.
-      if (p.secondary_bank_.size() >= 10000) {
+      if (p.secondary_bank().size() >= 10000) {
         fatal_error("The secondary particle bank appears to be growing without "
         "bound. You are likely running a subcritical multiplication problem "
         "with k-effective close to or greater than one.");
@@ -124,28 +124,30 @@ void sample_neutron_reaction(Particle& p)
   // If survival biasing is being used, the following subroutine adjusts the
   // weight of the particle. Otherwise, it checks to see if absorption occurs
 
-  if (p.neutron_xs_[i_nuclide].absorption > 0.0) {
+  if (p.neutron_xs(i_nuclide).absorption > 0.0) {
     absorption(p, i_nuclide);
   } else {
-    p.wgt_absorb_ = 0.0;
+    p.wgt_absorb() = 0.0;
   }
-  if (!p.alive_) return;
+  if (!p.alive())
+    return;
 
   // Sample a scattering reaction and determine the secondary energy of the
   // exiting neutron
   scatter(p, i_nuclide);
 
   // Advance URR seed stream 'N' times after energy changes
-  if (p.E_ != p.E_last_) {
-    p.stream_ = STREAM_URR_PTABLE;
+  if (p.E() != p.E_last()) {
+    p.stream() = STREAM_URR_PTABLE;
     advance_prn_seed(data::nuclides.size(), p.current_seed());
-    p.stream_ = STREAM_TRACKING;
+    p.stream() = STREAM_TRACKING;
   }
 
   // Play russian roulette if survival biasing is turned on
   if (settings::survival_biasing) {
     russian_roulette(p);
-    if (!p.alive_) return;
+    if (!p.alive())
+      return;
   }
 }
 
@@ -157,8 +159,9 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
   double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
 
   // Determine the expected number of neutrons produced
-  double nu_t = p.wgt_ / simulation::keff * weight * p.neutron_xs_[
-    i_nuclide].nu_fission / p.neutron_xs_[i_nuclide].total;
+  double nu_t = p.wgt() / simulation::keff * weight *
+                p.neutron_xs(i_nuclide).nu_fission /
+                p.neutron_xs(i_nuclide).total;
 
   // Sample the number of neutrons produced
   int nu = static_cast<int>(nu_t);
@@ -173,9 +176,9 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
   double nu_d[MAX_DELAYED_GROUPS] = {0.};
 
   // Clear out particle's nu fission bank
-  p.nu_bank_.clear();
+  p.nu_bank().clear();
 
-  p.fission_ = true;
+  p.fission() = true;
   int skipped = 0;
 
   // Determine whether to place fission sites into the shared fission bank
@@ -184,16 +187,16 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
 
   for (int i = 0; i < nu; ++i) {
     // Initialize fission site object with particle data
-    Particle::Bank site;
+    SourceSite site;
     site.r = p.r();
-    site.particle = Particle::Type::neutron;
+    site.particle = ParticleType::neutron;
     site.wgt = 1. / weight;
-    site.parent_id = p.id_;
-    site.progeny_id = p.n_progeny_++;
+    site.parent_id = p.id();
+    site.progeny_id = p.n_progeny()++;
     site.surf_id = 0;
 
     // Sample delayed group and angle/energy for fission reaction
-    sample_fission_neutron(i_nuclide, rx, p.E_, &site, p.current_seed());
+    sample_fission_neutron(i_nuclide, rx, p.E(), &site, p.current_seed());
 
     // Store fission site in bank
     if (use_fission_bank) {
@@ -205,20 +208,20 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
         break;
       }
     } else {
-      p.secondary_bank_.push_back(site);
+      p.secondary_bank().push_back(site);
     }
 
     // Set the delayed group on the particle as well
-    p.delayed_group_ = site.delayed_group;
+    p.delayed_group() = site.delayed_group;
 
     // Increment the number of neutrons born delayed
-    if (p.delayed_group_ > 0) {
-      nu_d[p.delayed_group_-1]++;
+    if (p.delayed_group() > 0) {
+      nu_d[p.delayed_group() - 1]++;
     }
 
     // Write fission particles to nuBank
-    p.nu_bank_.emplace_back();
-    Particle::NuBank* nu_bank_entry = &p.nu_bank_.back();
+    p.nu_bank().emplace_back();
+    NuBank* nu_bank_entry = &p.nu_bank().back();
     nu_bank_entry->wgt              = site.wgt;
     nu_bank_entry->E                = site.E;
     nu_bank_entry->delayed_group    = site.delayed_group;
@@ -227,7 +230,7 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
   // If shared fission bank was full, and no fissions could be added,
   // set the particle fission flag to false.
   if (nu == skipped) {
-    p.fission_ = false;
+    p.fission() = false;
     return;
   }
 
@@ -236,10 +239,10 @@ create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
   nu -= skipped;
 
   // Store the total weight banked for analog fission tallies
-  p.n_bank_ = nu;
-  p.wgt_bank_ = nu / weight;
+  p.n_bank() = nu;
+  p.wgt_bank() = nu / weight;
   for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {
-    p.n_delayed_bank_[d] = nu_d[d];
+    p.n_delayed_bank(d) = nu_d[d];
   }
 }
 
@@ -248,20 +251,20 @@ void sample_photon_reaction(Particle& p)
   // Kill photon if below energy cutoff -- an extra check is made here because
   // photons with energy below the cutoff may have been produced by neutrons
   // reactions or atomic relaxation
-  int photon = static_cast<int>(Particle::Type::photon);
-  if (p.E_ < settings::energy_cutoff[photon]) {
-    p.E_ = 0.0;
-    p.alive_ = false;
+  int photon = static_cast<int>(ParticleType::photon);
+  if (p.E() < settings::energy_cutoff[photon]) {
+    p.E() = 0.0;
+    p.alive() = false;
     return;
   }
 
   // Sample element within material
   int i_element = sample_element(p);
-  const auto& micro {p.photon_xs_[i_element]};
+  const auto& micro {p.photon_xs(i_element)};
   const auto& element {*data::elements[i_element]};
 
   // Calculate photon energy over electron rest mass equivalent
-  double alpha = p.E_/MASS_ELECTRON_EV;
+  double alpha = p.E() / MASS_ELECTRON_EV;
 
   // For tallying purposes, this routine might be called directly. In that
   // case, we need to sample a reaction via the cutoff variable
@@ -273,8 +276,8 @@ void sample_photon_reaction(Particle& p)
   if (prob > cutoff) {
     double mu = element.rayleigh_scatter(alpha, p.current_seed());
     p.u() = rotate_angle(p.u(), mu, nullptr, p.current_seed());
-    p.event_ = TallyEvent::SCATTER;
-    p.event_mt_ = COHERENT;
+    p.event() = TallyEvent::SCATTER;
+    p.event_mt() = COHERENT;
     return;
   }
 
@@ -297,12 +300,12 @@ void sample_photon_reaction(Particle& p)
     // Create Compton electron
     double phi = uniform_distribution(0., 2.0*PI, p.current_seed());
     double E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b;
-    int electron = static_cast<int>(Particle::Type::electron);
+    int electron = static_cast<int>(ParticleType::electron);
     if (E_electron >= settings::energy_cutoff[electron]) {
       double mu_electron = (alpha - alpha_out*mu)
         / std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu);
       Direction u = rotate_angle(p.u(), mu_electron, &phi, p.current_seed());
-      p.create_secondary(p.wgt_, u, E_electron, Particle::Type::electron);
+      p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron);
     }
 
     // TODO: Compton subshell data does not match atomic relaxation data
@@ -314,10 +317,10 @@ void sample_photon_reaction(Particle& p)
     }
 
     phi += PI;
-    p.E_ = alpha_out*MASS_ELECTRON_EV;
+    p.E() = alpha_out * MASS_ELECTRON_EV;
     p.u() = rotate_angle(p.u(), mu, &phi, p.current_seed());
-    p.event_ = TallyEvent::SCATTER;
-    p.event_mt_ = INCOHERENT;
+    p.event() = TallyEvent::SCATTER;
+    p.event_mt() = INCOHERENT;
     return;
   }
 
@@ -340,7 +343,7 @@ void sample_photon_reaction(Particle& p)
 
       prob += xs;
       if (prob > cutoff) {
-        double E_electron = p.E_ - shell.binding_energy;
+        double E_electron = p.E() - shell.binding_energy;
 
         // Sample mu using non-relativistic Sauter distribution.
         // See Eqns 3.19 and 3.20 in "Implementing a photon physics
@@ -363,15 +366,15 @@ void sample_photon_reaction(Particle& p)
         u.z = std::sqrt(1.0 - mu*mu)*std::sin(phi);
 
         // Create secondary electron
-        p.create_secondary(p.wgt_, u, E_electron, Particle::Type::electron);
+        p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron);
 
         // Allow electrons to fill orbital and produce auger electrons
         // and fluorescent photons
         element.atomic_relaxation(shell, p);
-        p.event_ = TallyEvent::ABSORB;
-        p.event_mt_ = 533 + shell.index_subshell;
-        p.alive_ = false;
-        p.E_ = 0.0;
+        p.event() = TallyEvent::ABSORB;
+        p.event_mt() = 533 + shell.index_subshell;
+        p.alive() = false;
+        p.E() = 0.0;
         return;
       }
     }
@@ -388,16 +391,16 @@ void sample_photon_reaction(Particle& p)
 
     // Create secondary electron
     Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed());
-    p.create_secondary(p.wgt_, u, E_electron, Particle::Type::electron);
+    p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron);
 
     // Create secondary positron
     u = rotate_angle(p.u(), mu_positron, nullptr, p.current_seed());
-    p.create_secondary(p.wgt_, u, E_positron, Particle::Type::positron);
+    p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron);
 
-    p.event_ = TallyEvent::ABSORB;
-    p.event_mt_ = PAIR_PROD;
-    p.alive_ = false;
-    p.E_ = 0.0;
+    p.event() = TallyEvent::ABSORB;
+    p.event_mt() = PAIR_PROD;
+    p.alive() = false;
+    p.E() = 0.0;
   }
 }
 
@@ -410,9 +413,9 @@ void sample_electron_reaction(Particle& p)
     thick_target_bremsstrahlung(p, &E_lost);
   }
 
-  p.E_ = 0.0;
-  p.alive_ = false;
-  p.event_ = TallyEvent::ABSORB;
+  p.E() = 0.0;
+  p.alive() = false;
+  p.event() = TallyEvent::ABSORB;
 }
 
 void sample_positron_reaction(Particle& p)
@@ -428,21 +431,21 @@ void sample_positron_reaction(Particle& p)
   Direction u = isotropic_direction(p.current_seed());
 
   // Create annihilation photon pair traveling in opposite directions
-  p.create_secondary(p.wgt_, u, MASS_ELECTRON_EV, Particle::Type::photon);
-  p.create_secondary(p.wgt_, -u, MASS_ELECTRON_EV, Particle::Type::photon);
+  p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon);
+  p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon);
 
-  p.E_ = 0.0;
-  p.alive_ = false;
-  p.event_ = TallyEvent::ABSORB;
+  p.E() = 0.0;
+  p.alive() = false;
+  p.event() = TallyEvent::ABSORB;
 }
 
 int sample_nuclide(Particle& p)
 {
   // Sample cumulative distribution function
-  double cutoff = prn(p.current_seed()) * p.macro_xs_.total;
+  double cutoff = prn(p.current_seed()) * p.macro_xs().total;
 
   // Get pointers to nuclide/density arrays
-  const auto& mat {model::materials[p.material_]};
+  const auto& mat {model::materials[p.material()]};
   int n = mat->nuclide_.size();
 
   double prob = 0.0;
@@ -452,7 +455,7 @@ int sample_nuclide(Particle& p)
     double atom_density = mat->atom_density_[i];
 
     // Increment probability to compare to cutoff
-    prob += atom_density * p.neutron_xs_[i_nuclide].total;
+    prob += atom_density * p.neutron_xs(i_nuclide).total;
     if (prob >= cutoff) return i_nuclide;
   }
 
@@ -464,10 +467,10 @@ int sample_nuclide(Particle& p)
 int sample_element(Particle& p)
 {
   // Sample cumulative distribution function
-  double cutoff = prn(p.current_seed()) * p.macro_xs_.total;
+  double cutoff = prn(p.current_seed()) * p.macro_xs().total;
 
   // Get pointers to elements, densities
-  const auto& mat {model::materials[p.material_]};
+  const auto& mat {model::materials[p.material()]};
 
   double prob = 0.0;
   for (int i = 0; i < mat->element_.size(); ++i) {
@@ -476,13 +479,13 @@ int sample_element(Particle& p)
     double atom_density = mat->atom_density_[i];
 
     // Determine microscopic cross section
-    double sigma = atom_density * p.photon_xs_[i_element].total;
+    double sigma = atom_density * p.photon_xs(i_element).total;
 
     // Increment probability to compare to cutoff
     prob += sigma;
     if (prob > cutoff) {
       // Save which nuclide particle had collision with for tally purpose
-      p.event_nuclide_ = mat->nuclide_[i];
+      p.event_nuclide() = mat->nuclide_[i];
 
       return i_element;
     }
@@ -501,23 +504,23 @@ Reaction& sample_fission(int i_nuclide, Particle& p)
   // If we're in the URR, by default use the first fission reaction. We also
   // default to the first reaction if we know that there are no partial fission
   // reactions
-  if (p.neutron_xs_[i_nuclide].use_ptable || !nuc->has_partial_fission_) {
+  if (p.neutron_xs(i_nuclide).use_ptable || !nuc->has_partial_fission_) {
     return *nuc->fission_rx_[0];
   }
 
   // Check to see if we are in a windowed multipole range.  WMP only supports
   // the first fission reaction.
   if (nuc->multipole_) {
-    if (p.E_ >= nuc->multipole_->E_min_ && p.E_ <= nuc->multipole_->E_max_) {
+    if (p.E() >= nuc->multipole_->E_min_ && p.E() <= nuc->multipole_->E_max_) {
       return *nuc->fission_rx_[0];
     }
   }
 
   // Get grid index and interpolatoin factor and sample fission cdf
-  int i_temp = p.neutron_xs_[i_nuclide].index_temp;
-  int i_grid = p.neutron_xs_[i_nuclide].index_grid;
-  double f = p.neutron_xs_[i_nuclide].interp_factor;
-  double cutoff = prn(p.current_seed()) * p.neutron_xs_[i_nuclide].fission;
+  int i_temp = p.neutron_xs(i_nuclide).index_temp;
+  int i_grid = p.neutron_xs(i_nuclide).index_grid;
+  double f = p.neutron_xs(i_nuclide).interp_factor;
+  double cutoff = prn(p.current_seed()) * p.neutron_xs(i_nuclide).fission;
   double prob = 0.0;
 
   // Loop through each partial fission reaction type
@@ -541,10 +544,10 @@ Reaction& sample_fission(int i_nuclide, Particle& p)
 void sample_photon_product(int i_nuclide, Particle& p, int* i_rx, int* i_product)
 {
   // Get grid index and interpolation factor and sample photon production cdf
-  int i_temp = p.neutron_xs_[i_nuclide].index_temp;
-  int i_grid = p.neutron_xs_[i_nuclide].index_grid;
-  double f = p.neutron_xs_[i_nuclide].interp_factor;
-  double cutoff = prn(p.current_seed()) * p.neutron_xs_[i_nuclide].photon_prod;
+  int i_temp = p.neutron_xs(i_nuclide).index_temp;
+  int i_grid = p.neutron_xs(i_nuclide).index_grid;
+  double f = p.neutron_xs(i_nuclide).interp_factor;
+  double cutoff = prn(p.current_seed()) * p.neutron_xs(i_nuclide).photon_prod;
   double prob = 0.0;
 
   // Loop through each reaction type
@@ -561,22 +564,22 @@ void sample_photon_product(int i_nuclide, Particle& p, int* i_rx, int* i_product
       + f*(rx->xs_[i_temp].value[i_grid - threshold + 1]));
 
     for (int j = 0; j < rx->products_.size(); ++j) {
-      if (rx->products_[j].particle_ == Particle::Type::photon) {
+      if (rx->products_[j].particle_ == ParticleType::photon) {
         // For fission, artificially increase the photon yield to account
         // for delayed photons
         double f = 1.0;
         if (settings::delayed_photon_scaling) {
           if (is_fission(rx->mt_)) {
             if (nuc->prompt_photons_ && nuc->delayed_photons_) {
-              double energy_prompt = (*nuc->prompt_photons_)(p.E_);
-              double energy_delayed = (*nuc->delayed_photons_)(p.E_);
+              double energy_prompt = (*nuc->prompt_photons_)(p.E());
+              double energy_delayed = (*nuc->delayed_photons_)(p.E());
               f = (energy_prompt + energy_delayed)/(energy_prompt);
             }
           }
         }
 
         // add to cumulative probability
-        prob += f * (*rx->products_[j].yield_)(p.E_) * xs;
+        prob += f * (*rx->products_[j].yield_)(p.E()) * xs;
 
         *i_rx = i;
         *i_product = j;
@@ -590,31 +593,33 @@ void absorption(Particle& p, int i_nuclide)
 {
   if (settings::survival_biasing) {
     // Determine weight absorbed in survival biasing
-    p.wgt_absorb_ = p.wgt_ * p.neutron_xs_[i_nuclide].absorption /
-          p.neutron_xs_[i_nuclide].total;
+    p.wgt_absorb() = p.wgt() * p.neutron_xs(i_nuclide).absorption /
+                     p.neutron_xs(i_nuclide).total;
 
     // Adjust weight of particle by probability of absorption
-    p.wgt_ -= p.wgt_absorb_;
-    p.wgt_last_ = p.wgt_;
+    p.wgt() -= p.wgt_absorb();
+    p.wgt_last() = p.wgt();
 
     // Score implicit absorption estimate of keff
     if (settings::run_mode == RunMode::EIGENVALUE) {
-      p.keff_tally_absorption_ += p.wgt_absorb_ * p.neutron_xs_[
-        i_nuclide].nu_fission / p.neutron_xs_[i_nuclide].absorption;
+      p.keff_tally_absorption() += p.wgt_absorb() *
+                                   p.neutron_xs(i_nuclide).nu_fission /
+                                   p.neutron_xs(i_nuclide).absorption;
     }
   } else {
     // See if disappearance reaction happens
-    if (p.neutron_xs_[i_nuclide].absorption >
-        prn(p.current_seed()) * p.neutron_xs_[i_nuclide].total) {
+    if (p.neutron_xs(i_nuclide).absorption >
+        prn(p.current_seed()) * p.neutron_xs(i_nuclide).total) {
       // Score absorption estimate of keff
       if (settings::run_mode == RunMode::EIGENVALUE) {
-        p.keff_tally_absorption_ += p.wgt_ * p.neutron_xs_[
-          i_nuclide].nu_fission / p.neutron_xs_[i_nuclide].absorption;
+        p.keff_tally_absorption() += p.wgt() *
+                                     p.neutron_xs(i_nuclide).nu_fission /
+                                     p.neutron_xs(i_nuclide).absorption;
       }
 
-      p.alive_ = false;
-      p.event_ = TallyEvent::ABSORB;
-      p.event_mt_ = N_DISAPPEAR;
+      p.alive() = false;
+      p.event() = TallyEvent::ABSORB;
+      p.event_mt() = N_DISAPPEAR;
     }
   }
 }
@@ -626,7 +631,7 @@ void scatter(Particle& p, int i_nuclide)
 
   // Get pointer to nuclide and grid index/interpolation factor
   const auto& nuc {data::nuclides[i_nuclide]};
-  const auto& micro {p.neutron_xs_[i_nuclide]};
+  const auto& micro {p.neutron_xs(i_nuclide)};
   int i_temp =  micro.index_temp;
   int i_grid =  micro.index_grid;
   double f = micro.interp_factor;
@@ -647,12 +652,12 @@ void scatter(Particle& p, int i_nuclide)
     // NON-S(A,B) ELASTIC SCATTERING
 
     // Determine temperature
-    double kT = nuc->multipole_ ? p.sqrtkT_*p.sqrtkT_ : nuc->kTs_[i_temp];
+    double kT = nuc->multipole_ ? p.sqrtkT() * p.sqrtkT() : nuc->kTs_[i_temp];
 
     // Perform collision physics for elastic scattering
     elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p);
 
-    p.event_mt_ = ELASTIC;
+    p.event_mt() = ELASTIC;
     sampled = true;
   }
 
@@ -663,7 +668,7 @@ void scatter(Particle& p, int i_nuclide)
 
     sab_scatter(i_nuclide, micro.index_sab, p);
 
-    p.event_mt_ = ELASTIC;
+    p.event_mt() = ELASTIC;
     sampled = true;
   }
 
@@ -695,20 +700,20 @@ void scatter(Particle& p, int i_nuclide)
     // Perform collision physics for inelastic scattering
     const auto& rx {nuc->reactions_[i]};
     inelastic_scatter(*nuc, *rx, p);
-    p.event_mt_ = rx->mt_;
+    p.event_mt() = rx->mt_;
   }
 
   // Set event component
-  p.event_ = TallyEvent::SCATTER;
+  p.event() = TallyEvent::SCATTER;
 
   // Sample new outgoing angle for isotropic-in-lab scattering
-  const auto& mat {model::materials[p.material_]};
+  const auto& mat {model::materials[p.material()]};
   if (!mat->p0_.empty()) {
     int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide];
     if (mat->p0_[i_nuc_mat]) {
       // Sample isotropic-in-lab outgoing direction
       p.u() = isotropic_direction(p.current_seed());
-      p.mu_ = u_old.dot(p.u());
+      p.mu() = u_old.dot(p.u());
     }
   }
 }
@@ -719,7 +724,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
   // get pointer to nuclide
   const auto& nuc {data::nuclides[i_nuclide]};
 
-  double vel = std::sqrt(p.E_);
+  double vel = std::sqrt(p.E());
   double awr = nuc->awr_;
 
   // Neutron velocity in LAB
@@ -727,9 +732,9 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
 
   // Sample velocity of target nucleus
   Direction v_t {};
-  if (!p.neutron_xs_[i_nuclide].use_ptable) {
-    v_t = sample_target_velocity(*nuc, p.E_, p.u(), v_n,
-      p.neutron_xs_[i_nuclide].elastic, kT, p.current_seed());
+  if (!p.neutron_xs(i_nuclide).use_ptable) {
+    v_t = sample_target_velocity(*nuc, p.E(), p.u(), v_n,
+      p.neutron_xs(i_nuclide).elastic, kT, p.current_seed());
   }
 
   // Velocity of center-of-mass
@@ -747,7 +752,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
   auto& d = rx.products_[0].distribution_[0];
   auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
   if (!d_->angle().empty()) {
-    mu_cm = d_->angle().sample(p.E_, p.current_seed());
+    mu_cm = d_->angle().sample(p.E(), p.current_seed());
   } else {
     mu_cm = uniform_distribution(-1., 1., p.current_seed());
   }
@@ -763,12 +768,12 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
   // Transform back to LAB frame
   v_n += v_cm;
 
-  p.E_ = v_n.dot(v_n);
-  vel = std::sqrt(p.E_);
+  p.E() = v_n.dot(v_n);
+  vel = std::sqrt(p.E());
 
   // compute cosine of scattering angle in LAB frame by taking dot product of
   // neutron's pre- and post-collision angle
-  p.mu_ = p.u().dot(v_n) / vel;
+  p.mu() = p.u().dot(v_n) / vel;
 
   // Set energy and direction of particle in LAB frame
   p.u() = v_n / vel;
@@ -776,22 +781,24 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
   // Because of floating-point roundoff, it may be possible for mu_lab to be
   // outside of the range [-1,1). In these cases, we just set mu_lab to exactly
   // -1 or 1
-  if (std::abs(p.mu_) > 1.0) p.mu_ = std::copysign(1.0, p.mu_);
+  if (std::abs(p.mu()) > 1.0)
+    p.mu() = std::copysign(1.0, p.mu());
 }
 
 void sab_scatter(int i_nuclide, int i_sab, Particle& p)
 {
   // Determine temperature index
-  const auto& micro {p.neutron_xs_[i_nuclide]};
+  const auto& micro {p.neutron_xs(i_nuclide)};
   int i_temp = micro.index_temp_sab;
 
   // Sample energy and angle
   double E_out;
-  data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, p.E_, &E_out, &p.mu_, p.current_seed());
+  data::thermal_scatt[i_sab]->data_[i_temp].sample(
+    micro, p.E(), &E_out, &p.mu(), p.current_seed());
 
   // Set energy to outgoing, change direction of particle
-  p.E_ = E_out;
-  p.u() = rotate_angle(p.u(), p.mu_, nullptr, p.current_seed());
+  p.E() = E_out;
+  p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed());
 }
 
 Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u,
@@ -996,7 +1003,8 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_
   return vt * rotate_angle(u, mu, nullptr, seed);
 }
 
-void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Particle::Bank* site, uint64_t* seed)
+void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in,
+  SourceSite* site, uint64_t* seed)
 {
   // Determine total nu, delayed nu, and delayed neutron fraction
   const auto& nuc {data::nuclides[i_nuclide]};
@@ -1044,7 +1052,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Part
     rx.products_[site->delayed_group].sample(E_in, site->E, mu, seed);
 
     // resample if energy is greater than maximum neutron energy
-    constexpr int neutron = static_cast<int>(Particle::Type::neutron);
+    constexpr int neutron = static_cast<int>(ParticleType::neutron);
     if (site->E < data::energy_max[neutron]) break;
 
     // check for large number of resamples
@@ -1065,7 +1073,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction& rx, double E_in, Part
 void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
 {
   // copy energy of neutron
-  double E_in = p.E_;
+  double E_in = p.E();
 
   // sample outgoing energy and scattering cosine
   double E;
@@ -1092,8 +1100,8 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
   if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu);
 
   // Set outgoing energy and scattering angle
-  p.E_ = E;
-  p.mu_ = mu;
+  p.E() = E;
+  p.mu() = mu;
 
   // change direction of particle
   p.u() = rotate_angle(p.u(), mu, nullptr, p.current_seed());
@@ -1103,19 +1111,19 @@ void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
   if (std::floor(yield) == yield) {
     // If yield is integral, create exactly that many secondary particles
     for (int i = 0; i < static_cast<int>(std::round(yield)) - 1; ++i) {
-      p.create_secondary(p.wgt_, p.u(), p.E_, Particle::Type::neutron);
+      p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron);
     }
   } else {
     // Otherwise, change weight of particle based on yield
-    p.wgt_ *= yield;
+    p.wgt() *= yield;
   }
 }
 
 void sample_secondary_photons(Particle& p, int i_nuclide)
 {
   // Sample the number of photons produced
-  double y_t = p.neutron_xs_[i_nuclide].photon_prod /
-    p.neutron_xs_[i_nuclide].total;
+  double y_t =
+    p.neutron_xs(i_nuclide).photon_prod / p.neutron_xs(i_nuclide).total;
   int y = static_cast<int>(y_t);
   if (prn(p.current_seed()) <= y_t - y) ++y;
 
@@ -1130,7 +1138,7 @@ void sample_secondary_photons(Particle& p, int i_nuclide)
     auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx];
     double E;
     double mu;
-    rx->products_[i_product].sample(p.E_, E, mu, p.current_seed());
+    rx->products_[i_product].sample(p.E(), E, mu, p.current_seed());
 
     // Sample the new direction
     Direction u = rotate_angle(p.u(), mu, nullptr, p.current_seed());
@@ -1142,14 +1150,13 @@ void sample_secondary_photons(Particle& p, int i_nuclide)
     // calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020.
     double wgt;
     if (settings::run_mode == RunMode::EIGENVALUE && !is_fission(rx->mt_)) {
-      wgt = simulation::keff * p.wgt_;
+      wgt = simulation::keff * p.wgt();
     } else {
-      wgt = p.wgt_;
+      wgt = p.wgt();
     }
 
     // Create the secondary photon
-    p.create_secondary(wgt, u, E, Particle::Type::photon);
-
+    p.create_secondary(wgt, u, E, ParticleType::photon);
   }
 }