diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h
index 9d680434ce..145fbcf058 100644
--- a/include/openmc/mgxs.h
+++ b/include/openmc/mgxs.h
@@ -134,7 +134,8 @@ class Mgxs {
     //! @param dg delayed group index; use nullptr if irrelevant.
     //! @return Requested cross section value.
     double
-    get_xs(int xstype, int gin, int* gout, double* mu, int* dg);
+    get_xs(int xstype, int gin, const int* gout, const double* mu,
+      const int* dg);
 
     //! \brief Samples the fission neutron energy and if prompt or delayed.
     //!
diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h
index 413cdad621..5f5f219e67 100644
--- a/include/openmc/mgxs_interface.h
+++ b/include/openmc/mgxs_interface.h
@@ -52,14 +52,16 @@ calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
      double& total_xs, double& abs_xs, double& nu_fiss_xs);
 
 double
-get_nuclide_xs(int index, int xstype, int gin, int* gout, double* mu, int* dg);
+get_nuclide_xs(int index, int xstype, int gin, const int* gout,
+  const double* mu, const int* dg);
 
 inline double
 get_nuclide_xs(int index, int xstype, int gin)
 {return get_nuclide_xs(index, xstype, gin, nullptr, nullptr, nullptr);}
 
 double
-get_macro_xs(int index, int xstype, int gin, int* gout, double* mu, int* dg);
+get_macro_xs(int index, int xstype, int gin, const int* gout,
+  const double* mu, const int* dg);
 
 inline double
 get_macro_xs(int index, int xstype, int gin)
diff --git a/include/openmc/tallies/derivative.h b/include/openmc/tallies/derivative.h
index da9fa60a6e..deee9d7aa9 100644
--- a/include/openmc/tallies/derivative.h
+++ b/include/openmc/tallies/derivative.h
@@ -14,7 +14,7 @@
 
 namespace openmc {
 
-extern "C" struct TallyDerivative {
+struct TallyDerivative {
   int id;  //!< User-defined identifier
   int variable;  //!< Independent variable (like temperature)
   int diff_material;  //!< Material this derivative is applied to
@@ -32,8 +32,8 @@ extern "C" struct TallyDerivative {
 //! Scale the given score by its logarithmic derivative
 
 void
-apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
-  double atom_density, int score_bin, double* score);
+apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
+  double atom_density, int score_bin, double& score);
 
 } // namespace openmc
 
@@ -46,18 +46,19 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
 extern template class std::vector<openmc::TallyDerivative>;
 
 namespace openmc {
-  namespace model {
-    extern std::vector<TallyDerivative> tally_derivs;
-    #pragma omp threadprivate(tally_derivs)
 
-    extern std::unordered_map<int, int> tally_deriv_map;
-  }
+namespace model {
+extern std::vector<TallyDerivative> tally_derivs;
+#pragma omp threadprivate(tally_derivs)
+extern std::unordered_map<int, int> tally_deriv_map;
+} // namespace model
 
-  // Independent variables
-  //TODO: convert to enum
-  constexpr int DIFF_DENSITY {1};
-  constexpr int DIFF_NUCLIDE_DENSITY {2};
-  constexpr int DIFF_TEMPERATURE {3};
-}
+// Independent variables
+//TODO: convert to enum
+constexpr int DIFF_DENSITY {1};
+constexpr int DIFF_NUCLIDE_DENSITY {2};
+constexpr int DIFF_TEMPERATURE {3};
+
+} // namespace openmc
 
 #endif // OPENMC_TALLIES_DERIVATIVE_H
diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h
index 0bb2127c0e..464ec9f3d9 100644
--- a/include/openmc/tallies/tally.h
+++ b/include/openmc/tallies/tally.h
@@ -104,15 +104,17 @@ private:
 extern "C" double total_weight;
 
 namespace model {
-  extern std::vector<std::unique_ptr<Tally>> tallies;
 
-  extern std::vector<int> active_tallies;
-  extern std::vector<int> active_analog_tallies;
-  extern std::vector<int> active_tracklength_tallies;
-  extern std::vector<int> active_collision_tallies;
-  extern std::vector<int> active_meshsurf_tallies;
-  extern std::vector<int> active_surface_tallies;
-}
+extern std::vector<std::unique_ptr<Tally>> tallies;
+
+extern std::vector<int> active_tallies;
+extern std::vector<int> active_analog_tallies;
+extern std::vector<int> active_tracklength_tallies;
+extern std::vector<int> active_collision_tallies;
+extern std::vector<int> active_meshsurf_tallies;
+extern std::vector<int> active_surface_tallies;
+
+} // namespace model
 
 // Threadprivate variables
 extern "C" double global_tally_absorption;
diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h
index ea415c514b..8aa659325a 100644
--- a/include/openmc/tallies/tally_scoring.h
+++ b/include/openmc/tallies/tally_scoring.h
@@ -22,7 +22,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, const Particle* p);
 
   //! Construct an iterator over all filter bin combinations.
   //
diff --git a/include/openmc/tallies/trigger.h b/include/openmc/tallies/trigger.h
index e670297ed1..a11759114a 100644
--- a/include/openmc/tallies/trigger.h
+++ b/include/openmc/tallies/trigger.h
@@ -17,8 +17,7 @@ enum class TriggerMetric {
 
 //! Stops the simulation early if a desired tally uncertainty is reached.
 
-struct Trigger
-{
+struct Trigger {
   TriggerMetric metric;  //!< The type of uncertainty (e.g. std dev) measured
   double threshold;  //!< Uncertainty value below which trigger is satisfied
   int score_index;  //!< Index of the relevant score in the tally's arrays
diff --git a/src/mgxs.cpp b/src/mgxs.cpp
index 08c4f78941..5b37d8ad54 100644
--- a/src/mgxs.cpp
+++ b/src/mgxs.cpp
@@ -424,7 +424,8 @@ Mgxs::combine(const std::vector<Mgxs*>& micros, const std::vector<double>& scala
 //==============================================================================
 
 double
-Mgxs::get_xs(int xstype, int gin, int* gout, double* mu, int* dg)
+Mgxs::get_xs(int xstype, int gin, const int* gout, const double* mu,
+  const int* dg)
 {
   // This method assumes that the temperature and angle indices are set
 #ifdef _OPENMP
diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp
index 546e921005..f239dcedcb 100644
--- a/src/mgxs_interface.cpp
+++ b/src/mgxs_interface.cpp
@@ -231,12 +231,13 @@ calculate_xs_c(int i_mat, int gin, double sqrtkT, const double uvw[3],
 //==============================================================================
 
 double
-get_nuclide_xs(int index, int xstype, int gin, int* gout, double* mu, int* dg)
+get_nuclide_xs(int index, int xstype, int gin, const int* gout,
+  const double* mu, const int* dg)
 {
   int gout_c;
-  int* gout_c_p;
+  const int* gout_c_p;
   int dg_c;
-  int* dg_c_p;
+  const int* dg_c_p;
   if (gout != nullptr) {
     gout_c = *gout - 1;
     gout_c_p = &gout_c;
@@ -255,12 +256,13 @@ get_nuclide_xs(int index, int xstype, int gin, int* gout, double* mu, int* dg)
 //==============================================================================
 
 double
-get_macro_xs(int index, int xstype, int gin, int* gout, double* mu, int* dg)
+get_macro_xs(int index, int xstype, int gin, const int* gout,
+  const double* mu, const int* dg)
 {
   int gout_c;
-  int* gout_c_p;
+  const int* gout_c_p;
   int dg_c;
-  int* dg_c_p;
+  const int* dg_c_p;
   if (gout != nullptr) {
     gout_c = *gout - 1;
     gout_c_p = &gout_c;
diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp
index df65b0751a..b82db28e64 100644
--- a/src/tallies/derivative.cpp
+++ b/src/tallies/derivative.cpp
@@ -1,9 +1,10 @@
+#include "openmc/tallies/derivative.h"
+
 #include "openmc/error.h"
 #include "openmc/material.h"
 #include "openmc/nuclide.h"
 #include "openmc/settings.h"
 #include "openmc/tallies/tally.h"
-#include "openmc/tallies/derivative.h"
 #include "openmc/xml_interface.h"
 
 #include <sstream>
@@ -107,12 +108,12 @@ read_tally_derivatives(pugi::xml_node* node)
 }
 
 void
-apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
-  double atom_density, int score_bin, double* score)
+apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
+  double atom_density, int score_bin, double& score)
 {
   const Tally& tally {*model::tallies[i_tally]};
 
-  if (*score == 0) return;
+  if (score == 0.0) return;
 
   // If our score was previously c then the new score is
   // c * (1/f * d_f/d_p + 1/c * d_c/d_p)
@@ -124,16 +125,16 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
 
   // Handle special cases where we know that d_c/d_p must be zero.
   if (score_bin == SCORE_FLUX) {
-    *score *= flux_deriv;
+    score *= flux_deriv;
     return;
   } else if (p->material == MATERIAL_VOID) {
-    *score *= flux_deriv;
+    score *= flux_deriv;
     return;
   }
   //TODO: off-by-one
   const Material& material {*model::materials[p->material-1]};
   if (material.id_ != deriv.diff_material) {
-    *score *= flux_deriv;
+    score *= flux_deriv;
     return;
   }
 
@@ -159,7 +160,7 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
       case SCORE_ABSORPTION:
       case SCORE_FISSION:
       case SCORE_NU_FISSION:
-        *score *= flux_deriv + 1. / material.density_gpcc_;
+        score *= flux_deriv + 1. / material.density_gpcc_;
         break;
 
       default:
@@ -193,7 +194,7 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
 
     case ESTIMATOR_ANALOG:
       if (p->event_nuclide != deriv.diff_nuclide) {
-        *score *= flux_deriv;
+        score *= flux_deriv;
         return;
       }
 
@@ -209,7 +210,7 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
           int i;
           for (i = 0; i < material.nuclide_.size(); ++i)
             if (material.nuclide_[i] == deriv.diff_nuclide - 1) break;
-          *score *= flux_deriv + 1. / material.atom_density_(i);
+          score *= flux_deriv + 1. / material.atom_density_(i);
         }
         break;
 
@@ -223,69 +224,69 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
       switch (score_bin) {
 
       case SCORE_TOTAL:
-        if (i_nuclide == -1 && simulation::material_xs.total) {
-          *score *= flux_deriv
+        if (i_nuclide == -1 && simulation::material_xs.total > 0.0) {
+          score *= flux_deriv
             + simulation::micro_xs[deriv.diff_nuclide-1].total
             / simulation::material_xs.total;
         } else if (i_nuclide == deriv.diff_nuclide-1
                    && simulation::micro_xs[i_nuclide].total) {
-          *score *= flux_deriv + 1. / atom_density;
+          score *= flux_deriv + 1. / atom_density;
         } else {
-          *score *= flux_deriv;
+          score *= flux_deriv;
         }
         break;
 
       case SCORE_SCATTER:
         if (i_nuclide == -1 && (simulation::material_xs.total
-                                - simulation::material_xs.absorption)) {
-          *score *= flux_deriv
+                                - 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::material_xs.total
             - simulation::material_xs.absorption);
         } else if (i_nuclide == deriv.diff_nuclide-1) {
-          *score *= flux_deriv + 1. / atom_density;
+          score *= flux_deriv + 1. / atom_density;
         } else {
-          *score *= flux_deriv;
+          score *= flux_deriv;
         }
         break;
 
       case SCORE_ABSORPTION:
-        if (i_nuclide == -1 && simulation::material_xs.absorption) {
-          *score *= flux_deriv
+        if (i_nuclide == -1 && simulation::material_xs.absorption > 0.0) {
+          score *= flux_deriv
             + simulation::micro_xs[deriv.diff_nuclide-1].absorption
             / simulation::material_xs.absorption;
         } else if (i_nuclide == deriv.diff_nuclide-1
                    && simulation::micro_xs[i_nuclide].absorption) {
-          *score *= flux_deriv + 1. / atom_density;
+          score *= flux_deriv + 1. / atom_density;
         } else {
-          *score *= flux_deriv;
+          score *= flux_deriv;
         }
         break;
 
       case SCORE_FISSION:
-        if (i_nuclide == -1 && simulation::material_xs.fission) {
-          *score *= flux_deriv
+        if (i_nuclide == -1 && simulation::material_xs.fission > 0.0) {
+          score *= flux_deriv
             + simulation::micro_xs[deriv.diff_nuclide-1].fission
             / simulation::material_xs.fission;
         } else if (i_nuclide == deriv.diff_nuclide-1
                    && simulation::micro_xs[i_nuclide].fission) {
-          *score *= flux_deriv + 1. / atom_density;
+          score *= flux_deriv + 1. / atom_density;
         } else {
-          *score *= flux_deriv;
+          score *= flux_deriv;
         }
         break;
 
       case SCORE_NU_FISSION:
-        if (i_nuclide == -1 && simulation::material_xs.nu_fission) {
-          *score *= flux_deriv
+        if (i_nuclide == -1 && simulation::material_xs.nu_fission > 0.0) {
+          score *= flux_deriv
             + simulation::micro_xs[deriv.diff_nuclide-1].nu_fission
             / simulation::material_xs.nu_fission;
         } else if (i_nuclide == deriv.diff_nuclide-1
                    && simulation::micro_xs[i_nuclide].nu_fission) {
-          *score *= flux_deriv + 1. / atom_density;
+          score *= flux_deriv + 1. / atom_density;
         } else {
-          *score *= flux_deriv;
+          score *= flux_deriv;
         }
         break;
 
@@ -326,7 +327,7 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
 
         const auto& nuc {*data::nuclides[p->event_nuclide-1]};
         if (!multipole_in_range(&nuc, p->last_E)) {
-          *score *= flux_deriv;
+          score *= flux_deriv;
           break;
         }
 
@@ -337,10 +338,10 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
             double dsig_s, dsig_a, dsig_f;
             std::tie(dsig_s, dsig_a, dsig_f)
               = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
-            *score *= flux_deriv + (dsig_s + dsig_a) * material.atom_density_(i)
+            score *= flux_deriv + (dsig_s + dsig_a) * material.atom_density_(i)
               / simulation::material_xs.total;
           } else {
-            *score *= flux_deriv;
+            score *= flux_deriv;
           }
           break;
 
@@ -350,11 +351,11 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
             double dsig_s, dsig_a, dsig_f;
             std::tie(dsig_s, dsig_a, dsig_f)
               = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
-            *score *= flux_deriv + dsig_s * material.atom_density_(i)
+            score *= flux_deriv + dsig_s * material.atom_density_(i)
               / (simulation::material_xs.total
               - simulation::material_xs.absorption);
           } else {
-            *score *= flux_deriv;
+            score *= flux_deriv;
           }
           break;
 
@@ -363,10 +364,10 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
             double dsig_s, dsig_a, dsig_f;
             std::tie(dsig_s, dsig_a, dsig_f)
               = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
-            *score *= flux_deriv + dsig_a * material.atom_density_(i)
+            score *= flux_deriv + dsig_a * material.atom_density_(i)
               / simulation::material_xs.absorption;
           } else {
-            *score *= flux_deriv;
+            score *= flux_deriv;
           }
           break;
 
@@ -375,10 +376,10 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
             double dsig_s, dsig_a, dsig_f;
             std::tie(dsig_s, dsig_a, dsig_f)
               = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
-            *score *= flux_deriv + dsig_f * material.atom_density_(i)
+            score *= flux_deriv + dsig_f * material.atom_density_(i)
               / simulation::material_xs.fission;
           } else {
-            *score *= flux_deriv;
+            score *= flux_deriv;
           }
           break;
 
@@ -389,10 +390,10 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
             double dsig_s, dsig_a, dsig_f;
             std::tie(dsig_s, dsig_a, dsig_f)
               = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
-            *score *= flux_deriv + nu * dsig_f * material.atom_density_(i)
+            score *= flux_deriv + nu * dsig_f * material.atom_density_(i)
               / simulation::material_xs.nu_fission;
           } else {
-            *score *= flux_deriv;
+            score *= flux_deriv;
           }
           break;
 
@@ -405,9 +406,9 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
 
     case ESTIMATOR_COLLISION:
       if (i_nuclide != -1) {
-        const auto& nuc {*data::nuclides[i_nuclide]};
-        if (!multipole_in_range(&nuc, p->last_E)) {
-          *score *= flux_deriv;
+        const auto& nuc {data::nuclides[i_nuclide]};
+        if (!multipole_in_range(nuc.get(), p->last_E)) {
+          score *= flux_deriv;
           return;
         }
       }
@@ -415,7 +416,7 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
       switch (score_bin) {
 
       case SCORE_TOTAL:
-        if (i_nuclide == -1 && simulation::material_xs.total) {
+        if (i_nuclide == -1 && simulation::material_xs.total > 0.0) {
           double cum_dsig = 0;
           for (auto i = 0; i < material.nuclide_.size(); ++i) {
             auto i_nuc = material.nuclide_[i];
@@ -428,16 +429,16 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
               cum_dsig += (dsig_s + dsig_a) * material.atom_density_(i);
             }
           }
-          *score *= flux_deriv + cum_dsig / simulation::material_xs.total;
+          score *= flux_deriv + cum_dsig / simulation::material_xs.total;
         } else if (simulation::micro_xs[i_nuclide].total) {
           const auto& nuc {*data::nuclides[i_nuclide]};
           double dsig_s, dsig_a, dsig_f;
           std::tie(dsig_s, dsig_a, dsig_f)
             = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
-          *score *= flux_deriv
+          score *= flux_deriv
             + (dsig_s + dsig_a) / simulation::micro_xs[i_nuclide].total;
         } else {
-          *score *= flux_deriv;
+          score *= flux_deriv;
         }
         break;
 
@@ -457,7 +458,7 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
               cum_dsig += dsig_s * material.atom_density_(i);
             }
           }
-          *score *= flux_deriv + cum_dsig / (simulation::material_xs.total
+          score *= flux_deriv + cum_dsig / (simulation::material_xs.total
             - simulation::material_xs.absorption);
         } else if (simulation::micro_xs[i_nuclide].total
                    - simulation::micro_xs[i_nuclide].absorption) {
@@ -465,15 +466,15 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
           double dsig_s, dsig_a, dsig_f;
           std::tie(dsig_s, dsig_a, dsig_f)
             = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
-          *score *= flux_deriv + dsig_s / (simulation::micro_xs[i_nuclide].total
+          score *= flux_deriv + dsig_s / (simulation::micro_xs[i_nuclide].total
             - simulation::micro_xs[i_nuclide].absorption);
         } else {
-          *score *= flux_deriv;
+          score *= flux_deriv;
         }
         break;
 
       case SCORE_ABSORPTION:
-        if (i_nuclide == -1 && simulation::material_xs.absorption) {
+        if (i_nuclide == -1 && simulation::material_xs.absorption > 0.0) {
           double cum_dsig = 0;
           for (auto i = 0; i < material.nuclide_.size(); ++i) {
             auto i_nuc = material.nuclide_[i];
@@ -486,21 +487,21 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
               cum_dsig += dsig_a * material.atom_density_(i);
             }
           }
-          *score *= flux_deriv + cum_dsig / simulation::material_xs.absorption;
+          score *= flux_deriv + cum_dsig / simulation::material_xs.absorption;
         } else if (simulation::micro_xs[i_nuclide].absorption) {
           const auto& nuc {*data::nuclides[i_nuclide]};
           double dsig_s, dsig_a, dsig_f;
           std::tie(dsig_s, dsig_a, dsig_f)
             = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
-          *score *= flux_deriv
+          score *= flux_deriv
             + dsig_a / simulation::micro_xs[i_nuclide].absorption;
         } else {
-          *score *= flux_deriv;
+          score *= flux_deriv;
         }
         break;
 
       case SCORE_FISSION:
-        if (i_nuclide == -1 && simulation::material_xs.fission) {
+        if (i_nuclide == -1 && simulation::material_xs.fission > 0.0) {
           double cum_dsig = 0;
           for (auto i = 0; i < material.nuclide_.size(); ++i) {
             auto i_nuc = material.nuclide_[i];
@@ -513,21 +514,21 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
               cum_dsig += dsig_f * material.atom_density_(i);
             }
           }
-          *score *= flux_deriv + cum_dsig / simulation::material_xs.fission;
+          score *= flux_deriv + cum_dsig / simulation::material_xs.fission;
         } else if (simulation::micro_xs[i_nuclide].fission) {
           const auto& nuc {*data::nuclides[i_nuclide]};
           double dsig_s, dsig_a, dsig_f;
           std::tie(dsig_s, dsig_a, dsig_f)
             = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
-          *score *= flux_deriv
+          score *= flux_deriv
             + dsig_f / simulation::micro_xs[i_nuclide].fission;
         } else {
-          *score *= flux_deriv;
+          score *= flux_deriv;
         }
         break;
 
       case SCORE_NU_FISSION:
-        if (i_nuclide == -1 && simulation::material_xs.nu_fission) {
+        if (i_nuclide == -1 && simulation::material_xs.nu_fission > 0.0) {
           double cum_dsig = 0;
           for (auto i = 0; i < material.nuclide_.size(); ++i) {
             auto i_nuc = material.nuclide_[i];
@@ -542,16 +543,16 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
               cum_dsig += nu * dsig_f * material.atom_density_(i);
             }
           }
-          *score *= flux_deriv + cum_dsig / simulation::material_xs.nu_fission;
+          score *= flux_deriv + cum_dsig / simulation::material_xs.nu_fission;
         } else if (simulation::micro_xs[i_nuclide].fission) {
           const auto& nuc {*data::nuclides[i_nuclide]};
           double dsig_s, dsig_a, dsig_f;
           std::tie(dsig_s, dsig_a, dsig_f)
             = nuc.multipole_->evaluate_deriv(p->last_E, p->sqrtkT);
-          *score *= flux_deriv
+          score *= flux_deriv
             + dsig_f / simulation::micro_xs[i_nuclide].fission;
         } else {
-          *score *= flux_deriv;
+          score *= flux_deriv;
         }
         break;
 
@@ -571,7 +572,7 @@ apply_derivative_to_score(Particle* p, int i_tally, int i_nuclide,
 //! Adjust diff tally flux derivatives for a particle tracking event.
 
 extern "C" void
-score_track_derivative(Particle* p, double distance)
+score_track_derivative(const Particle* p, double distance)
 {
   // A void material cannot be perturbed so it will not affect flux derivatives.
   if (p->material == MATERIAL_VOID) return;
@@ -627,7 +628,7 @@ score_track_derivative(Particle* p, double distance)
 //! further tallies are scored.
 
 extern "C" void
-score_collision_derivative(Particle* p)
+score_collision_derivative(const Particle* p)
 {
   // A void material cannot be perturbed so it will not affect flux derivatives.
   if (p->material == MATERIAL_VOID) return;
@@ -713,7 +714,7 @@ extern "C" int n_tally_derivs() {return model::tally_derivs.size();}
 extern "C" TallyDerivative*
 tally_deriv_c(int i)
 {
-  return &(model::tally_derivs[i]);
+  return &model::tally_derivs[i];
 }
 
 }// namespace openmc
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index e9bee044df..4bfb15effb 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -637,7 +637,6 @@ free_memory_tally_c()
 {
   #pragma omp parallel
   {
-    simulation::filter_matches.clear();
     model::tally_derivs.clear();
   }
 
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index 77bfcd2738..8c509a0aae 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -23,7 +23,7 @@ namespace openmc {
 // FilterBinIter implementation
 //==============================================================================
 
-FilterBinIter::FilterBinIter(const Tally& tally, Particle* p)
+FilterBinIter::FilterBinIter(const Tally& tally, const Particle* p)
   : tally_{tally}
 {
   // Find all valid bins in each relevant filter if they have not already been
@@ -179,7 +179,7 @@ score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index)
 //! neutrons produced with different energies.
 
 void
-score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin)
+score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
 {
   const Tally& tally {*model::tallies[i_tally]};
   auto results = tally_results(i_tally);
@@ -210,7 +210,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_) {
 
@@ -325,8 +325,8 @@ 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,
-  int i_nuclide, double atom_density, double flux)
+score_general_ce(const Particle* p, int i_tally, int start_index,
+  int filter_index, int i_nuclide, double atom_density, double flux)
 {
   const Tally& tally {*model::tallies[i_tally]};
   auto results = tally_results(i_tally);
@@ -1218,7 +1218,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, int filter_index,
     // Add derivative information on score for differnetial tallies.
     if (tally.deriv_ != C_NONE)
       apply_derivative_to_score(p, i_tally, i_nuclide, atom_density, score_bin,
-        &score);
+        score);
 
     // Update tally results
     #pragma omp atomic
@@ -1232,8 +1232,8 @@ 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,
-  int i_nuclide, double atom_density, double flux)
+score_general_mg(const Particle* p, int i_tally, int start_index,
+  int filter_index, int i_nuclide, double atom_density, double flux)
 {
   const Tally& tally {*model::tallies[i_tally]};
   auto results = tally_results(i_tally);
@@ -1241,7 +1241,7 @@ score_general_mg(Particle* p, int i_tally, int start_index, int filter_index,
   //TODO: off-by-one throughout on p->material
 
   // Set the direction and group to use with get_xs
-  double* p_uvw;
+  const double* p_uvw;
   int p_g;
   if (tally.estimator_ == ESTIMATOR_ANALOG
     || tally.estimator_ == ESTIMATOR_COLLISION) {
@@ -1937,7 +1937,8 @@ 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, int filter_index)
+score_all_nuclides(const Particle* p, int i_tally, double flux,
+  int filter_index)
 {
   const Tally& tally {*model::tallies[i_tally]};
   const Material& material {*model::materials[p->material-1]};
@@ -1976,7 +1977,7 @@ score_all_nuclides(Particle* p, int i_tally, double flux, int filter_index)
 //! Analog tallies ar etriggered at every collision, not every event.
 
 extern "C" void
-score_analog_tally_ce(Particle* p)
+score_analog_tally_ce(const Particle* p)
 {
   for (auto i_tally : model::active_analog_tallies) {
     const Tally& tally {*model::tallies[i_tally]};
@@ -2040,7 +2041,7 @@ score_analog_tally_ce(Particle* p)
 //! Analog tallies ar etriggered at every collision, not every event.
 
 extern "C" void
-score_analog_tally_mg(Particle* p)
+score_analog_tally_mg(const Particle* p)
 {
   for (auto i_tally : model::active_analog_tallies) {
     const Tally& tally {*model::tallies[i_tally]};
@@ -2096,7 +2097,7 @@ score_analog_tally_mg(Particle* p)
 //! information.
 
 extern "C" void
-score_tracklength_tally(Particle* p, double distance)
+score_tracklength_tally(const Particle* p, double distance)
 {
   // Determine the tracklength estimate of the flux
   double flux = p->wgt * distance;
@@ -2172,7 +2173,7 @@ score_tracklength_tally(Particle* p, double distance)
 //! since collisions do not occur in voids.
 
 extern "C" void
-score_collision_tally(Particle* p)
+score_collision_tally(const Particle* p)
 {
   // Determine the collision estimate of the flux
   double flux;
@@ -2243,7 +2244,7 @@ score_collision_tally(Particle* p)
 //! Score surface or mesh-surface tallies for particle currents.
 
 static void
-score_surface_tally_inner(Particle* p, const std::vector<int>& tallies)
+score_surface_tally_inner(const Particle* p, const std::vector<int>& tallies)
 {
   // No collision, so no weight change when survival biasing
   double flux = p->wgt;
@@ -2291,7 +2292,7 @@ score_surface_tally_inner(Particle* p, const std::vector<int>& tallies)
 //! Score mesh-surface tallies for particle currents.
 
 extern "C" void
-score_meshsurface_tally(Particle* p)
+score_meshsurface_tally(const Particle* p)
 {
   score_surface_tally_inner(p, model::active_meshsurf_tallies);
 }
@@ -2299,7 +2300,7 @@ score_meshsurface_tally(Particle* p)
 //! Score surface tallies for particle currents.
 
 extern "C" void
-score_surface_tally(Particle* p)
+score_surface_tally(const Particle* p)
 {
   score_surface_tally_inner(p, model::active_surface_tallies);
 }
diff --git a/src/tallies/trigger.cpp b/src/tallies/trigger.cpp
index 6cd99f1330..a53cef817e 100644
--- a/src/tallies/trigger.cpp
+++ b/src/tallies/trigger.cpp
@@ -26,7 +26,7 @@ namespace settings {
 // Non-member functions
 //==============================================================================
 
-static std::pair<double, double>
+std::pair<double, double>
 get_tally_uncertainty(int i_tally, int score_index, int filter_index)
 {
   int n;
@@ -87,6 +87,7 @@ check_tally_triggers(double& ratio, int& tally_id, int& score)
               break;
             case TriggerMetric::relative_error:
               uncertainty = rel_err;
+              break;
           }
 
           // Compute the uncertainty / threshold ratio.