diff --git a/src/output.cpp b/src/output.cpp
index 35cd9f9aaa..97436dc244 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -657,8 +657,7 @@ write_tallies()
       case DIFF_NUCLIDE_DENSITY:
         tallies_out << " Nuclide density derivative  Material "
           << std::to_string(deriv.diff_material) << "  Nuclide "
-          // TODO: off-by-one
-          << data::nuclides[deriv.diff_nuclide-1]->name_ << "\n";
+          << data::nuclides[deriv.diff_nuclide]->name_ << "\n";
         break;
       case DIFF_TEMPERATURE:
         tallies_out << " Temperature derivative  Material "
diff --git a/src/physics.cpp b/src/physics.cpp
index b93c2c232e..23ce8f5603 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -63,10 +63,10 @@ void collision(Particle* p)
     std::stringstream msg;
     if (static_cast<ParticleType>(p->type) == ParticleType::neutron) {
       msg << "    " << reaction_name(p->event_MT) << " with " <<
-        data::nuclides[p->event_nuclide-1]->name_ << ". Energy = " << p->E << " eV.";
+        data::nuclides[p->event_nuclide]->name_ << ". Energy = " << p->E << " eV.";
     } else {
       msg << "    " << reaction_name(p->event_MT) << " with " <<
-        data::elements[p->event_nuclide-1].name_ << ". Energy = " << p->E << " eV.";
+        data::elements[p->event_nuclide].name_ << ". Energy = " << p->E << " eV.";
     }
     write_message(msg, 1);
   }
@@ -78,8 +78,7 @@ void sample_neutron_reaction(Particle* p)
   int i_nuclide = sample_nuclide(p);
 
   // Save which nuclide particle had collision with
-  // TODO: off-by-one
-  p->event_nuclide = i_nuclide + 1;
+  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
@@ -228,8 +227,7 @@ void sample_photon_reaction(Particle* p)
 
   // Sample element within material
   int i_element = sample_element(p);
-  // TODO: off-by-one
-  p->event_nuclide = i_element + 1;
+  p->event_nuclide = i_element;
   const auto& micro {simulation::micro_photon_xs[i_element]};
   const auto& element {data::elements[i_element]};
 
diff --git a/src/state_point.cpp b/src/state_point.cpp
index 2a3f9d2bbc..3730c32163 100644
--- a/src/state_point.cpp
+++ b/src/state_point.cpp
@@ -125,7 +125,7 @@ openmc_statepoint_write(const char* filename, bool* write_source)
           write_dataset(deriv_group, "independent variable", "nuclide_density");
           //TODO: off-by-one
           write_dataset(deriv_group, "nuclide",
-            data::nuclides[deriv.diff_nuclide-1]->name_);
+            data::nuclides[deriv.diff_nuclide]->name_);
         } else if (deriv.variable == DIFF_TEMPERATURE) {
           write_dataset(deriv_group, "independent variable", "temperature");
         } else {
diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp
index 1944b39ddc..1098adc550 100644
--- a/src/tallies/derivative.cpp
+++ b/src/tallies/derivative.cpp
@@ -51,8 +51,7 @@ TallyDerivative::TallyDerivative(pugi::xml_node node)
     for (auto i = 0; i < data::nuclides.size(); ++i) {
       if (data::nuclides[i]->name_ == nuclide_name) {
         found = true;
-        //TODO: off-by-one
-        diff_nuclide = i + 1;
+        diff_nuclide = i;
       }
     }
     if (!found) {
@@ -188,7 +187,6 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
   // where i is the perturbed nuclide.
 
   case DIFF_NUCLIDE_DENSITY:
-    //TODO: off-by-one throughout on diff_nuclide
     switch (tally.estimator_) {
 
     case ESTIMATOR_ANALOG:
@@ -208,7 +206,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
           // Find the index of the perturbed nuclide.
           int i;
           for (i = 0; i < material.nuclide_.size(); ++i)
-            if (material.nuclide_[i] == deriv.diff_nuclide - 1) break;
+            if (material.nuclide_[i] == deriv.diff_nuclide) break;
           score *= flux_deriv + 1. / material.atom_density_(i);
         }
         break;
@@ -225,9 +223,9 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
       case SCORE_TOTAL:
         if (i_nuclide == -1 && simulation::material_xs.total > 0.0) {
           score *= flux_deriv
-            + simulation::micro_xs[deriv.diff_nuclide-1].total
+            + simulation::micro_xs[deriv.diff_nuclide].total
             / simulation::material_xs.total;
-        } else if (i_nuclide == deriv.diff_nuclide-1
+        } else if (i_nuclide == deriv.diff_nuclide
                    && simulation::micro_xs[i_nuclide].total) {
           score *= flux_deriv + 1. / atom_density;
         } else {
@@ -239,11 +237,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
         if (i_nuclide == -1 && (simulation::material_xs.total
                                 - simulation::material_xs.absorption) > 0.0) {
           score *= flux_deriv
-            + (simulation::micro_xs[deriv.diff_nuclide-1].total
-            - simulation::micro_xs[deriv.diff_nuclide-1].absorption)
+            + (simulation::micro_xs[deriv.diff_nuclide].total
+            - simulation::micro_xs[deriv.diff_nuclide].absorption)
             / (simulation::material_xs.total
             - simulation::material_xs.absorption);
-        } else if (i_nuclide == deriv.diff_nuclide-1) {
+        } else if (i_nuclide == deriv.diff_nuclide) {
           score *= flux_deriv + 1. / atom_density;
         } else {
           score *= flux_deriv;
@@ -253,9 +251,9 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
       case SCORE_ABSORPTION:
         if (i_nuclide == -1 && simulation::material_xs.absorption > 0.0) {
           score *= flux_deriv
-            + simulation::micro_xs[deriv.diff_nuclide-1].absorption
+            + simulation::micro_xs[deriv.diff_nuclide].absorption
             / simulation::material_xs.absorption;
-        } else if (i_nuclide == deriv.diff_nuclide-1
+        } else if (i_nuclide == deriv.diff_nuclide
                    && simulation::micro_xs[i_nuclide].absorption) {
           score *= flux_deriv + 1. / atom_density;
         } else {
@@ -266,9 +264,9 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
       case SCORE_FISSION:
         if (i_nuclide == -1 && simulation::material_xs.fission > 0.0) {
           score *= flux_deriv
-            + simulation::micro_xs[deriv.diff_nuclide-1].fission
+            + simulation::micro_xs[deriv.diff_nuclide].fission
             / simulation::material_xs.fission;
-        } else if (i_nuclide == deriv.diff_nuclide-1
+        } else if (i_nuclide == deriv.diff_nuclide
                    && simulation::micro_xs[i_nuclide].fission) {
           score *= flux_deriv + 1. / atom_density;
         } else {
@@ -279,9 +277,9 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
       case SCORE_NU_FISSION:
         if (i_nuclide == -1 && simulation::material_xs.nu_fission > 0.0) {
           score *= flux_deriv
-            + simulation::micro_xs[deriv.diff_nuclide-1].nu_fission
+            + simulation::micro_xs[deriv.diff_nuclide].nu_fission
             / simulation::material_xs.nu_fission;
-        } else if (i_nuclide == deriv.diff_nuclide-1
+        } else if (i_nuclide == deriv.diff_nuclide
                    && simulation::micro_xs[i_nuclide].nu_fission) {
           score *= flux_deriv + 1. / atom_density;
         } else {
@@ -322,9 +320,9 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
         // Find the index of the event nuclide.
         int i;
         for (i = 0; i < material.nuclide_.size(); ++i)
-          if (material.nuclide_[i] == p->event_nuclide-1) break;
+          if (material.nuclide_[i] == p->event_nuclide) break;
 
-        const auto& nuc {*data::nuclides[p->event_nuclide-1]};
+        const auto& nuc {*data::nuclides[p->event_nuclide]};
         if (!multipole_in_range(&nuc, p->last_E)) {
           score *= flux_deriv;
           break;
@@ -333,7 +331,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
         switch (score_bin) {
 
         case SCORE_TOTAL:
-          if (simulation::micro_xs[p->event_nuclide-1].total) {
+          if (simulation::micro_xs[p->event_nuclide].total) {
             double dsig_s, dsig_a, dsig_f;
             std::tie(dsig_s, dsig_a, dsig_f)
               = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
@@ -345,8 +343,8 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
           break;
 
         case SCORE_SCATTER:
-          if (simulation::micro_xs[p->event_nuclide-1].total
-              - simulation::micro_xs[p->event_nuclide-1].absorption) {
+          if (simulation::micro_xs[p->event_nuclide].total
+              - simulation::micro_xs[p->event_nuclide].absorption) {
             double dsig_s, dsig_a, dsig_f;
             std::tie(dsig_s, dsig_a, dsig_f)
               = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
@@ -359,7 +357,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
           break;
 
         case SCORE_ABSORPTION:
-          if (simulation::micro_xs[p->event_nuclide-1].absorption) {
+          if (simulation::micro_xs[p->event_nuclide].absorption) {
             double dsig_s, dsig_a, dsig_f;
             std::tie(dsig_s, dsig_a, dsig_f)
               = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
@@ -371,7 +369,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
           break;
 
         case SCORE_FISSION:
-          if (simulation::micro_xs[p->event_nuclide-1].fission) {
+          if (simulation::micro_xs[p->event_nuclide].fission) {
             double dsig_s, dsig_a, dsig_f;
             std::tie(dsig_s, dsig_a, dsig_f)
               = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
@@ -383,9 +381,9 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
           break;
 
         case SCORE_NU_FISSION:
-          if (simulation::micro_xs[p->event_nuclide-1].fission) {
-            double nu = simulation::micro_xs[p->event_nuclide-1].nu_fission
-              / simulation::micro_xs[p->event_nuclide-1].fission;
+          if (simulation::micro_xs[p->event_nuclide].fission) {
+            double nu = simulation::micro_xs[p->event_nuclide].nu_fission
+              / simulation::micro_xs[p->event_nuclide].fission;
             double dsig_s, dsig_a, dsig_f;
             std::tie(dsig_s, dsig_a, dsig_f)
               = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
@@ -592,9 +590,8 @@ score_track_derivative(const Particle* p, double distance)
       // phi is proportional to e^(-Sigma_tot * dist)
       // (1 / phi) * (d_phi / d_N) = - (d_Sigma_tot / d_N) * dist
       // (1 / phi) * (d_phi / d_N) = - sigma_tot * dist
-      //TODO: off-by-one
       deriv.flux_deriv -= distance
-        * simulation::micro_xs[deriv.diff_nuclide-1].total;
+        * simulation::micro_xs[deriv.diff_nuclide].total;
       break;
 
     case DIFF_TEMPERATURE:
@@ -635,17 +632,16 @@ void score_collision_derivative(const Particle* p)
       break;
 
     case DIFF_NUCLIDE_DENSITY:
-      //TODO: off-by-one throughout on diff_nuclide
       if (p->event_nuclide != deriv.diff_nuclide) continue;
       // Find the index in this material for the diff_nuclide.
       int i;
       for (i = 0; i < material.nuclide_.size(); ++i)
-        if (material.nuclide_[i] == deriv.diff_nuclide - 1) break;
+        if (material.nuclide_[i] == deriv.diff_nuclide) break;
       // Make sure we found the nuclide.
-      if (material.nuclide_[i] != deriv.diff_nuclide - 1) {
+      if (material.nuclide_[i] != deriv.diff_nuclide) {
         std::stringstream err_msg;
         err_msg << "Could not find nuclide "
-          << data::nuclides[deriv.diff_nuclide-1]->name_ << " in material "
+          << data::nuclides[deriv.diff_nuclide]->name_ << " in material "
           << material.id_ << " for tally derivative " << deriv.id;
         fatal_error(err_msg);
       }
@@ -661,8 +657,7 @@ void score_collision_derivative(const Particle* p)
       for (auto i_nuc : material.nuclide_) {
         const auto& nuc {*data::nuclides[i_nuc]};
         //TODO: off-by-one
-        if (i_nuc == p->event_nuclide - 1
-            && multipole_in_range(&nuc, p->last_E)) {
+        if (i_nuc == p->event_nuclide && multipole_in_range(&nuc, p->last_E)) {
           // phi is proportional to Sigma_s
           // (1 / phi) * (d_phi / d_T) = (d_Sigma_s / d_T) / Sigma_s
           // (1 / phi) * (d_phi / d_T) = (d_sigma_s / d_T) / sigma_s
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index a0c2883b00..96aa1abaa3 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -337,9 +337,6 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
 
     double score;
 
-    //TODO: off-by-one throughout on p->event_nuclide
-    //TODO: off-by-one throughout on p->material
-
     switch (score_bin) {
 
 
@@ -445,8 +442,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
       } else {
         // Get yield and apply to score
         auto m =
-          data::nuclides[p->event_nuclide-1]->reaction_index_[p->event_MT];
-        const auto& rxn {*data::nuclides[p->event_nuclide-1]->reactions_[m]};
+          data::nuclides[p->event_nuclide]->reaction_index_[p->event_MT];
+        const auto& rxn {*data::nuclides[p->event_nuclide]->reactions_[m]};
         score = p->last_wgt * flux * (*rxn.products_[0].yield_)(E);
       }
       break;
@@ -483,10 +480,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
           // No fission events occur if survival biasing is on -- need to
           // calculate fraction of absorptions that would have resulted in
           // fission
-          if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) {
+          if (simulation::micro_xs[p->event_nuclide].absorption > 0) {
             score = p->absorb_wgt
-              * simulation::micro_xs[p->event_nuclide-1].fission
-              / simulation::micro_xs[p->event_nuclide-1].absorption * flux;
+              * simulation::micro_xs[p->event_nuclide].fission
+              / simulation::micro_xs[p->event_nuclide].absorption * flux;
           } else {
             score = 0.;
           }
@@ -497,8 +494,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
           // weight entering the collision as the estimate for the fission
           // reaction rate
           score = p->last_wgt
-            * simulation::micro_xs[p->event_nuclide-1].fission
-            / simulation::micro_xs[p->event_nuclide-1].absorption * flux;
+            * simulation::micro_xs[p->event_nuclide].fission
+            / simulation::micro_xs[p->event_nuclide].absorption * flux;
         }
       } else {
         if (i_nuclide >= 0) {
@@ -525,10 +522,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
           // No fission events occur if survival biasing is on -- need to
           // calculate fraction of absorptions that would have resulted in
           // nu-fission
-          if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) {
+          if (simulation::micro_xs[p->event_nuclide].absorption > 0) {
             score = p->absorb_wgt
-              * simulation::micro_xs[p->event_nuclide-1].nu_fission
-              / simulation::micro_xs[p->event_nuclide-1].absorption * flux;
+              * simulation::micro_xs[p->event_nuclide].nu_fission
+              / simulation::micro_xs[p->event_nuclide].absorption * flux;
           } else {
             score = 0.;
           }
@@ -568,12 +565,12 @@ score_general_ce(const Particle* p, int i_tally, int start_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 (simulation::micro_xs[p->event_nuclide-1].absorption > 0) {
+          if (simulation::micro_xs[p->event_nuclide].absorption > 0) {
             score = p->absorb_wgt
-              * simulation::micro_xs[p->event_nuclide-1].fission
-              * data::nuclides[p->event_nuclide-1]
+              * simulation::micro_xs[p->event_nuclide].fission
+              * data::nuclides[p->event_nuclide]
               ->nu(E, ReactionProduct::EmissionMode::prompt)
-              / simulation::micro_xs[p->event_nuclide-1].absorption * flux;
+              / simulation::micro_xs[p->event_nuclide].absorption * flux;
           } else {
             score = 0.;
           }
@@ -630,8 +627,8 @@ score_general_ce(const Particle* p, int i_tally, int start_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 (simulation::micro_xs[p->event_nuclide-1].absorption > 0
-            && data::nuclides[p->event_nuclide-1]->fissionable_) {
+          if (simulation::micro_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
@@ -640,11 +637,11 @@ score_general_ce(const Particle* p, int i_tally, int start_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-1]
+                auto yield = data::nuclides[p->event_nuclide]
                   ->nu(E, ReactionProduct::EmissionMode::delayed, d);
                 score = p->absorb_wgt * yield
-                  * simulation::micro_xs[p->event_nuclide-1].fission
-                  / simulation::micro_xs[p->event_nuclide-1].absorption * flux;
+                  * simulation::micro_xs[p->event_nuclide].fission
+                  / simulation::micro_xs[p->event_nuclide].absorption * flux;
                 score_fission_delayed_dg(i_tally, d_bin+1, score,
                   score_index);
               }
@@ -654,10 +651,10 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
               // by multiplying the absorbed weight by the fraction of the
               // delayed-nu-fission xs to the absorption xs
               score = p->absorb_wgt
-                * simulation::micro_xs[p->event_nuclide-1].fission
-                * data::nuclides[p->event_nuclide-1]
+                * simulation::micro_xs[p->event_nuclide].fission
+                * data::nuclides[p->event_nuclide]
                 ->nu(E, ReactionProduct::EmissionMode::delayed)
-                / simulation::micro_xs[p->event_nuclide-1].absorption *flux;
+                / simulation::micro_xs[p->event_nuclide].absorption *flux;
             }
           }
         } else {
@@ -766,8 +763,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
           // 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-1]};
-          if (simulation::micro_xs[p->event_nuclide-1].absorption > 0
+          const auto& nuc {*data::nuclides[p->event_nuclide]};
+          if (simulation::micro_xs[p->event_nuclide].absorption > 0
             && nuc.fissionable_) {
             const auto& rxn {*nuc.fission_rx_[0]};
             if (tally.delayedgroup_filter_ != C_NONE) {
@@ -782,8 +779,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
                   = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d);
                 auto rate = rxn.products_[d].decay_rate_;
                 score = p->absorb_wgt * yield
-                  * simulation::micro_xs[p->event_nuclide-1].fission
-                  / simulation::micro_xs[p->event_nuclide-1].absorption
+                  * simulation::micro_xs[p->event_nuclide].fission
+                  / simulation::micro_xs[p->event_nuclide].absorption
                   * rate * flux;
                 score_fission_delayed_dg(i_tally, d_bin+1, score,
                   score_index);
@@ -805,8 +802,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
                   = nuc.nu(E, ReactionProduct::EmissionMode::delayed, d+1);
                 auto rate = rxn.products_[d+1].decay_rate_;
                 score += rate * p->absorb_wgt
-                  * simulation::micro_xs[p->event_nuclide-1].fission * yield
-                  / simulation::micro_xs[p->event_nuclide-1].absorption * flux;
+                  * simulation::micro_xs[p->event_nuclide].fission * yield
+                  / simulation::micro_xs[p->event_nuclide].absorption * flux;
               }
             }
           }
@@ -826,7 +823,7 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
             const auto& bank = simulation::fission_bank[i_bank];
             auto g = bank.delayed_group;
             if (g != 0) {
-              const auto& nuc {*data::nuclides[p->event_nuclide-1]};
+              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;
@@ -953,13 +950,13 @@ score_general_ce(const Particle* p, int i_tally, int start_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-1]};
-          if (simulation::micro_xs[p->event_nuclide-1].absorption > 0
+          const auto& nuc {*data::nuclides[p->event_nuclide]};
+          if (simulation::micro_xs[p->event_nuclide].absorption > 0
             && nuc.fissionable_) {
             const auto& rxn {*nuc.fission_rx_[0]};
             score = p->absorb_wgt * rxn.q_value_
-              * simulation::micro_xs[p->event_nuclide-1].fission
-              / simulation::micro_xs[p->event_nuclide-1].absorption * flux;
+              * simulation::micro_xs[p->event_nuclide].fission
+              / simulation::micro_xs[p->event_nuclide].absorption * flux;
           }
         } else {
           // Skip any non-absorption events
@@ -967,13 +964,13 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
           // 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-1]};
-          if (simulation::micro_xs[p->event_nuclide-1].absorption > 0
+          const auto& nuc {*data::nuclides[p->event_nuclide]};
+          if (simulation::micro_xs[p->event_nuclide].absorption > 0
             && nuc.fissionable_) {
             const auto& rxn {*nuc.fission_rx_[0]};
             score = p->last_wgt * rxn.q_value_
-              * simulation::micro_xs[p->event_nuclide-1].fission
-              / simulation::micro_xs[p->event_nuclide-1].absorption * flux;
+              * simulation::micro_xs[p->event_nuclide].fission
+              / simulation::micro_xs[p->event_nuclide].absorption * flux;
           }
         }
       } else {
@@ -1046,8 +1043,8 @@ score_general_ce(const Particle* p, int i_tally, int start_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-1]};
-          if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) {
+          const auto& nuc {*data::nuclides[p->event_nuclide]};
+          if (simulation::micro_xs[p->event_nuclide].absorption > 0) {
             double q_value = 0.;
             if (score_bin == SCORE_FISS_Q_PROMPT) {
               if (nuc.fission_q_prompt_)
@@ -1057,8 +1054,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
                 q_value = (*nuc.fission_q_recov_)(p->last_E);
             }
             score = p->absorb_wgt * q_value
-              * simulation::micro_xs[p->event_nuclide-1].fission
-              / simulation::micro_xs[p->event_nuclide-1].absorption * flux;
+              * simulation::micro_xs[p->event_nuclide].fission
+              / simulation::micro_xs[p->event_nuclide].absorption * flux;
           }
         } else {
           // Skip any non-absorption events
@@ -1066,8 +1063,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
           // 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-1]};
-          if (simulation::micro_xs[p->event_nuclide-1].absorption > 0) {
+          const auto& nuc {*data::nuclides[p->event_nuclide]};
+          if (simulation::micro_xs[p->event_nuclide].absorption > 0) {
             double q_value = 0.;
             if (score_bin == SCORE_FISS_Q_PROMPT) {
               if (nuc.fission_q_prompt_)
@@ -1077,8 +1074,8 @@ score_general_ce(const Particle* p, int i_tally, int start_index,
                 q_value = (*nuc.fission_q_recov_)(p->last_E);
             }
             score = p->last_wgt * q_value
-              * simulation::micro_xs[p->event_nuclide-1].fission
-              / simulation::micro_xs[p->event_nuclide-1].absorption * flux;
+              * simulation::micro_xs[p->event_nuclide].fission
+              / simulation::micro_xs[p->event_nuclide].absorption * flux;
           }
         }
       } else {
@@ -1992,8 +1989,7 @@ void score_analog_tally_ce(const Particle* p)
           // the event nuclide or the total material.  Note that the i_nuclide
           // and flux arguments for score_general are not used for analog
           // tallies.
-          //TODO: off-by-one
-          if (i_nuclide == p->event_nuclide-1 || i_nuclide == -1)
+          if (i_nuclide == p->event_nuclide || i_nuclide == -1)
             score_general_ce(p, i_tally, i*tally.scores_.size(), filter_index,
               -1, -1., filter_weight);
         }