diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index 3df2f40d5..e5797872e 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -25,10 +25,11 @@ namespace openmc {
 // Constants
 //==============================================================================
 
-// TODO: Convert to enum
-constexpr int FILL_MATERIAL {1};
-constexpr int FILL_UNIVERSE {2};
-constexpr int FILL_LATTICE {3};
+enum class Fill {
+  MATERIAL,
+  UNIVERSE,
+  LATTICE
+};
 
 // TODO: Convert to enum
 constexpr int32_t OP_LEFT_PAREN   {std::numeric_limits<int32_t>::max()};
@@ -149,7 +150,7 @@ public:
 
   int32_t id_;                //!< Unique ID
   std::string name_;          //!< User-defined name
-  int type_;                  //!< Material, universe, or lattice
+  Fill type_;                  //!< Material, universe, or lattice
   int32_t universe_;          //!< Universe # this cell is in
   int32_t fill_;              //!< Universe # filling this cell
   int32_t n_instances_{0};    //!< Number of instances of this cell
diff --git a/include/openmc/constants.h b/include/openmc/constants.h
index d55744117..28700bb9b 100644
--- a/include/openmc/constants.h
+++ b/include/openmc/constants.h
@@ -113,192 +113,170 @@ constexpr int NUCLIDE_NONE  {-1};
 // ============================================================================
 // CROSS SECTION RELATED CONSTANTS
 
-// Angular distribution type
-// TODO: Convert to enum
-constexpr int ANGLE_ISOTROPIC {1};
-constexpr int ANGLE_32_EQUI {2};
-constexpr int ANGLE_TABULAR {3};
-constexpr int ANGLE_LEGENDRE {4};
-constexpr int ANGLE_HISTOGRAM {5};
-
 // Temperature treatment method
-// TODO: Convert to enum?
-constexpr int TEMPERATURE_NEAREST {1};
-constexpr int TEMPERATURE_INTERPOLATION {2};
+enum class TemperatureInterpolationType {
+ NEAREST,
+ INTERPOLATION
+};
 
 // Reaction types
-// TODO: Convert to enum
-constexpr int REACTION_NONE {0};
-constexpr int TOTAL_XS {1};
-constexpr int ELASTIC {2};
-constexpr int N_NONELASTIC {3};
-constexpr int N_LEVEL {4};
-constexpr int MISC    {5};
-constexpr int N_2ND   {11};
-constexpr int N_2N    {16};
-constexpr int N_3N    {17};
-constexpr int N_FISSION {18};
-constexpr int N_F     {19};
-constexpr int N_NF    {20};
-constexpr int N_2NF   {21};
-constexpr int N_NA    {22};
-constexpr int N_N3A   {23};
-constexpr int N_2NA   {24};
-constexpr int N_3NA   {25};
-constexpr int N_NP    {28};
-constexpr int N_N2A   {29};
-constexpr int N_2N2A  {30};
-constexpr int N_ND    {32};
-constexpr int N_NT    {33};
-constexpr int N_N3HE  {34};
-constexpr int N_ND2A  {35};
-constexpr int N_NT2A  {36};
-constexpr int N_4N    {37};
-constexpr int N_3NF   {38};
-constexpr int N_2NP   {41};
-constexpr int N_3NP   {42};
-constexpr int N_N2P   {44};
-constexpr int N_NPA   {45};
-constexpr int N_N1    {51};
-constexpr int N_N40   {90};
-constexpr int N_NC    {91};
-constexpr int N_DISAPPEAR {101};
-constexpr int N_GAMMA {102};
-constexpr int N_P     {103};
-constexpr int N_D     {104};
-constexpr int N_T     {105};
-constexpr int N_3HE   {106};
-constexpr int N_A     {107};
-constexpr int N_2A    {108};
-constexpr int N_3A    {109};
-constexpr int N_2P    {111};
-constexpr int N_PA    {112};
-constexpr int N_T2A   {113};
-constexpr int N_D2A   {114};
-constexpr int N_PD    {115};
-constexpr int N_PT    {116};
-constexpr int N_DA    {117};
-constexpr int N_5N    {152};
-constexpr int N_6N    {153};
-constexpr int N_2NT   {154};
-constexpr int N_TA    {155};
-constexpr int N_4NP   {156};
-constexpr int N_3ND   {157};
-constexpr int N_NDA   {158};
-constexpr int N_2NPA  {159};
-constexpr int N_7N    {160};
-constexpr int N_8N    {161};
-constexpr int N_5NP   {162};
-constexpr int N_6NP   {163};
-constexpr int N_7NP   {164};
-constexpr int N_4NA   {165};
-constexpr int N_5NA   {166};
-constexpr int N_6NA   {167};
-constexpr int N_7NA   {168};
-constexpr int N_4ND   {169};
-constexpr int N_5ND   {170};
-constexpr int N_6ND   {171};
-constexpr int N_3NT   {172};
-constexpr int N_4NT   {173};
-constexpr int N_5NT   {174};
-constexpr int N_6NT   {175};
-constexpr int N_2N3HE {176};
-constexpr int N_3N3HE {177};
-constexpr int N_4N3HE {178};
-constexpr int N_3N2P  {179};
-constexpr int N_3N3A  {180};
-constexpr int N_3NPA  {181};
-constexpr int N_DT    {182};
-constexpr int N_NPD   {183};
-constexpr int N_NPT   {184};
-constexpr int N_NDT   {185};
-constexpr int N_NP3HE {186};
-constexpr int N_ND3HE {187};
-constexpr int N_NT3HE {188};
-constexpr int N_NTA   {189};
-constexpr int N_2N2P  {190};
-constexpr int N_P3HE  {191};
-constexpr int N_D3HE  {192};
-constexpr int N_3HEA  {193};
-constexpr int N_4N2P  {194};
-constexpr int N_4N2A  {195};
-constexpr int N_4NPA  {196};
-constexpr int N_3P    {197};
-constexpr int N_N3P   {198};
-constexpr int N_3N2PA {199};
-constexpr int N_5N2P  {200};
-constexpr int N_XP    {203};
-constexpr int N_XD    {204};
-constexpr int N_XT    {205};
-constexpr int N_X3HE  {206};
-constexpr int N_XA    {207};
-constexpr int HEATING {301};
-constexpr int DAMAGE_ENERGY {444};
-constexpr int COHERENT {502};
-constexpr int INCOHERENT {504};
-constexpr int PAIR_PROD_ELEC {515};
-constexpr int PAIR_PROD {516};
-constexpr int PAIR_PROD_NUC {517};
-constexpr int PHOTOELECTRIC {522};
-constexpr int N_P0    {600};
-constexpr int N_PC    {649};
-constexpr int N_D0    {650};
-constexpr int N_DC    {699};
-constexpr int N_T0    {700};
-constexpr int N_TC    {749};
-constexpr int N_3HE0  {750};
-constexpr int N_3HEC  {799};
-constexpr int N_A0    {800};
-constexpr int N_AC    {849};
-constexpr int N_2N0   {875};
-constexpr int N_2NC   {891};
-constexpr int HEATING_LOCAL {901};
+enum ReactionType {
+  REACTION_NONE = 0,
+  TOTAL_XS = 1,
+  ELASTIC = 2,
+  N_NONELASTIC = 3,
+  N_LEVEL = 4,
+  MISC    = 5,
+  N_2ND   = 11,
+  N_2N    = 16,
+  N_3N    = 17,
+  N_FISSION = 18,
+  N_F     = 19,
+  N_NF    = 20,
+  N_2NF   = 21,
+  N_NA    = 22,
+  N_N3A   = 23,
+  N_2NA   = 24,
+  N_3NA   = 25,
+  N_NP    = 28,
+  N_N2A   = 29,
+  N_2N2A  = 30,
+  N_ND    = 32,
+  N_NT    = 33,
+  N_N3HE  = 34,
+  N_ND2A  = 35,
+  N_NT2A  = 36,
+  N_4N    = 37,
+  N_3NF   = 38,
+  N_2NP   = 41,
+  N_3NP   = 42,
+  N_N2P   = 44,
+  N_NPA   = 45,
+  N_N1    = 51,
+  N_N40   = 90,
+  N_NC    = 91,
+  N_DISAPPEAR = 101,
+  N_GAMMA = 102,
+  N_P     = 103,
+  N_D     = 104,
+  N_T     = 105,
+  N_3HE   = 106,
+  N_A     = 107,
+  N_2A    = 108,
+  N_3A    = 109,
+  N_2P    = 111,
+  N_PA    = 112,
+  N_T2A   = 113,
+  N_D2A   = 114,
+  N_PD    = 115,
+  N_PT    = 116,
+  N_DA    = 117,
+  N_5N    = 152,
+  N_6N    = 153,
+  N_2NT   = 154,
+  N_TA    = 155,
+  N_4NP   = 156,
+  N_3ND   = 157,
+  N_NDA   = 158,
+  N_2NPA  = 159,
+  N_7N    = 160,
+  N_8N    = 161,
+  N_5NP   = 162,
+  N_6NP   = 163,
+  N_7NP   = 164,
+  N_4NA   = 165,
+  N_5NA   = 166,
+  N_6NA   = 167,
+  N_7NA   = 168,
+  N_4ND   = 169,
+  N_5ND   = 170,
+  N_6ND   = 171,
+  N_3NT   = 172,
+  N_4NT   = 173,
+  N_5NT   = 174,
+  N_6NT   = 175,
+  N_2N3HE = 176,
+  N_3N3HE = 177,
+  N_4N3HE = 178,
+  N_3N2P  = 179,
+  N_3N3A  = 180,
+  N_3NPA  = 181,
+  N_DT    = 182,
+  N_NPD   = 183,
+  N_NPT   = 184,
+  N_NDT   = 185,
+  N_NP3HE = 186,
+  N_ND3HE = 187,
+  N_NT3HE = 188,
+  N_NTA   = 189,
+  N_2N2P  = 190,
+  N_P3HE  = 191,
+  N_D3HE  = 192,
+  N_3HEA  = 193,
+  N_4N2P  = 194,
+  N_4N2A  = 195,
+  N_4NPA  = 196,
+  N_3P    = 197,
+  N_N3P   = 198,
+  N_3N2PA = 199,
+  N_5N2P  = 200,
+  N_XP    = 203,
+  N_XD    = 204,
+  N_XT    = 205,
+  N_X3HE  = 206,
+  N_XA    = 207,
+  HEATING = 301,
+  DAMAGE_ENERGY = 444,
+  COHERENT = 502,
+  INCOHERENT = 504,
+  PAIR_PROD_ELEC = 515,
+  PAIR_PROD = 516,
+  PAIR_PROD_NUC = 517,
+  PHOTOELECTRIC = 522,
+  N_P0    = 600,
+  N_PC    = 649,
+  N_D0    = 650,
+  N_DC    = 699,
+  N_T0    = 700,
+  N_TC    = 749,
+  N_3HE0  = 750,
+  N_3HEC  = 799,
+  N_A0    = 800,
+  N_AC    = 849,
+  N_2N0   = 875,
+  N_2NC   = 891,
+  HEATING_LOCAL = 901
+};
 
 constexpr std::array<int, 6> DEPLETION_RX {N_GAMMA, N_P, N_A, N_2N, N_3N, N_4N};
 
-// Fission neutron emission (nu) type
-constexpr int NU_NONE       {0}; // No nu values (non-fissionable)
-constexpr int NU_POLYNOMIAL {1}; // Nu values given by polynomial
-constexpr int NU_TABULAR    {2}; // Nu values given by tabular distribution
-
-// Library types
-constexpr int LIBRARY_NEUTRON {1};
-constexpr int LIBRARY_THERMAL {2};
-constexpr int LIBRARY_PHOTON {3};
-constexpr int LIBRARY_MULTIGROUP {4};
-
-// Probability table parameters
-constexpr int URR_CUM_PROB {0};
-constexpr int URR_TOTAL    {1};
-constexpr int URR_ELASTIC  {2};
-constexpr int URR_FISSION  {3};
-constexpr int URR_N_GAMMA  {4};
-constexpr int URR_HEATING  {5};
+enum class UnresolvProbTableParam {
+  CUM_PROB,
+  TOTAL,
+  ELASTIC,
+  FISSION,
+  N_GAMMA,
+  HEATING 
+};
 
 // Maximum number of partial fission reactions
 constexpr int PARTIAL_FISSION_MAX {4};
 
 // Resonance elastic scattering methods
-// TODO: Convert to enum
 enum class ResScatMethod {
   rvs, // Relative velocity sampling
   dbrc, // Doppler broadening rejection correction
   cxs // Constant cross section
 };
 
-// Electron treatments
-// TODO: Convert to enum
-constexpr int ELECTRON_LED {1}; // Local Energy Deposition
-constexpr int ELECTRON_TTB {2}; // Thick Target Bremsstrahlung
+enum class ElectronTreatment {
+  LED, // Local Energy Deposition
+  TTB  // Thick Target Bremsstrahlung
+};
 
 // ============================================================================
 // MULTIGROUP RELATED
 
-// MGXS Table Types
-// TODO: Convert to enum
-constexpr int MGXS_ISOTROPIC {1}; // Isotroically weighted data
-constexpr int MGXS_ANGLE {2};     // Data by angular bins
-
 // Flag to denote this was a macroscopic data object
 constexpr double MACROSCOPIC_AWR {-2.};
 
@@ -306,103 +284,86 @@ constexpr double MACROSCOPIC_AWR {-2.};
 constexpr int DEFAULT_NMU {33};
 
 // Mgxs::get_xs enumerated types
-// TODO: Convert to enum
-constexpr int MG_GET_XS_TOTAL              {0};
-constexpr int MG_GET_XS_ABSORPTION         {1};
-constexpr int MG_GET_XS_INVERSE_VELOCITY   {2};
-constexpr int MG_GET_XS_DECAY_RATE         {3};
-constexpr int MG_GET_XS_SCATTER            {4};
-constexpr int MG_GET_XS_SCATTER_MULT       {5};
-constexpr int MG_GET_XS_SCATTER_FMU_MULT   {6};
-constexpr int MG_GET_XS_SCATTER_FMU        {7};
-constexpr int MG_GET_XS_FISSION            {8};
-constexpr int MG_GET_XS_KAPPA_FISSION      {9};
-constexpr int MG_GET_XS_PROMPT_NU_FISSION  {10};
-constexpr int MG_GET_XS_DELAYED_NU_FISSION {11};
-constexpr int MG_GET_XS_NU_FISSION         {12};
-constexpr int MG_GET_XS_CHI_PROMPT         {13};
-constexpr int MG_GET_XS_CHI_DELAYED        {14};
+enum class MgxsType {
+  TOTAL,
+  ABSORPTION,
+  INVERSE_VELOCITY,
+  DECAY_RATE,
+  SCATTER,
+  SCATTER_MULT,
+  SCATTER_FMU_MULT,
+  SCATTER_FMU,
+  FISSION,
+  KAPPA_FISSION,
+  PROMPT_NU_FISSION,
+  DELAYED_NU_FISSION,
+  NU_FISSION,
+  CHI_PROMPT,
+  CHI_DELAYED
+};
 
 // ============================================================================
 // TALLY-RELATED CONSTANTS
 
-// Tally result entries
-constexpr int RESULT_VALUE  {0};
-constexpr int RESULT_SUM    {1};
-constexpr int RESULT_SUM_SQ {2};
+enum class TallyResult {
+  VALUE,
+  SUM,
+  SUM_SQ
+};
 
-// Tally type
-// TODO: Convert to enum
-constexpr int TALLY_VOLUME          {1};
-constexpr int TALLY_MESH_SURFACE    {2};
-constexpr int TALLY_SURFACE         {3};
+enum class TallyType {
+  VOLUME,
+  MESH_SURFACE,
+  SURFACE
+};
 
-// Tally estimator types
-// TODO: Convert to enum
-constexpr int ESTIMATOR_ANALOG      {1};
-constexpr int ESTIMATOR_TRACKLENGTH {2};
-constexpr int ESTIMATOR_COLLISION   {3};
+enum class TallyEstimator {
+  ANALOG,
+  TRACKLENGTH,
+  COLLISION
+};
 
-// Event types for tallies
-// TODO: Convert to enum
-constexpr int EVENT_SURFACE {-2};
-constexpr int EVENT_LATTICE {-1};
-constexpr int EVENT_KILL    {0};
-constexpr int EVENT_SCATTER {1};
-constexpr int EVENT_ABSORB  {2};
+enum class TallyEvent {
+  SURFACE,
+  LATTICE,
+  KILL,
+  SCATTER,
+  ABSORB 
+};
 
 // Tally score type -- if you change these, make sure you also update the
 // _SCORES dictionary in openmc/capi/tally.py
-// TODO: Convert to enum
-constexpr int SCORE_FLUX               {-1}; // flux
-constexpr int SCORE_TOTAL              {-2}; // total reaction rate
-constexpr int SCORE_SCATTER            {-3}; // scattering rate
-constexpr int SCORE_NU_SCATTER         {-4}; // scattering production rate
-constexpr int SCORE_ABSORPTION         {-5}; // absorption rate
-constexpr int SCORE_FISSION            {-6}; // fission rate
-constexpr int SCORE_NU_FISSION         {-7}; // neutron production rate
-constexpr int SCORE_KAPPA_FISSION      {-8}; // fission energy production rate
-constexpr int SCORE_CURRENT            {-9}; // current
-constexpr int SCORE_EVENTS             {-10}; // number of events
-constexpr int SCORE_DELAYED_NU_FISSION {-11}; // delayed neutron production rate
-constexpr int SCORE_PROMPT_NU_FISSION  {-12}; // prompt neutron production rate
-constexpr int SCORE_INVERSE_VELOCITY   {-13}; // flux-weighted inverse velocity
-constexpr int SCORE_FISS_Q_PROMPT      {-14}; // prompt fission Q-value
-constexpr int SCORE_FISS_Q_RECOV       {-15}; // recoverable fission Q-value
-constexpr int SCORE_DECAY_RATE         {-16}; // delayed neutron precursor decay rate
-
-// Tally map bin finding
-constexpr int NO_BIN_FOUND {-1};
-
-// Tally filter and map types
-// TODO: Refactor to remove or convert to enum
-constexpr int FILTER_UNIVERSE       {1};
-constexpr int FILTER_MATERIAL       {2};
-constexpr int FILTER_CELL           {3};
-
-// Mesh types
-constexpr int MESH_REGULAR {1};
-
-// Tally surface current directions
-constexpr int OUT_LEFT   {1};  // x min
-constexpr int IN_LEFT    {2};  // x min
-constexpr int OUT_RIGHT  {3};  // x max
-constexpr int IN_RIGHT   {4};  // x max
-constexpr int OUT_BACK   {5};  // y min
-constexpr int IN_BACK    {6};  // y min
-constexpr int OUT_FRONT  {7};  // y max
-constexpr int IN_FRONT   {8};  // y max
-constexpr int OUT_BOTTOM {9};  // z min
-constexpr int IN_BOTTOM  {10}; // z min
-constexpr int OUT_TOP    {11}; // z max
-constexpr int IN_TOP     {12}; // z max
+//
+// These are kept as a normal enum and made negative, since variables which
+// store one of these enum values usually also may be responsible for storing
+// MT numbers from the long enum above.
+enum TallyScore {
+  SCORE_FLUX = -1, // flux
+  SCORE_TOTAL = -2, // total reaction rate
+  SCORE_SCATTER = -3, // scattering rate
+  SCORE_NU_SCATTER = -4, // scattering production rate
+  SCORE_ABSORPTION = -5, // absorption rate
+  SCORE_FISSION = -6, // fission rate
+  SCORE_NU_FISSION = -7, // neutron production rate
+  SCORE_KAPPA_FISSION = -8, // fission energy production rate
+  SCORE_CURRENT = -9, // current
+  SCORE_EVENTS = -10, // number of events
+  SCORE_DELAYED_NU_FISSION = -11, // delayed neutron production rate
+  SCORE_PROMPT_NU_FISSION = -12, // prompt neutron production rate
+  SCORE_INVERSE_VELOCITY = -13, // flux-weighted inverse velocity
+  SCORE_FISS_Q_PROMPT = -14, // prompt fission Q-value
+  SCORE_FISS_Q_RECOV = -15, // recoverable fission Q-value
+  SCORE_DECAY_RATE = -16 // delayed neutron precursor decay rate
+};
 
 // Global tally parameters
 constexpr int N_GLOBAL_TALLIES {4};
-constexpr int K_COLLISION   {0};
-constexpr int K_ABSORPTION  {1};
-constexpr int K_TRACKLENGTH {2};
-constexpr int LEAKAGE       {3};
+enum class GlobalTally {
+  K_COLLISION, 
+  K_ABSORPTION, 
+  K_TRACKLENGTH, 
+  LEAKAGE       
+};
 
 // Miscellaneous
 constexpr int C_NONE {-1};
@@ -413,12 +374,14 @@ enum class Interpolation {
   histogram = 1, lin_lin = 2, lin_log = 3, log_lin = 4, log_log = 5
 };
 
-// Run modes
-constexpr int RUN_MODE_FIXEDSOURCE {1};
-constexpr int RUN_MODE_EIGENVALUE  {2};
-constexpr int RUN_MODE_PLOTTING    {3};
-constexpr int RUN_MODE_PARTICLE    {4};
-constexpr int RUN_MODE_VOLUME      {5};
+enum class RunMode {
+  UNSET, // default value, OpenMC throws error if left to this
+  FIXEDSOURCE,
+  EIGENVALUE,
+  PLOTTING,
+  PARTICLE,
+  VOLUME
+};
 
 // ============================================================================
 // CMFD CONSTANTS
diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h
index b7dc374fd..8d8ca4d8e 100644
--- a/include/openmc/mgxs.h
+++ b/include/openmc/mgxs.h
@@ -39,7 +39,7 @@ class Mgxs {
   private:
 
     xt::xtensor<double, 1> kTs;   // temperature in eV (k * T)
-    int scatter_format; // flag for if this is legendre, histogram, or tabular
+    AngleDistributionType scatter_format; // flag for if this is legendre, histogram, or tabular
     int num_delayed_groups; // number of delayed neutron groups
     int num_groups;     // number of energy groups
     std::vector<XsData> xs; // Cross section data
@@ -64,7 +64,7 @@ class Mgxs {
     //! @param in_azimuthal Azimuthal angle grid.
     void
     init(const std::string& in_name, double in_awr, const std::vector<double>& in_kTs,
-         bool in_fissionable, int in_scatter_format, bool in_is_isotropic,
+         bool in_fissionable, AngleDistributionType in_scatter_format, bool in_is_isotropic,
          const std::vector<double>& in_polar, const std::vector<double>& in_azimuthal);
 
     //! \brief Initializes the Mgxs object metadata from the HDF5 file
@@ -141,11 +141,11 @@ class Mgxs {
     //! @param dg delayed group index; use nullptr if irrelevant.
     //! @return Requested cross section value.
     double
-    get_xs(int xstype, int gin, const int* gout, const double* mu,
+    get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
       const int* dg);
 
     inline double
-    get_xs(int xstype, int gin)
+    get_xs(MgxsType xstype, int gin)
     {return get_xs(xstype, gin, nullptr, nullptr, nullptr);}
 
 
diff --git a/include/openmc/particle.h b/include/openmc/particle.h
index 568c925e9..1b049c2e2 100644
--- a/include/openmc/particle.h
+++ b/include/openmc/particle.h
@@ -263,7 +263,7 @@ public:
 
   // What event took place
   bool fission_ {false}; //!< did particle cause implicit fission
-  int event_;          //!< scatter, absorption
+  TallyEvent event_;          //!< scatter, absorption
   int event_nuclide_;  //!< index in nuclides array
   int event_mt_;       //!< reaction MT
   int delayed_group_ {0};  //!< delayed group
diff --git a/include/openmc/physics.h b/include/openmc/physics.h
index fb43775a6..048ba7998 100644
--- a/include/openmc/physics.h
+++ b/include/openmc/physics.h
@@ -27,15 +27,15 @@ void sample_neutron_reaction(Particle* p);
 void sample_photon_reaction(Particle* p);
 
 //! Terminates the particle and either deposits all energy locally
-//! (electron_treatment = ELECTRON_LED) or creates secondary bremsstrahlung
+//! (electron_treatment = ElectronTreatment::LED) or creates secondary bremsstrahlung
 //! photons from electron deflections with charged particles (electron_treatment
-//! = ELECTRON_TTB).
+//! = ElectronTreatment::TTB).
 void sample_electron_reaction(Particle* p);
 
 //! Terminates the particle and either deposits all energy locally
-//! (electron_treatment = ELECTRON_LED) or creates secondary bremsstrahlung
+//! (electron_treatment = ElectronTreatment::LED) or creates secondary bremsstrahlung
 //! photons from electron deflections with charged particles (electron_treatment
-//! = ELECTRON_TTB). Two annihilation photons of energy MASS_ELECTRON_EV (0.511
+//! = ElectronTreatment::TTB). Two annihilation photons of energy MASS_ELECTRON_EV (0.511
 //! MeV) are created and travel in opposite directions.
 void sample_positron_reaction(Particle* p);
 
diff --git a/include/openmc/scattdata.h b/include/openmc/scattdata.h
index f17e689bf..a7ab243fa 100644
--- a/include/openmc/scattdata.h
+++ b/include/openmc/scattdata.h
@@ -125,7 +125,7 @@ class ScattData {
     //!   use nullptr if irrelevant.
     //! @return Requested cross section value.
     double
-    get_xs(int xstype, int gin, const int* gout, const double* mu);
+    get_xs(MgxsType xstype, int gin, const int* gout, const double* mu);
 };
 
 //==============================================================================
diff --git a/include/openmc/settings.h b/include/openmc/settings.h
index 6570a25f7..c603f4cd0 100644
--- a/include/openmc/settings.h
+++ b/include/openmc/settings.h
@@ -66,7 +66,7 @@ extern "C" int32_t n_inactive;     //!< number of inactive batches
 extern "C" int32_t gen_per_batch;  //!< number of generations per batch
 extern "C" int64_t n_particles;    //!< number of particles per generation
 
-extern int electron_treatment;       //!< how to treat secondary electrons
+extern ElectronTreatment electron_treatment;       //!< how to treat secondary electrons
 extern std::array<double, 4> energy_cutoff;  //!< Energy cutoff in [eV] for each particle type
 extern int legendre_to_tabular_points; //!< number of points to convert Legendres
 extern int max_order;                //!< Maximum Legendre order for multigroup data
@@ -76,10 +76,10 @@ extern ResScatMethod res_scat_method; //!< resonance upscattering method
 extern double res_scat_energy_min;   //!< Min energy in [eV] for res. upscattering
 extern double res_scat_energy_max;   //!< Max energy in [eV] for res. upscattering
 extern std::vector<std::string> res_scat_nuclides;  //!< Nuclides using res. upscattering treatment
-extern "C" int run_mode;                 //!< Run mode (eigenvalue, fixed src, etc.)
+extern RunMode run_mode;                 //!< Run mode (eigenvalue, fixed src, etc.)
 extern std::unordered_set<int> sourcepoint_batch; //!< Batches when source should be written
 extern std::unordered_set<int> statepoint_batch; //!< Batches when state should be written
-extern int temperature_method;           //!< method for choosing temperatures
+extern TemperatureInterpolationType temperature_method;           //!< method for choosing temperatures
 extern double temperature_tolerance;     //!< Tolerance in [K] on choosing temperatures
 extern double temperature_default;       //!< Default T in [K]
 extern std::array<double, 2> temperature_range;  //!< Min/max T in [K] over which to load xs
diff --git a/include/openmc/surface.h b/include/openmc/surface.h
index 2d5000d28..f0c76e873 100644
--- a/include/openmc/surface.h
+++ b/include/openmc/surface.h
@@ -16,17 +16,6 @@
 
 namespace openmc {
 
-//==============================================================================
-// Module constant declarations (defined in .cpp)
-//==============================================================================
-
-// TODO: Convert to enum
-extern "C" const int BC_TRANSMIT;
-extern "C" const int BC_VACUUM;
-extern "C" const int BC_REFLECT;
-extern "C" const int BC_PERIODIC;
-extern "C" const int BC_WHITE;
-
 //==============================================================================
 // Global variables
 //==============================================================================
@@ -93,8 +82,18 @@ struct BoundingBox
 class Surface
 {
 public:
+
+  // Types of available boundary conditions on a surface
+  enum class Bc {
+    TRANSMIT,
+    VACUUM,
+    REFLECT,
+    PERIODIC,
+    WHITE
+  };
+
   int id_;                    //!< Unique ID
-  int bc_;                    //!< Boundary condition
+  Bc bc_;                    //!< Boundary condition
   std::string name_;          //!< User-defined name
 
   explicit Surface(pugi::xml_node surf_node);
diff --git a/include/openmc/tallies/derivative.h b/include/openmc/tallies/derivative.h
index 673afb0c9..e208fc432 100644
--- a/include/openmc/tallies/derivative.h
+++ b/include/openmc/tallies/derivative.h
@@ -14,9 +14,17 @@
 
 namespace openmc {
 
+// Different independent variables
+enum class WithRespectTo {
+  DENSITY,
+  NUCLIDE_DENSITY,
+  TEMPERATURE
+};
+
 struct TallyDerivative {
+
+  WithRespectTo variable;  //!< Independent variable (like temperature)
   int id;  //!< User-defined identifier
-  int variable;  //!< Independent variable (like temperature)
   int diff_material;  //!< Material this derivative is applied to
   int diff_nuclide;  //!< Nuclide this material is applied to
   double flux_deriv;  //!< Derivative of the current particle's weight
@@ -75,12 +83,6 @@ extern std::vector<TallyDerivative> 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};
-
 } // namespace openmc
 
 #endif // OPENMC_TALLIES_DERIVATIVE_H
diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h
index 7c46a24c2..b95d95bbb 100644
--- a/include/openmc/tallies/filter.h
+++ b/include/openmc/tallies/filter.h
@@ -79,7 +79,7 @@ public:
   //! \param[out] match will contain the matching bins and corresponding
   //!   weights; note that there may be zero matching bins
   virtual void
-  get_all_bins(const Particle* p, int estimator, FilterMatch& match) const = 0;
+  get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match) const = 0;
 
   //! Writes data describing this filter to an HDF5 statepoint group.
   virtual void
diff --git a/include/openmc/tallies/filter_azimuthal.h b/include/openmc/tallies/filter_azimuthal.h
index 84c2a2d8c..71fa5baa6 100644
--- a/include/openmc/tallies/filter_azimuthal.h
+++ b/include/openmc/tallies/filter_azimuthal.h
@@ -29,7 +29,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_cell.h b/include/openmc/tallies/filter_cell.h
index b570ee027..a5cc71685 100644
--- a/include/openmc/tallies/filter_cell.h
+++ b/include/openmc/tallies/filter_cell.h
@@ -30,7 +30,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_cell_instance.h b/include/openmc/tallies/filter_cell_instance.h
index 1b2cdd425..889881557 100644
--- a/include/openmc/tallies/filter_cell_instance.h
+++ b/include/openmc/tallies/filter_cell_instance.h
@@ -33,7 +33,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_cellborn.h b/include/openmc/tallies/filter_cellborn.h
index 706f6c47a..6300b9e1a 100644
--- a/include/openmc/tallies/filter_cellborn.h
+++ b/include/openmc/tallies/filter_cellborn.h
@@ -19,7 +19,7 @@ public:
 
   std::string type() const override {return "cellborn";}
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   std::string text_label(int bin) const override;
diff --git a/include/openmc/tallies/filter_cellfrom.h b/include/openmc/tallies/filter_cellfrom.h
index e86e34854..d64b80466 100644
--- a/include/openmc/tallies/filter_cellfrom.h
+++ b/include/openmc/tallies/filter_cellfrom.h
@@ -19,7 +19,7 @@ public:
 
   std::string type() const override {return "cellfrom";}
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   std::string text_label(int bin) const override;
diff --git a/include/openmc/tallies/filter_delayedgroup.h b/include/openmc/tallies/filter_delayedgroup.h
index 8a2bbeeaa..9f82cfbed 100644
--- a/include/openmc/tallies/filter_delayedgroup.h
+++ b/include/openmc/tallies/filter_delayedgroup.h
@@ -31,7 +31,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_distribcell.h b/include/openmc/tallies/filter_distribcell.h
index 9430a0906..145292eb2 100644
--- a/include/openmc/tallies/filter_distribcell.h
+++ b/include/openmc/tallies/filter_distribcell.h
@@ -26,7 +26,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_energy.h b/include/openmc/tallies/filter_energy.h
index af7601721..7b2c6ec75 100644
--- a/include/openmc/tallies/filter_energy.h
+++ b/include/openmc/tallies/filter_energy.h
@@ -28,7 +28,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
@@ -68,7 +68,7 @@ public:
 
   std::string type() const override {return "energyout";}
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   std::string text_label(int bin) const override;
diff --git a/include/openmc/tallies/filter_energyfunc.h b/include/openmc/tallies/filter_energyfunc.h
index df82e659e..7810a651f 100644
--- a/include/openmc/tallies/filter_energyfunc.h
+++ b/include/openmc/tallies/filter_energyfunc.h
@@ -33,7 +33,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_legendre.h b/include/openmc/tallies/filter_legendre.h
index 3a14ec3cf..24952e60b 100644
--- a/include/openmc/tallies/filter_legendre.h
+++ b/include/openmc/tallies/filter_legendre.h
@@ -26,7 +26,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_material.h b/include/openmc/tallies/filter_material.h
index 65cf832a3..064da2293 100644
--- a/include/openmc/tallies/filter_material.h
+++ b/include/openmc/tallies/filter_material.h
@@ -30,7 +30,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_mesh.h b/include/openmc/tallies/filter_mesh.h
index 98dec5d50..3d4ff2fa2 100644
--- a/include/openmc/tallies/filter_mesh.h
+++ b/include/openmc/tallies/filter_mesh.h
@@ -28,7 +28,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_meshsurface.h b/include/openmc/tallies/filter_meshsurface.h
index 19d178c7f..699405225 100644
--- a/include/openmc/tallies/filter_meshsurface.h
+++ b/include/openmc/tallies/filter_meshsurface.h
@@ -13,7 +13,7 @@ public:
 
   std::string type() const override {return "meshsurface";}
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   std::string text_label(int bin) const override;
@@ -22,6 +22,21 @@ public:
   // Accessors
 
   void set_mesh(int32_t mesh) override;
+
+  enum class MeshDir {
+    OUT_LEFT,  // x min
+    IN_LEFT,  // x min
+    OUT_RIGHT,  // x max
+    IN_RIGHT,  // x max
+    OUT_BACK,  // y min
+    IN_BACK,  // y min
+    OUT_FRONT,  // y max
+    IN_FRONT,  // y max
+    OUT_BOTTOM,  // z min
+    IN_BOTTOM, // z min
+    OUT_TOP, // z max
+    IN_TOP // z max
+  };
 };
 
 } // namespace openmc
diff --git a/include/openmc/tallies/filter_mu.h b/include/openmc/tallies/filter_mu.h
index ae9c3e06a..810f33368 100644
--- a/include/openmc/tallies/filter_mu.h
+++ b/include/openmc/tallies/filter_mu.h
@@ -29,7 +29,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_particle.h b/include/openmc/tallies/filter_particle.h
index 61aa3106b..c2a27aaaa 100644
--- a/include/openmc/tallies/filter_particle.h
+++ b/include/openmc/tallies/filter_particle.h
@@ -27,7 +27,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_polar.h b/include/openmc/tallies/filter_polar.h
index 965950456..ef6ceaec1 100644
--- a/include/openmc/tallies/filter_polar.h
+++ b/include/openmc/tallies/filter_polar.h
@@ -29,7 +29,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_sph_harm.h b/include/openmc/tallies/filter_sph_harm.h
index 80f998b2f..aec21e117 100644
--- a/include/openmc/tallies/filter_sph_harm.h
+++ b/include/openmc/tallies/filter_sph_harm.h
@@ -32,7 +32,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_sptl_legendre.h b/include/openmc/tallies/filter_sptl_legendre.h
index fea727515..c39bd0100 100644
--- a/include/openmc/tallies/filter_sptl_legendre.h
+++ b/include/openmc/tallies/filter_sptl_legendre.h
@@ -30,7 +30,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_surface.h b/include/openmc/tallies/filter_surface.h
index 23b9be33d..417e770d0 100644
--- a/include/openmc/tallies/filter_surface.h
+++ b/include/openmc/tallies/filter_surface.h
@@ -30,7 +30,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_universe.h b/include/openmc/tallies/filter_universe.h
index fc2ace18f..ce93b79e5 100644
--- a/include/openmc/tallies/filter_universe.h
+++ b/include/openmc/tallies/filter_universe.h
@@ -30,7 +30,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
diff --git a/include/openmc/tallies/filter_zernike.h b/include/openmc/tallies/filter_zernike.h
index e3ae89dec..3108bd487 100644
--- a/include/openmc/tallies/filter_zernike.h
+++ b/include/openmc/tallies/filter_zernike.h
@@ -26,7 +26,7 @@ public:
 
   void from_xml(pugi::xml_node node) override;
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   void to_statepoint(hid_t filter_group) const override;
@@ -76,7 +76,7 @@ public:
 
   std::string type() const override {return "zernikeradial";}
 
-  void get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+  void get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
   const override;
 
   std::string text_label(int bin) const override;
diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h
index aebade144..e8f2ff5e5 100644
--- a/include/openmc/tallies/tally.h
+++ b/include/openmc/tallies/tally.h
@@ -79,10 +79,10 @@ public:
 
   std::string name_; //!< User-defined name
 
-  int type_ {TALLY_VOLUME}; //!< e.g. volume, surface current
+  TallyType type_ {TallyType::VOLUME}; //!< e.g. volume, surface current
 
   //! Event type that contributes to this tally
-  int estimator_ {ESTIMATOR_TRACKLENGTH};
+  TallyEstimator estimator_ {TallyEstimator::TRACKLENGTH};
 
   //! Whether this tally is currently being updated
   bool active_ {false};
diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h
index 16be335b6..ca2f78562 100644
--- a/include/openmc/thermal.h
+++ b/include/openmc/thermal.h
@@ -16,21 +16,6 @@
 
 namespace openmc {
 
-//==============================================================================
-// Constants
-//==============================================================================
-
-// Secondary energy mode for S(a,b) inelastic scattering
-// TODO: Convert to enum
-constexpr int SAB_SECONDARY_EQUAL  {0}; // Equally-likely outgoing energy bins
-constexpr int SAB_SECONDARY_SKEWED {1}; // Skewed outgoing energy bins
-constexpr int SAB_SECONDARY_CONT   {2}; // Continuous, linear-linear interpolation
-
-// Elastic mode for S(a,b) elastic scattering
-// TODO: Convert to enum
-constexpr int SAB_ELASTIC_INCOHERENT {3}; // Incoherent elastic scattering
-constexpr int SAB_ELASTIC_COHERENT   {4}; // Coherent elastic scattering (Bragg edges)
-
 //==============================================================================
 // Global variables
 //==============================================================================
diff --git a/include/openmc/volume_calc.h b/include/openmc/volume_calc.h
index 1b49da814..e025ab3b0 100644
--- a/include/openmc/volume_calc.h
+++ b/include/openmc/volume_calc.h
@@ -46,8 +46,15 @@ public:
   //! \param[in] results   Vector of results for each domain
   void to_hdf5(const std::string& filename, const std::vector<Result>& results) const;
 
+  // Tally filter and map types
+  enum class TallyDomain {
+    UNIVERSE,
+    MATERIAL,
+    CELL
+  };
+
   // Data members
-  int domain_type_; //!< Type of domain (cell, material, etc.)
+  TallyDomain domain_type_; //!< Type of domain (cell, material, etc.)
   size_t n_samples_; //!< Number of samples to use
   double threshold_ {-1.0}; //!< Error threshold for domain volumes
   TriggerMetric trigger_type_ {TriggerMetric::not_active}; //!< Trigger metric for the volume calculation
diff --git a/include/openmc/xsdata.h b/include/openmc/xsdata.h
index 3e9fc51cc..64d0c01cf 100644
--- a/include/openmc/xsdata.h
+++ b/include/openmc/xsdata.h
@@ -15,6 +15,15 @@
 
 namespace openmc {
 
+// Angular distribution type
+enum class AngleDistributionType {
+  ISOTROPIC,
+  EQUI_32,
+  TABULAR,
+  LEGENDRE,
+  HISTOGRAM
+};
+
 //==============================================================================
 // XSDATA contains the temperature-independent cross section data for an MGXS
 //==============================================================================
@@ -25,8 +34,8 @@ class XsData {
 
     //! \brief Reads scattering data from the HDF5 file
     void
-    scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
-     int scatter_format, int final_scatter_format, int order_data);
+    scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, AngleDistributionType scatter_format,
+        AngleDistributionType final_scatter_format, int order_data);
 
     //! \brief Reads fission data from the HDF5 file
     void
@@ -104,8 +113,8 @@ class XsData {
     //! @param n_azi Number of azimuthal angles.
     //! @param n_groups Number of energy groups.
     //! @param n_d_groups Number of delayed neutron groups.
-    XsData(bool fissionable, int scatter_format, int n_pol, int n_azi,
-           size_t n_groups, size_t n_d_groups);
+    XsData(bool fissionable, AngleDistributionType scatter_format, int n_pol,
+        int n_azi, size_t n_groups, size_t n_d_groups);
 
     //! \brief Loads the XsData object from the HDF5 file
     //!
@@ -121,8 +130,8 @@ class XsData {
     //! @param n_pol Number of polar angles.
     //! @param n_azi Number of azimuthal angles.
     void
-    from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
-         int final_scatter_format, int order_data, bool is_isotropic, int n_pol,
+    from_hdf5(hid_t xsdata_grp, bool fissionable, AngleDistributionType scatter_format,
+         AngleDistributionType final_scatter_format, int order_data, bool is_isotropic, int n_pol,
          int n_azi);
 
     //! \brief Combines the microscopic data to a macroscopic object.
diff --git a/openmc/lib/tally.py b/openmc/lib/tally.py
index 2e8d43ae8..df28cb5dc 100644
--- a/openmc/lib/tally.py
+++ b/openmc/lib/tally.py
@@ -98,10 +98,10 @@ _SCORES = {
     -15: 'fission-q-recoverable', -16: 'decay-rate'
 }
 _ESTIMATORS = {
-    1: 'analog', 2: 'tracklength', 3: 'collision'
+    0: 'analog', 1: 'tracklength', 2: 'collision'
 }
 _TALLY_TYPES = {
-    1: 'volume', 2: 'mesh-surface', 3: 'surface'
+    0: 'volume', 1: 'mesh-surface', 2: 'surface'
 }
 
 
diff --git a/src/cell.cpp b/src/cell.cpp
index 9496dbba2..fbcdab42c 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -247,7 +247,7 @@ Cell::temperature(int32_t instance) const
 void
 Cell::set_temperature(double T, int32_t instance)
 {
-  if (settings::temperature_method == TEMPERATURE_INTERPOLATION) {
+  if (settings::temperature_method == TemperatureInterpolationType::INTERPOLATION) {
     if (T < data::temperature_min) {
       throw std::runtime_error{"Temperature is below minimum temperature at "
         "which data is available."};
@@ -565,7 +565,7 @@ CSGCell::to_hdf5(hid_t cell_group) const
   }
 
   // Write fill information.
-  if (type_ == FILL_MATERIAL) {
+  if (type_ == Fill::MATERIAL) {
     write_dataset(group, "fill_type", "material");
     std::vector<int32_t> mat_ids;
     for (auto i_mat : material_) {
@@ -586,7 +586,7 @@ CSGCell::to_hdf5(hid_t cell_group) const
       temps.push_back(sqrtkT_val * sqrtkT_val / K_BOLTZMANN);
     write_dataset(group, "temperature", temps);
 
-  } else if (type_ == FILL_UNIVERSE) {
+  } else if (type_ == Fill::UNIVERSE) {
     write_dataset(group, "fill_type", "universe");
     write_dataset(group, "fill", model::universes[fill_]->id_);
     if (translation_ != Position(0, 0, 0)) {
@@ -601,7 +601,7 @@ CSGCell::to_hdf5(hid_t cell_group) const
       }
     }
 
-  } else if (type_ == FILL_LATTICE) {
+  } else if (type_ == Fill::LATTICE) {
     write_dataset(group, "fill_type", "lattice");
     write_dataset(group, "lattice", model::lattices[fill_]->id_);
   }
@@ -1047,8 +1047,8 @@ openmc_cell_get_fill(int32_t index, int* type, int32_t** indices, int32_t* n)
 {
   if (index >= 0 && index < model::cells.size()) {
     Cell& c {*model::cells[index]};
-    *type = c.type_;
-    if (c.type_ == FILL_MATERIAL) {
+    *type = static_cast<int>(c.type_);
+    if (c.type_ == Fill::MATERIAL) {
       *indices = c.material_.data();
       *n = c.material_.size();
     } else {
@@ -1066,10 +1066,11 @@ extern "C" int
 openmc_cell_set_fill(int32_t index, int type, int32_t n,
                      const int32_t* indices)
 {
+  Fill filltype = static_cast<Fill>(type);
   if (index >= 0 && index < model::cells.size()) {
     Cell& c {*model::cells[index]};
-    if (type == FILL_MATERIAL) {
-      c.type_ = FILL_MATERIAL;
+    if (filltype == Fill::MATERIAL) {
+      c.type_ = Fill::MATERIAL;
       c.material_.clear();
       for (int i = 0; i < n; i++) {
         int i_mat = indices[i];
@@ -1083,10 +1084,10 @@ openmc_cell_set_fill(int32_t index, int type, int32_t n,
         }
       }
       c.material_.shrink_to_fit();
-    } else if (type == FILL_UNIVERSE) {
-      c.type_ = FILL_UNIVERSE;
+    } else if (filltype == Fill::UNIVERSE) {
+      c.type_ = Fill::UNIVERSE;
     } else {
-      c.type_ = FILL_LATTICE;
+      c.type_ = Fill::LATTICE;
     }
   } else {
     set_errmsg("Index in cells array is out of bounds.");
diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp
index dca57b451..22bf9c1ed 100644
--- a/src/cross_sections.cpp
+++ b/src/cross_sections.cpp
@@ -324,7 +324,7 @@ read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
   simulation::log_spacing = std::log(data::energy_max[neutron] /
     data::energy_min[neutron]) / settings::n_log_bins;
 
-  if (settings::photon_transport && settings::electron_treatment == ELECTRON_TTB) {
+  if (settings::photon_transport && settings::electron_treatment == ElectronTreatment::TTB) {
     // Determine if minimum/maximum energy for bremsstrahlung is greater/less
     // than the current minimum/maximum
     if (data::ttb_e_grid.size() >= 1) {
diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp
index 110f0100e..1fc1d4ca9 100644
--- a/src/eigenvalue.cpp
+++ b/src/eigenvalue.cpp
@@ -55,7 +55,7 @@ void calculate_generation_keff()
   const auto& gt = simulation::global_tallies;
 
   // Get keff for this generation by subtracting off the starting value
-  simulation::keff_generation = gt(K_TRACKLENGTH, RESULT_VALUE) - simulation::keff_generation;
+  simulation::keff_generation = gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) - simulation::keff_generation;
 
   double keff_reduced;
 #ifdef OPENMC_MPI
@@ -409,12 +409,12 @@ int openmc_get_keff(double* k_combined)
 
   std::array<double, 3> kv {};
   xt::xtensor<double, 2> cov = xt::zeros<double>({3, 3});
-  kv[0] = gt(K_COLLISION, RESULT_SUM) / n;
-  kv[1] = gt(K_ABSORPTION, RESULT_SUM) / n;
-  kv[2] = gt(K_TRACKLENGTH, RESULT_SUM) / n;
-  cov(0, 0) = (gt(K_COLLISION, RESULT_SUM_SQ) - n*kv[0]*kv[0]) / (n - 1);
-  cov(1, 1) = (gt(K_ABSORPTION, RESULT_SUM_SQ) - n*kv[1]*kv[1]) / (n - 1);
-  cov(2, 2) = (gt(K_TRACKLENGTH, RESULT_SUM_SQ) - n*kv[2]*kv[2]) / (n - 1);
+  kv[0] = gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n;
+  kv[1] = gt(GlobalTally::K_ABSORPTION, TallyResult::SUM) / n;
+  kv[2] = gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM) / n;
+  cov(0, 0) = (gt(GlobalTally::K_COLLISION, TallyResult::SUM_SQ) - n*kv[0]*kv[0]) / (n - 1);
+  cov(1, 1) = (gt(GlobalTally::K_ABSORPTION, TallyResult::SUM_SQ) - n*kv[1]*kv[1]) / (n - 1);
+  cov(2, 2) = (gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM_SQ) - n*kv[2]*kv[2]) / (n - 1);
 
   // Calculate covariances based on sums with Bessel's correction
   cov(0, 1) = (simulation::k_col_abs - n * kv[0] * kv[1]) / (n - 1);
diff --git a/src/finalize.cpp b/src/finalize.cpp
index d359f3885..d45cd1a7b 100644
--- a/src/finalize.cpp
+++ b/src/finalize.cpp
@@ -69,7 +69,7 @@ int openmc_finalize()
   settings::check_overlaps = false;
   settings::confidence_intervals = false;
   settings::create_fission_neutrons = true;
-  settings::electron_treatment = ELECTRON_LED;
+  settings::electron_treatment = ElectronTreatment::LED;
   settings::energy_cutoff = {0.0, 1000.0, 0.0, 0.0};
   settings::entropy_on = false;
   settings::gen_per_batch = 1;
@@ -88,14 +88,14 @@ int openmc_finalize()
   settings::res_scat_energy_max = 1000.0;
   settings::restart_run = false;
   settings::run_CE = true;
-  settings::run_mode = -1;
+  settings::run_mode = RunMode::UNSET;
   settings::dagmc = false;
   settings::source_latest = false;
   settings::source_separate = false;
   settings::source_write = true;
   settings::survival_biasing = false;
   settings::temperature_default = 293.6;
-  settings::temperature_method = TEMPERATURE_NEAREST;
+  settings::temperature_method = TemperatureInterpolationType::NEAREST;
   settings::temperature_multipole = false;
   settings::temperature_range = {0.0, 0.0};
   settings::temperature_tolerance = 10.0;
diff --git a/src/geometry.cpp b/src/geometry.cpp
index 6bcfc9c58..48a75af46 100644
--- a/src/geometry.cpp
+++ b/src/geometry.cpp
@@ -127,7 +127,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
 
   if (found) {
     Cell& c {*model::cells[i_cell]};
-    if (c.type_ == FILL_MATERIAL) {
+    if (c.type_ == Fill::MATERIAL) {
       //=======================================================================
       //! Found a material cell which means this is the lowest coord level.
 
@@ -136,9 +136,9 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
       if (c.distribcell_index_ >= 0) {
         for (int i = 0; i < p->n_coord_; i++) {
           const auto& c_i {*model::cells[p->coord_[i].cell]};
-          if (c_i.type_ == FILL_UNIVERSE) {
+          if (c_i.type_ == Fill::UNIVERSE) {
             offset += c_i.offset_[c.distribcell_index_];
-          } else if (c_i.type_ == FILL_LATTICE) {
+          } else if (c_i.type_ == Fill::LATTICE) {
             auto& lat {*model::lattices[p->coord_[i+1].lattice]};
             int i_xyz[3] {p->coord_[i+1].lattice_x,
                           p->coord_[i+1].lattice_y,
@@ -167,7 +167,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
 
       return true;
 
-    } else if (c.type_ == FILL_UNIVERSE) {
+    } else if (c.type_ == Fill::UNIVERSE) {
       //========================================================================
       //! Found a lower universe, update this coord level then search the next.
 
@@ -191,7 +191,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list)
       ++p->n_coord_;
       return find_cell_inner(p, nullptr);
 
-    } else if (c.type_ == FILL_LATTICE) {
+    } else if (c.type_ == Fill::LATTICE) {
       //========================================================================
       //! Found a lower lattice, update this coord level then search the next.
 
diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp
index 011dc2075..6ee289949 100644
--- a/src/geometry_aux.cpp
+++ b/src/geometry_aux.cpp
@@ -88,10 +88,10 @@ adjust_indices()
       auto search_univ = model::universe_map.find(id);
       auto search_lat = model::lattice_map.find(id);
       if (search_univ != model::universe_map.end()) {
-        c->type_ = FILL_UNIVERSE;
+        c->type_ = Fill::UNIVERSE;
         c->fill_ = search_univ->second;
       } else if (search_lat != model::lattice_map.end()) {
-        c->type_ = FILL_LATTICE;
+        c->type_ = Fill::LATTICE;
         c->fill_ = search_lat->second;
       } else {
         std::stringstream err_msg;
@@ -100,7 +100,7 @@ adjust_indices()
         fatal_error(err_msg);
       }
     } else {
-      c->type_ = FILL_MATERIAL;
+      c->type_ = Fill::MATERIAL;
       for (auto& mat_id : c->material_) {
         if (mat_id != MATERIAL_VOID) {
           auto search = model::material_map.find(mat_id);
@@ -334,7 +334,7 @@ prepare_distribcell()
   // By default, add material cells to the list of distributed cells
   if (settings::material_cell_offsets) {
     for (gsl::index i = 0; i < model::cells.size(); ++i) {
-      if (model::cells[i]->type_ == FILL_MATERIAL) distribcells.insert(i);
+      if (model::cells[i]->type_ == Fill::MATERIAL) distribcells.insert(i);
     }
   }
 
@@ -382,7 +382,7 @@ prepare_distribcell()
   // Allocate the cell and lattice offset tables.
   int n_maps = target_univ_ids.size();
   for (auto& c : model::cells) {
-    if (c->type_ != FILL_MATERIAL) {
+    if (c->type_ != Fill::MATERIAL) {
       c->offset_.resize(n_maps, C_NONE);
     }
   }
@@ -398,12 +398,12 @@ prepare_distribcell()
       for (int32_t cell_indx : univ->cells_) {
         Cell& c = *model::cells[cell_indx];
 
-        if (c.type_ == FILL_UNIVERSE) {
+        if (c.type_ == Fill::UNIVERSE) {
           c.offset_[map] = offset;
           int32_t search_univ = c.fill_;
           offset += count_universe_instances(search_univ, target_univ_id);
 
-        } else if (c.type_ == FILL_LATTICE) {
+        } else if (c.type_ == Fill::LATTICE) {
           Lattice& lat = *model::lattices[c.fill_];
           offset = lat.fill_offset_table(offset, target_univ_id, map);
         }
@@ -429,11 +429,11 @@ count_cell_instances(int32_t univ_indx)
       ++c.n_instances_;
       model::universe_cell_counts[univ_indx][cell_indx] += 1;
 
-      if (c.type_ == FILL_UNIVERSE) {
+      if (c.type_ == Fill::UNIVERSE) {
         // This cell contains another universe.  Recurse into that universe.
         count_cell_instances(c.fill_);
         update_universe_cell_count(univ_indx, c.fill_);
-      } else if (c.type_ == FILL_LATTICE) {
+      } else if (c.type_ == Fill::LATTICE) {
         // This cell contains a lattice.  Recurse into the lattice universes.
         Lattice& lat = *model::lattices[c.fill_];
         for (auto it = lat.begin(); it != lat.end(); ++it) {
@@ -459,11 +459,11 @@ count_universe_instances(int32_t search_univ, int32_t target_univ_id)
   for (int32_t cell_indx : model::universes[search_univ]->cells_) {
     Cell& c = *model::cells[cell_indx];
 
-    if (c.type_ == FILL_UNIVERSE) {
+    if (c.type_ == Fill::UNIVERSE) {
       int32_t next_univ = c.fill_;
       count += count_universe_instances(next_univ, target_univ_id);
 
-    } else if (c.type_ == FILL_LATTICE) {
+    } else if (c.type_ == Fill::LATTICE) {
       Lattice& lat = *model::lattices[c.fill_];
       for (auto it = lat.begin(); it != lat.end(); ++it) {
         int32_t next_univ = *it;
@@ -504,9 +504,9 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
     Cell& c = *model::cells[*cell_it];
 
     // Material cells don't contain other cells so ignore them.
-    if (c.type_ != FILL_MATERIAL) {
+    if (c.type_ != Fill::MATERIAL) {
       int32_t temp_offset;
-      if (c.type_ == FILL_UNIVERSE) {
+      if (c.type_ == Fill::UNIVERSE) {
         temp_offset = offset + c.offset_[map];
       } else {
         Lattice& lat = *model::lattices[c.fill_];
@@ -524,7 +524,7 @@ distribcell_path_inner(int32_t target_cell, int32_t map, int32_t target_offset,
   Cell& c = *model::cells[*cell_it];
   path << "c" << c.id_ << "->";
 
-  if (c.type_ == FILL_UNIVERSE) {
+  if (c.type_ == Fill::UNIVERSE) {
     // Recurse into the fill cell.
     offset += c.offset_[map];
     path << distribcell_path_inner(target_cell, map, target_offset,
@@ -571,10 +571,10 @@ maximum_levels(int32_t univ)
 
   for (int32_t cell_indx : model::universes[univ]->cells_) {
     Cell& c = *model::cells[cell_indx];
-    if (c.type_ == FILL_UNIVERSE) {
+    if (c.type_ == Fill::UNIVERSE) {
       int32_t next_univ = c.fill_;
       levels_below = std::max(levels_below, maximum_levels(next_univ));
-    } else if (c.type_ == FILL_LATTICE) {
+    } else if (c.type_ == Fill::LATTICE) {
       Lattice& lat = *model::lattices[c.fill_];
       for (auto it = lat.begin(); it != lat.end(); ++it) {
         int32_t next_univ = *it;
diff --git a/src/initialize.cpp b/src/initialize.cpp
index 460db2436..38b8e5c29 100644
--- a/src/initialize.cpp
+++ b/src/initialize.cpp
@@ -74,7 +74,7 @@ int openmc_init(int argc, char* argv[], const void* intracomm)
   read_input_xml();
 
   // Check for particle restart run
-  if (settings::particle_restart_run) settings::run_mode = RUN_MODE_PARTICLE;
+  if (settings::particle_restart_run) settings::run_mode = RunMode::PARTICLE;
 
   // Stop initialization timer
   simulation::time_initialize.stop();
@@ -126,7 +126,7 @@ parse_command_line(int argc, char* argv[])
     std::string arg {argv[i]};
     if (arg[0] == '-') {
       if (arg == "-p" || arg == "--plot") {
-        settings::run_mode = RUN_MODE_PLOTTING;
+        settings::run_mode = RunMode::PLOTTING;
         settings::check_overlaps = true;
 
       } else if (arg == "-n" || arg == "--particles") {
@@ -188,7 +188,7 @@ parse_command_line(int argc, char* argv[])
       } else if (arg == "-g" || arg == "--geometry-debug") {
       settings::check_overlaps = true;
       } else if (arg == "-c" || arg == "--volume") {
-        settings::run_mode = RUN_MODE_VOLUME;
+        settings::run_mode = RunMode::VOLUME;
       } else if (arg == "-s" || arg == "--threads") {
         // Read number of threads
         i += 1;
@@ -253,7 +253,7 @@ void read_input_xml()
   double_2dvec thermal_temps(data::thermal_scatt_map.size());
   finalize_geometry(nuc_temps, thermal_temps);
 
-  if (settings::run_mode != RUN_MODE_PLOTTING) {
+  if (settings::run_mode != RunMode::PLOTTING) {
     simulation::time_read_xs.start();
     if (settings::run_CE) {
       // Read continuous-energy cross sections
@@ -272,7 +272,7 @@ void read_input_xml()
   // Initialize distribcell_filters
   prepare_distribcell();
 
-  if (settings::run_mode == RUN_MODE_PLOTTING) {
+  if (settings::run_mode == RunMode::PLOTTING) {
     // Read plots.xml if it exists
     read_plots_xml();
     if (mpi::master && settings::verbosity >= 5) print_plot();
diff --git a/src/main.cpp b/src/main.cpp
index 936c89b78..36f7211a8 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -29,18 +29,18 @@ int main(int argc, char* argv[]) {
 
   // start problem based on mode
   switch (settings::run_mode) {
-    case RUN_MODE_FIXEDSOURCE:
-    case RUN_MODE_EIGENVALUE:
+    case RunMode::FIXEDSOURCE:
+    case RunMode::EIGENVALUE:
       err = openmc_run();
       break;
-    case RUN_MODE_PLOTTING:
+    case RunMode::PLOTTING:
       err = openmc_plot_geometry();
       break;
-    case RUN_MODE_PARTICLE:
+    case RunMode::PARTICLE:
       if (mpi::master) run_particle_restart();
       err = 0;
       break;
-    case RUN_MODE_VOLUME:
+    case RunMode::VOLUME:
       err = openmc_calculate_volumes();
       break;
   }
diff --git a/src/material.cpp b/src/material.cpp
index 38eef82a7..e7dd04c28 100644
--- a/src/material.cpp
+++ b/src/material.cpp
@@ -348,7 +348,7 @@ void Material::finalize()
   }
 
   // Generate material bremsstrahlung data for electrons and positrons
-  if (settings::photon_transport && settings::electron_treatment == ELECTRON_TTB) {
+  if (settings::photon_transport && settings::electron_treatment == ElectronTreatment::TTB) {
     this->init_bremsstrahlung();
   }
 
diff --git a/src/mgxs.cpp b/src/mgxs.cpp
index bba58f668..9dd3cb394 100644
--- a/src/mgxs.cpp
+++ b/src/mgxs.cpp
@@ -30,8 +30,8 @@ namespace openmc {
 
 void
 Mgxs::init(const std::string& in_name, double in_awr,
-     const std::vector<double>& in_kTs, bool in_fissionable, int in_scatter_format,
-     bool in_is_isotropic,
+     const std::vector<double>& in_kTs, bool in_fissionable,
+     AngleDistributionType in_scatter_format, bool in_is_isotropic,
      const std::vector<double>& in_polar, const std::vector<double>& in_azimuthal)
 {
   // Set the metadata
@@ -102,15 +102,15 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, const std::vector<double>& temperature,
 
   // If only one temperature is available, lets just use nearest temperature
   // interpolation
-  if ((num_temps == 1) && (settings::temperature_method == TEMPERATURE_INTERPOLATION)) {
+  if ((num_temps == 1) && (settings::temperature_method == TemperatureInterpolationType::INTERPOLATION)) {
     warning("Cross sections for " + strtrim(name) + " are only available " +
             "at one temperature.  Reverting to the nearest temperature " +
             "method.");
-    settings::temperature_method = TEMPERATURE_NEAREST;
+    settings::temperature_method = TemperatureInterpolationType::NEAREST;
   }
 
   switch(settings::temperature_method) {
-    case TEMPERATURE_NEAREST:
+    case TemperatureInterpolationType::NEAREST:
       // Determine actual temperatures to read
       for (const auto& T : temperature) {
         auto i_closest = xt::argmin(xt::abs(available_temps - T))[0];
@@ -129,7 +129,7 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, const std::vector<double>& temperature,
       }
       break;
 
-    case TEMPERATURE_INTERPOLATION:
+    case TemperatureInterpolationType::INTERPOLATION:
       for (int i = 0; i < temperature.size(); i++) {
         for (int j = 0; j < num_temps - 1; j++) {
           if ((available_temps[j] <= temperature[i]) &&
@@ -167,22 +167,22 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, const std::vector<double>& temperature,
   close_group(kT_group);
 
   // Load the remaining metadata
-  int in_scatter_format;
+  AngleDistributionType in_scatter_format;
   if (attribute_exists(xs_id, "scatter_format")) {
     std::string temp_str(MAX_WORD_LEN, ' ');
     read_attr_string(xs_id, "scatter_format", MAX_WORD_LEN, &temp_str[0]);
     to_lower(strtrim(temp_str));
     if (temp_str.compare(0, 8, "legendre") == 0) {
-      in_scatter_format = ANGLE_LEGENDRE;
+      in_scatter_format = AngleDistributionType::LEGENDRE;
     } else if (temp_str.compare(0, 9, "histogram") == 0) {
-      in_scatter_format = ANGLE_HISTOGRAM;
+      in_scatter_format = AngleDistributionType::HISTOGRAM;
     } else if (temp_str.compare(0, 7, "tabular") == 0) {
-      in_scatter_format = ANGLE_TABULAR;
+      in_scatter_format = AngleDistributionType::TABULAR;
     } else {
       fatal_error("Invalid scatter_format option!");
     }
   } else {
-    in_scatter_format = ANGLE_LEGENDRE;
+    in_scatter_format = AngleDistributionType::LEGENDRE;
   }
 
   if (attribute_exists(xs_id, "scatter_shape")) {
@@ -215,7 +215,7 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, const std::vector<double>& temperature,
   // However Pn has n+1 sets of points (since you need to count the P0
   // moment). Adjust for that. Histogram and Tabular formats dont need this
   // adjustment.
-  if (in_scatter_format == ANGLE_LEGENDRE) {
+  if (in_scatter_format == AngleDistributionType::LEGENDRE) {
     ++order_dim;
   }
 
@@ -281,9 +281,9 @@ Mgxs::Mgxs(hid_t xs_id, const std::vector<double>& temperature,
   metadata_from_hdf5(xs_id, temperature, temps_to_read, order_data);
 
   // Set number of energy and delayed groups
-  int final_scatter_format = scatter_format;
+  AngleDistributionType final_scatter_format = scatter_format;
   if (settings::legendre_to_tabular) {
-    if (scatter_format == ANGLE_LEGENDRE) final_scatter_format = ANGLE_TABULAR;
+    if (scatter_format == AngleDistributionType::LEGENDRE) final_scatter_format = AngleDistributionType::TABULAR;
   }
 
   // Load the more specific XsData information
@@ -322,7 +322,7 @@ Mgxs::Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
   }
   // Force all of the following data to be the same; these will be verified
   // to be true later
-  int in_scatter_format = micros[0]->scatter_format;
+  AngleDistributionType in_scatter_format = micros[0]->scatter_format;
   bool in_is_isotropic = micros[0]->is_isotropic;
   std::vector<double> in_polar = micros[0]->polar;
   std::vector<double> in_azimuthal = micros[0]->azimuthal;
@@ -344,7 +344,7 @@ Mgxs::Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
     std::vector<double> micro_t_interp(micros.size(), 0.);
     for (int m = 0; m < micros.size(); m++) {
       switch(settings::temperature_method) {
-      case TEMPERATURE_NEAREST:
+      case TemperatureInterpolationType::NEAREST:
         {
           micro_t[m] = xt::argmin(xt::abs(micros[m]->kTs - temp_desired))[0];
           auto temp_actual = micros[m]->kTs[micro_t[m]];
@@ -357,7 +357,7 @@ Mgxs::Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
           }
         }
         break;
-      case TEMPERATURE_INTERPOLATION:
+      case TemperatureInterpolationType::INTERPOLATION:
         // Get a list of bounding temperatures for each actual temperature
         // present in the model
         for (int k = 0; k < micros[m]->kTs.shape()[0] - 1; k++) {
@@ -381,7 +381,7 @@ Mgxs::Mgxs(const std::string& in_name, const std::vector<double>& mat_kTs,
     // If we are doing nearest temperature interpolation, then we don't need
     // to do the 2nd temperature
     int num_interp_points = 2;
-    if (settings::temperature_method == TEMPERATURE_NEAREST) num_interp_points = 1;
+    if (settings::temperature_method == TemperatureInterpolationType::NEAREST) num_interp_points = 1;
     for (int interp_point = 0; interp_point < num_interp_points; interp_point++) {
       std::vector<double> interp(micros.size());
       std::vector<double> temp_indices(micros.size());
@@ -416,7 +416,7 @@ Mgxs::combine(const std::vector<Mgxs*>& micros, const std::vector<double>& scala
 //==============================================================================
 
 double
-Mgxs::get_xs(int xstype, int gin, const int* gout, const double* mu,
+Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
   const int* dg)
 {
   // This method assumes that the temperature and angle indices are set
@@ -430,31 +430,31 @@ Mgxs::get_xs(int xstype, int gin, const int* gout, const double* mu,
 #endif
   double val;
   switch(xstype) {
-  case MG_GET_XS_TOTAL:
+  case MgxsType::TOTAL:
     val = xs_t->total(a, gin);
     break;
-  case MG_GET_XS_NU_FISSION:
+  case MgxsType::NU_FISSION:
     val = fissionable ? xs_t->nu_fission(a, gin) : 0.;
     break;
-  case MG_GET_XS_ABSORPTION:
+  case MgxsType::ABSORPTION:
     val = xs_t->absorption(a, gin);;
     break;
-  case MG_GET_XS_FISSION:
+  case MgxsType::FISSION:
     val = fissionable ? xs_t->fission(a, gin) : 0.;
     break;
-  case MG_GET_XS_KAPPA_FISSION:
+  case MgxsType::KAPPA_FISSION:
     val = fissionable ? xs_t->kappa_fission(a, gin) : 0.;
     break;
-  case MG_GET_XS_SCATTER:
-  case MG_GET_XS_SCATTER_MULT:
-  case MG_GET_XS_SCATTER_FMU_MULT:
-  case MG_GET_XS_SCATTER_FMU:
+  case MgxsType::SCATTER:
+  case MgxsType::SCATTER_MULT:
+  case MgxsType::SCATTER_FMU_MULT:
+  case MgxsType::SCATTER_FMU:
     val = xs_t->scatter[a]->get_xs(xstype, gin, gout, mu);
     break;
-  case MG_GET_XS_PROMPT_NU_FISSION:
+  case MgxsType::PROMPT_NU_FISSION:
     val = fissionable ? xs_t->prompt_nu_fission(a, gin) : 0.;
     break;
-  case MG_GET_XS_DELAYED_NU_FISSION:
+  case MgxsType::DELAYED_NU_FISSION:
     if (fissionable) {
       if (dg != nullptr) {
         val = xs_t->delayed_nu_fission(a, *dg, gin);
@@ -468,7 +468,7 @@ Mgxs::get_xs(int xstype, int gin, const int* gout, const double* mu,
       val = 0.;
     }
     break;
-  case MG_GET_XS_CHI_PROMPT:
+  case MgxsType::CHI_PROMPT:
     if (fissionable) {
       if (gout != nullptr) {
         val = xs_t->chi_prompt(a, gin, *gout);
@@ -483,7 +483,7 @@ Mgxs::get_xs(int xstype, int gin, const int* gout, const double* mu,
       val = 0.;
     }
     break;
-  case MG_GET_XS_CHI_DELAYED:
+  case MgxsType::CHI_DELAYED:
     if (fissionable) {
       if (gout != nullptr) {
         if (dg != nullptr) {
@@ -510,10 +510,10 @@ Mgxs::get_xs(int xstype, int gin, const int* gout, const double* mu,
       val = 0.;
     }
     break;
-  case MG_GET_XS_INVERSE_VELOCITY:
+  case MgxsType::INVERSE_VELOCITY:
     val = xs_t->inverse_velocity(a, gin);
     break;
-  case MG_GET_XS_DECAY_RATE:
+  case MgxsType::DECAY_RATE:
     if (dg != nullptr) {
       val = xs_t->decay_rate(a, *dg);
     } else {
diff --git a/src/nuclide.cpp b/src/nuclide.cpp
index 86f2ee140..62f1d1b05 100644
--- a/src/nuclide.cpp
+++ b/src/nuclide.cpp
@@ -69,12 +69,12 @@ Nuclide::Nuclide(hid_t group, const std::vector<double>& temperature, int i_nucl
   std::sort(temps_available.begin(), temps_available.end());
 
   // If only one temperature is available, revert to nearest temperature
-  if (temps_available.size() == 1 && settings::temperature_method == TEMPERATURE_INTERPOLATION) {
+  if (temps_available.size() == 1 && settings::temperature_method == TemperatureInterpolationType::INTERPOLATION) {
     if (mpi::master) {
       warning("Cross sections for " + name_ + " are only available at one "
         "temperature. Reverting to nearest temperature method.");
     }
-    settings::temperature_method = TEMPERATURE_NEAREST;
+    settings::temperature_method = TemperatureInterpolationType::NEAREST;
   }
 
   // Determine actual temperatures to read -- start by checking whether a
@@ -93,7 +93,7 @@ Nuclide::Nuclide(hid_t group, const std::vector<double>& temperature, int i_nucl
   }
 
   switch (settings::temperature_method) {
-  case TEMPERATURE_NEAREST:
+  case TemperatureInterpolationType::NEAREST:
     // Find nearest temperatures
     for (double T_desired : temperature) {
 
@@ -126,7 +126,7 @@ Nuclide::Nuclide(hid_t group, const std::vector<double>& temperature, int i_nucl
     }
     break;
 
-  case TEMPERATURE_INTERPOLATION:
+  case TemperatureInterpolationType::INTERPOLATION:
     // If temperature interpolation or multipole is selected, get a list of
     // bounding temperatures for each actual temperature present in the model
     for (double T_desired : temperature) {
@@ -560,7 +560,7 @@ void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle
     double f;
     int i_temp = -1;
     switch (settings::temperature_method) {
-    case TEMPERATURE_NEAREST:
+    case TemperatureInterpolationType::NEAREST:
       {
         double max_diff = INFTY;
         for (int t = 0; t < kTs_.size(); ++t) {
@@ -573,7 +573,7 @@ void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle
       }
       break;
 
-    case TEMPERATURE_INTERPOLATION:
+    case TemperatureInterpolationType::INTERPOLATION:
       // Find temperatures that bound the actual temperature
       for (i_temp = 0; i_temp < kTs_.size() - 1; ++i_temp) {
         if (kTs_[i_temp] <= kT && kT < kTs_[i_temp + 1]) break;
@@ -763,10 +763,10 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const
   p.stream_ = STREAM_TRACKING;
 
   int i_low = 0;
-  while (urr.prob_(i_energy, URR_CUM_PROB, i_low) <= r) {++i_low;};
+  while (urr.prob_(i_energy, UnresolvProbTableParam::CUM_PROB, i_low) <= r) {++i_low;};
 
   int i_up = 0;
-  while (urr.prob_(i_energy + 1, URR_CUM_PROB, i_up) <= r) {++i_up;};
+  while (urr.prob_(i_energy + 1, UnresolvProbTableParam::CUM_PROB, i_up) <= r) {++i_up;};
 
   // Determine elastic, fission, and capture cross sections from the
   // probability table
@@ -779,46 +779,46 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const
     f = (p.E_ - urr.energy_(i_energy)) /
          (urr.energy_(i_energy + 1) - urr.energy_(i_energy));
 
-    elastic = (1. - f) * urr.prob_(i_energy, URR_ELASTIC, i_low) +
-         f * urr.prob_(i_energy + 1, URR_ELASTIC, i_up);
-    fission = (1. - f) * urr.prob_(i_energy, URR_FISSION, i_low) +
-         f * urr.prob_(i_energy + 1, URR_FISSION, i_up);
-    capture = (1. - f) * urr.prob_(i_energy, URR_N_GAMMA, i_low) +
-         f * urr.prob_(i_energy + 1, URR_N_GAMMA, i_up);
+    elastic = (1. - f) * urr.prob_(i_energy, UnresolvProbTableParam::ELASTIC, i_low) +
+         f * urr.prob_(i_energy + 1, UnresolvProbTableParam::ELASTIC, i_up);
+    fission = (1. - f) * urr.prob_(i_energy, UnresolvProbTableParam::FISSION, i_low) +
+         f * urr.prob_(i_energy + 1, UnresolvProbTableParam::FISSION, i_up);
+    capture = (1. - f) * urr.prob_(i_energy, UnresolvProbTableParam::N_GAMMA, i_low) +
+         f * urr.prob_(i_energy + 1, UnresolvProbTableParam::N_GAMMA, i_up);
   } else if (urr.interp_ == Interpolation::log_log) {
     // Determine interpolation factor on the table
     f = std::log(p.E_ / urr.energy_(i_energy)) /
          std::log(urr.energy_(i_energy + 1) / urr.energy_(i_energy));
 
     // Calculate the elastic cross section/factor
-    if ((urr.prob_(i_energy, URR_ELASTIC, i_low) > 0.) &&
-        (urr.prob_(i_energy + 1, URR_ELASTIC, i_up) > 0.)) {
+    if ((urr.prob_(i_energy, UnresolvProbTableParam::ELASTIC, i_low) > 0.) &&
+        (urr.prob_(i_energy + 1, UnresolvProbTableParam::ELASTIC, i_up) > 0.)) {
       elastic =
            std::exp((1. - f) *
-                    std::log(urr.prob_(i_energy, URR_ELASTIC, i_low)) +
-                    f * std::log(urr.prob_(i_energy + 1, URR_ELASTIC, i_up)));
+                    std::log(urr.prob_(i_energy, UnresolvProbTableParam::ELASTIC, i_low)) +
+                    f * std::log(urr.prob_(i_energy + 1, UnresolvProbTableParam::ELASTIC, i_up)));
     } else {
       elastic = 0.;
     }
 
     // Calculate the fission cross section/factor
-    if ((urr.prob_(i_energy, URR_FISSION, i_low) > 0.) &&
-        (urr.prob_(i_energy + 1, URR_FISSION, i_up) > 0.)) {
+    if ((urr.prob_(i_energy, UnresolvProbTableParam::FISSION, i_low) > 0.) &&
+        (urr.prob_(i_energy + 1, UnresolvProbTableParam::FISSION, i_up) > 0.)) {
       fission =
            std::exp((1. - f) *
-                    std::log(urr.prob_(i_energy, URR_FISSION, i_low)) +
-                    f * std::log(urr.prob_(i_energy + 1, URR_FISSION, i_up)));
+                    std::log(urr.prob_(i_energy, UnresolvProbTableParam::FISSION, i_low)) +
+                    f * std::log(urr.prob_(i_energy + 1, UnresolvProbTableParam::FISSION, i_up)));
     } else {
       fission = 0.;
     }
 
     // Calculate the capture cross section/factor
-    if ((urr.prob_(i_energy, URR_N_GAMMA, i_low) > 0.) &&
-        (urr.prob_(i_energy + 1, URR_N_GAMMA, i_up) > 0.)) {
+    if ((urr.prob_(i_energy, UnresolvProbTableParam::N_GAMMA, i_low) > 0.) &&
+        (urr.prob_(i_energy + 1, UnresolvProbTableParam::N_GAMMA, i_up) > 0.)) {
       capture =
            std::exp((1. - f) *
-                    std::log(urr.prob_(i_energy, URR_N_GAMMA, i_low)) +
-                    f * std::log(urr.prob_(i_energy + 1, URR_N_GAMMA, i_up)));
+                    std::log(urr.prob_(i_energy, UnresolvProbTableParam::N_GAMMA, i_low)) +
+                    f * std::log(urr.prob_(i_energy + 1, UnresolvProbTableParam::N_GAMMA, i_up)));
     } else {
       capture = 0.;
     }
diff --git a/src/output.cpp b/src/output.cpp
index c8b91c4dc..76003ad9b 100644
--- a/src/output.cpp
+++ b/src/output.cpp
@@ -466,11 +466,11 @@ void print_runtime()
   show_time("Total time in simulation", time_inactive.elapsed() +
     time_active.elapsed());
   show_time("Time in transport only", time_transport.elapsed(), 1);
-  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+  if (settings::run_mode == RunMode::EIGENVALUE) {
     show_time("Time in inactive batches", time_inactive.elapsed(), 1);
   }
   show_time("Time in active batches", time_active.elapsed(), 1);
-  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+  if (settings::run_mode == RunMode::EIGENVALUE) {
     show_time("Time synchronizing fission bank", time_bank.elapsed(), 1);
     show_time("Sampling source sites", time_bank_sample.elapsed(), 2);
     show_time("SEND/RECV source sites", time_bank_sendrecv.elapsed(), 2);
@@ -524,8 +524,8 @@ void print_runtime()
 std::pair<double, double>
 mean_stdev(const double* x, int n)
 {
-  double mean = x[RESULT_SUM] / n;
-  double stdev = n > 1 ? std::sqrt((x[RESULT_SUM_SQ]/n
+  double mean = x[static_cast<int>(TallyResult::SUM)] / n;
+  double stdev = n > 1 ? std::sqrt((x[static_cast<int>(TallyResult::SUM_SQ)]/n
     - mean*mean)/(n - 1)) : 0.0;
   return {mean, stdev};
 }
@@ -560,14 +560,14 @@ void print_results()
   const auto& gt = simulation::global_tallies;
   double mean, stdev;
   if (n > 1) {
-    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
-      std::tie(mean, stdev) = mean_stdev(&gt(K_COLLISION, 0), n);
+    if (settings::run_mode == RunMode::EIGENVALUE) {
+      std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::K_COLLISION, 0), n);
       std::cout << " k-effective (Collision)     = "
         << mean << " +/- " << t_n1 * stdev << '\n';
-      std::tie(mean, stdev) = mean_stdev(&gt(K_TRACKLENGTH, 0), n);
+      std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::K_TRACKLENGTH, 0), n);
       std::cout << " k-effective (Track-length)  = "
         << mean << " +/- " << t_n1 * stdev << '\n';
-      std::tie(mean, stdev) = mean_stdev(&gt(K_ABSORPTION, 0), n);
+      std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::K_ABSORPTION, 0), n);
       std::cout << " k-effective (Absorption)    = "
         << mean << " +/- " << t_n1 * stdev << '\n';
       if (n > 3) {
@@ -577,23 +577,23 @@ void print_results()
           << k_combined[0] << " +/- " << t_n3 * k_combined[1] << '\n';
       }
     }
-    std::tie(mean, stdev) = mean_stdev(&gt(LEAKAGE, 0), n);
+    std::tie(mean, stdev) = mean_stdev(&gt(GlobalTally::LEAKAGE, 0), n);
     std::cout << " Leakage Fraction            = "
       << mean << " +/- " << t_n1 * stdev << '\n';
   } else {
     if (mpi::master) warning("Could not compute uncertainties -- only one "
       "active batch simulated!");
 
-    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+    if (settings::run_mode == RunMode::EIGENVALUE) {
       std::cout << " k-effective (Collision)    = "
-        << gt(K_COLLISION, RESULT_SUM) / n << '\n';
+        << gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n << '\n';
       std::cout << " k-effective (Track-length) = "
-        << gt(K_TRACKLENGTH, RESULT_SUM) / n << '\n';
+        << gt(GlobalTally::K_TRACKLENGTH, TallyResult::SUM) / n << '\n';
       std::cout << " k-effective (Absorption)   = "
-        << gt(K_ABSORPTION, RESULT_SUM) / n << '\n';
+        << gt(GlobalTally::K_ABSORPTION, TallyResult::SUM) / n << '\n';
     }
     std::cout << " Leakage Fraction           = "
-      << gt(LEAKAGE, RESULT_SUM) / n << '\n';
+      << gt(GlobalTally::LEAKAGE, TallyResult::SUM) / n << '\n';
   }
   std::cout << '\n';
 
@@ -659,16 +659,16 @@ write_tallies()
     if (tally.deriv_ != C_NONE) {
       const auto& deriv {model::tally_derivs[tally.deriv_]};
       switch (deriv.variable) {
-      case DIFF_DENSITY:
+      case WithRespectTo::DENSITY:
         tallies_out << " Density derivative Material "
           << std::to_string(deriv.diff_material) << "\n";
         break;
-      case DIFF_NUCLIDE_DENSITY:
+      case WithRespectTo::NUCLIDE_DENSITY:
         tallies_out << " Nuclide density derivative Material "
           << std::to_string(deriv.diff_material) << "  Nuclide "
           << data::nuclides[deriv.diff_nuclide]->name_ << "\n";
         break;
-      case DIFF_TEMPERATURE:
+      case WithRespectTo::TEMPERATURE:
         tallies_out << " Temperature derivative Material "
           << std::to_string(deriv.diff_material) << "\n";
         break;
diff --git a/src/particle.cpp b/src/particle.cpp
index 25f595082..724821824 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -169,7 +169,7 @@ Particle::transport()
     r_last_ = this->r();
 
     // Reset event variables
-    event_ = EVENT_KILL;
+    event_ = TallyEvent::KILL;
     event_nuclide_ = NUCLIDE_NONE;
     event_mt_ = REACTION_NONE;
 
@@ -245,7 +245,7 @@ Particle::transport()
     }
 
     // Score track-length estimate of k-eff
-    if (settings::run_mode == RUN_MODE_EIGENVALUE &&
+    if (settings::run_mode == RunMode::EIGENVALUE &&
         type_ == Particle::Type::neutron) {
       global_tally_tracklength += wgt_ * distance * macro_xs_.nu_fission;
     }
@@ -274,11 +274,11 @@ Particle::transport()
           boundary.lattice_translation[2] != 0) {
         // Particle crosses lattice boundary
         cross_lattice(this, boundary);
-        event_ = EVENT_LATTICE;
+        event_ = TallyEvent::LATTICE;
       } else {
         // Particle crosses surface
         this->cross_surface();
-        event_ = EVENT_SURFACE;
+        event_ = TallyEvent::SURFACE;
       }
       // Score cell to cell partial currents
       if (!model::active_surface_tallies.empty()) {
@@ -289,7 +289,7 @@ Particle::transport()
       // PARTICLE HAS COLLISION
 
       // Score collision estimate of keff
-      if (settings::run_mode == RUN_MODE_EIGENVALUE &&
+      if (settings::run_mode == RunMode::EIGENVALUE &&
           type_ == Particle::Type::neutron) {
         global_tally_collision += wgt_ * macro_xs_.nu_fission
           / macro_xs_.total;
@@ -400,7 +400,7 @@ Particle::cross_surface()
     write_message("    Crossing surface " + std::to_string(surf->id_));
   }
 
-  if (surf->bc_ == BC_VACUUM && (settings::run_mode != RUN_MODE_PLOTTING)) {
+  if (surf->bc_ == Surface::Bc::VACUUM && (settings::run_mode != RunMode::PLOTTING)) {
     // =======================================================================
     // PARTICLE LEAKS OUT OF PROBLEM
 
@@ -428,8 +428,8 @@ Particle::cross_surface()
     }
     return;
 
-  } else if ((surf->bc_ == BC_REFLECT || surf->bc_ == BC_WHITE)
-                                    && (settings::run_mode != RUN_MODE_PLOTTING)) {
+  } else if ((surf->bc_ == Surface::Bc::REFLECT || surf->bc_ == Surface::Bc::WHITE)
+                                    && (settings::run_mode != RunMode::PLOTTING)) {
     // =======================================================================
     // PARTICLE REFLECTS FROM SURFACE
 
@@ -459,7 +459,7 @@ Particle::cross_surface()
       this->r() = r;
     }
 
-    Direction u = (surf->bc_ == BC_REFLECT) ?
+    Direction u = (surf->bc_ == Surface::Bc::REFLECT) ?
       surf->reflect(this->r(), this->u()) :
       surf->diffuse_reflect(this->r(), this->u(), this->current_seed());
 
@@ -492,7 +492,7 @@ Particle::cross_surface()
     }
     return;
 
-  } else if (surf->bc_ == BC_PERIODIC && settings::run_mode != RUN_MODE_PLOTTING) {
+  } else if (surf->bc_ == Surface::Bc::PERIODIC && settings::run_mode != RunMode::PLOTTING) {
     // =======================================================================
     // PERIODIC BOUNDARY
 
@@ -579,7 +579,7 @@ Particle::cross_surface()
   n_coord_ = 1;
   bool found = find_cell(this, false);
 
-  if (settings::run_mode != RUN_MODE_PLOTTING && (!found)) {
+  if (settings::run_mode != RunMode::PLOTTING && (!found)) {
     // If a cell is still not found, there are two possible causes: 1) there is
     // a void in the model, and 2) the particle hit a surface at a tangent. If
     // the particle is really traveling tangent to a surface, if we move it
@@ -628,7 +628,7 @@ void
 Particle::write_restart() const
 {
   // Dont write another restart file if in particle restart mode
-  if (settings::run_mode == RUN_MODE_PARTICLE) return;
+  if (settings::run_mode == RunMode::PARTICLE) return;
 
   // Set up file name
   std::stringstream filename;
@@ -654,13 +654,13 @@ Particle::write_restart() const
     write_dataset(file_id, "current_generation", simulation::current_gen);
     write_dataset(file_id, "n_particles", settings::n_particles);
     switch (settings::run_mode) {
-      case RUN_MODE_FIXEDSOURCE:
+      case RunMode::FIXEDSOURCE:
         write_dataset(file_id, "run_mode", "fixed source");
         break;
-      case RUN_MODE_EIGENVALUE:
+      case RunMode::EIGENVALUE:
         write_dataset(file_id, "run_mode", "eigenvalue");
         break;
-      case RUN_MODE_PARTICLE:
+      case RunMode::PARTICLE:
         write_dataset(file_id, "run_mode", "particle restart");
         break;
     }
diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp
index 053193207..6b206c42c 100644
--- a/src/particle_restart.cpp
+++ b/src/particle_restart.cpp
@@ -19,7 +19,7 @@
 
 namespace openmc {
 
-void read_particle_restart(Particle& p, int& previous_run_mode)
+void read_particle_restart(Particle& p, RunMode& previous_run_mode)
 {
   // Write meessage
   write_message("Loading particle restart file " +
@@ -36,9 +36,9 @@ void read_particle_restart(Particle& p, int& previous_run_mode)
   std::string mode;
   read_dataset(file_id, "run_mode", mode);
   if (mode == "eigenvalue") {
-    previous_run_mode = RUN_MODE_EIGENVALUE;
+    previous_run_mode = RunMode::EIGENVALUE;
   } else if (mode == "fixed source") {
-    previous_run_mode = RUN_MODE_FIXEDSOURCE;
+    previous_run_mode = RunMode::FIXEDSOURCE;
   }
   read_dataset(file_id, "id", p.id_);
   int type;
@@ -76,7 +76,7 @@ void run_particle_restart()
   Particle p;
 
   // Read in the restart information
-  int previous_run_mode;
+  RunMode previous_run_mode;
   read_particle_restart(p, previous_run_mode);
 
   // write track if that was requested on command line
@@ -88,15 +88,15 @@ void run_particle_restart()
   // Compute random number seed
   int64_t particle_seed;
   switch (previous_run_mode) {
-  case RUN_MODE_EIGENVALUE:
+  case RunMode::EIGENVALUE:
     particle_seed = (simulation::total_gen + overall_generation() - 1)*settings::n_particles + p.id_;
     break;
-  case RUN_MODE_FIXEDSOURCE:
+  case RunMode::FIXEDSOURCE:
     particle_seed = p.id_;
     break;
   default:
     throw std::runtime_error{"Unexpected run mode: " +
-      std::to_string(previous_run_mode)};
+      std::to_string(static_cast<int>(previous_run_mode))};
   }
   init_particle_seeds(particle_seed, p.seeds_);
 
diff --git a/src/photon.cpp b/src/photon.cpp
index de90dee15..7af8687df 100644
--- a/src/photon.cpp
+++ b/src/photon.cpp
@@ -207,7 +207,7 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
   // Calculate total pair production
   pair_production_total_ = pair_production_nuclear_ + pair_production_electron_;
 
-  if (settings::electron_treatment == ELECTRON_TTB) {
+  if (settings::electron_treatment == ElectronTreatment::TTB) {
     // Read bremsstrahlung scaled DCS
     rgroup = open_group(group, "bremsstrahlung");
     read_dataset(rgroup, "dcs", dcs_);
diff --git a/src/physics.cpp b/src/physics.cpp
index 4b80d1039..fa31214c5 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -62,7 +62,7 @@ void collision(Particle* p)
   // Display information about collision
   if (settings::verbosity >= 10 || simulation::trace) {
     std::stringstream msg;
-    if (p->event_ == EVENT_KILL) {
+    if (p->event_ == TallyEvent::KILL) {
       msg << "    Killed. Energy = " << p->E_ << " eV.";
     } else if (p->type_ == Particle::Type::neutron) {
       msg << "    " << reaction_name(p->event_mt_) << " with " <<
@@ -94,9 +94,9 @@ void sample_neutron_reaction(Particle* p)
 
   if (nuc->fissionable_) {
     Reaction* rx = sample_fission(i_nuclide, p);
-    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+    if (settings::run_mode == RunMode::EIGENVALUE) {
       create_fission_sites(p, i_nuclide, rx, simulation::fission_bank);
-    } else if (settings::run_mode == RUN_MODE_FIXEDSOURCE &&
+    } else if (settings::run_mode == RunMode::FIXEDSOURCE &&
       settings::create_fission_neutrons) {
       create_fission_sites(p, i_nuclide, rx, simulation::secondary_bank);
 
@@ -230,7 +230,7 @@ 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_ = EVENT_SCATTER;
+    p->event_ = TallyEvent::SCATTER;
     p->event_mt_ = COHERENT;
     return;
   }
@@ -273,7 +273,7 @@ void sample_photon_reaction(Particle* p)
     phi += PI;
     p->E_ = alpha_out*MASS_ELECTRON_EV;
     p->u() = rotate_angle(p->u(), mu, &phi, p->current_seed());
-    p->event_ = EVENT_SCATTER;
+    p->event_ = TallyEvent::SCATTER;
     p->event_mt_ = INCOHERENT;
     return;
   }
@@ -325,7 +325,7 @@ void sample_photon_reaction(Particle* p)
         // Allow electrons to fill orbital and produce auger electrons
         // and fluorescent photons
         element.atomic_relaxation(shell, *p);
-        p->event_ = EVENT_ABSORB;
+        p->event_ = TallyEvent::ABSORB;
         p->event_mt_ = 533 + shell.index_subshell;
         p->alive_ = false;
         p->E_ = 0.0;
@@ -351,7 +351,7 @@ void sample_photon_reaction(Particle* p)
     u = rotate_angle(p->u(), mu_positron, nullptr, p->current_seed());
     p->create_secondary(u, E_positron, Particle::Type::positron);
 
-    p->event_ = EVENT_ABSORB;
+    p->event_ = TallyEvent::ABSORB;
     p->event_mt_ = PAIR_PROD;
     p->alive_ = false;
     p->E_ = 0.0;
@@ -362,21 +362,21 @@ void sample_electron_reaction(Particle* p)
 {
   // TODO: create reaction types
 
-  if (settings::electron_treatment == ELECTRON_TTB) {
+  if (settings::electron_treatment == ElectronTreatment::TTB) {
     double E_lost;
     thick_target_bremsstrahlung(*p, &E_lost);
   }
 
   p->E_ = 0.0;
   p->alive_ = false;
-  p->event_ = EVENT_ABSORB;
+  p->event_ = TallyEvent::ABSORB;
 }
 
 void sample_positron_reaction(Particle* p)
 {
   // TODO: create reaction types
 
-  if (settings::electron_treatment == ELECTRON_TTB) {
+  if (settings::electron_treatment == ElectronTreatment::TTB) {
     double E_lost;
     thick_target_bremsstrahlung(*p, &E_lost);
   }
@@ -395,7 +395,7 @@ void sample_positron_reaction(Particle* p)
 
   p->E_ = 0.0;
   p->alive_ = false;
-  p->event_ = EVENT_ABSORB;
+  p->event_ = TallyEvent::ABSORB;
 }
 
 int sample_nuclide(Particle* p)
@@ -547,7 +547,7 @@ void absorption(Particle* p, int i_nuclide)
     p->wgt_last_ = p->wgt_;
 
     // Score implicit absorption estimate of keff
-    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+    if (settings::run_mode == RunMode::EIGENVALUE) {
       global_tally_absorption += p->wgt_absorb_ * p->neutron_xs_[
         i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption;
     }
@@ -556,13 +556,13 @@ void absorption(Particle* p, int i_nuclide)
     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 == RUN_MODE_EIGENVALUE) {
+      if (settings::run_mode == RunMode::EIGENVALUE) {
         global_tally_absorption += p->wgt_ * p->neutron_xs_[
           i_nuclide].nu_fission / p->neutron_xs_[i_nuclide].absorption;
       }
 
       p->alive_ = false;
-      p->event_ = EVENT_ABSORB;
+      p->event_ = TallyEvent::ABSORB;
       p->event_mt_ = N_DISAPPEAR;
     }
   }
@@ -648,7 +648,7 @@ void scatter(Particle* p, int i_nuclide)
   }
 
   // Set event component
-  p->event_ = EVENT_SCATTER;
+  p->event_ = TallyEvent::SCATTER;
 
   // Sample new outgoing angle for isotropic-in-lab scattering
   const auto& mat {model::materials[p->material_]};
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index f52e6c5bf..91a02b87c 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -47,9 +47,9 @@ sample_reaction(Particle* p)
   // absorption (including fission)
 
   if (model::materials[p->material_]->fissionable_) {
-    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+    if (settings::run_mode == RunMode::EIGENVALUE) {
       create_fission_sites(p, simulation::fission_bank);
-    } else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
+    } else if ((settings::run_mode == RunMode::FIXEDSOURCE) &&
                (settings::create_fission_neutrons)) {
       create_fission_sites(p, simulation::secondary_bank);
     }
@@ -87,7 +87,7 @@ scatter(Particle* p)
   p->E_ = data::mg.energy_bin_avg_[p->g_];
 
   // Set event component
-  p->event_ = EVENT_SCATTER;
+  p->event_ = TallyEvent::SCATTER;
 }
 
 void
@@ -185,7 +185,7 @@ absorption(Particle* p)
       global_tally_absorption += p->wgt_ * p->macro_xs_.nu_fission /
            p->macro_xs_.absorption;
       p->alive_ = false;
-      p->event_ = EVENT_ABSORB;
+      p->event_ = TallyEvent::ABSORB;
     }
 
   }
diff --git a/src/plot.cpp b/src/plot.cpp
index 31b771a47..5aa6a61b8 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -44,7 +44,7 @@ IdData::set_value(size_t y, size_t x, const Particle& p, int level) {
   if (p.material_ == MATERIAL_VOID) {
     data_(y,x,1) = MATERIAL_VOID;
     return;
-  } else if (c->type_ != FILL_UNIVERSE) {
+  } else if (c->type_ != Fill::UNIVERSE) {
     Material* m = model::materials[p.material_].get();
     data_(y,x,1) = m->id_;
   }
@@ -62,7 +62,7 @@ void
 PropertyData::set_value(size_t y, size_t x, const Particle& p, int level) {
   Cell* c = model::cells[p.coord_[level].cell].get();
   data_(y,x,0) = (p.sqrtkT_ * p.sqrtkT_) / K_BOLTZMANN;
-  if (c->type_ != FILL_UNIVERSE && p.material_ != MATERIAL_VOID) {
+  if (c->type_ != Fill::UNIVERSE && p.material_ != MATERIAL_VOID) {
     Material* m = model::materials[p.material_].get();
     data_(y,x,1) = m->density_gpcc_;
   }
@@ -678,7 +678,7 @@ void Plot::set_overlap_color(pugi::xml_node plot_node) {
 
   // make sure we allocate the vector for counting overlap checks if
   // they're going to be plotted
-  if (color_overlaps_ && settings::run_mode == RUN_MODE_PLOTTING) {
+  if (color_overlaps_ && settings::run_mode == RunMode::PLOTTING) {
     settings::check_overlaps = true;
     model::overlap_check_count.resize(model::cells.size(), 0);
   }
diff --git a/src/scattdata.cpp b/src/scattdata.cpp
index c61d0249c..11e5c4c4c 100644
--- a/src/scattdata.cpp
+++ b/src/scattdata.cpp
@@ -183,7 +183,7 @@ ScattData::sample_energy(int gin, int& gout, int& i_gout, uint64_t* seed)
 //==============================================================================
 
 double
-ScattData::get_xs(int xstype, int gin, const int* gout, const double* mu)
+ScattData::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu)
 {
   // Set the outgoing group offset index as needed
   int i_gout = 0;
@@ -198,10 +198,10 @@ ScattData::get_xs(int xstype, int gin, const int* gout, const double* mu)
 
   double val = scattxs[gin];
   switch(xstype) {
-  case MG_GET_XS_SCATTER:
+  case MgxsType::SCATTER:
     if (gout != nullptr) val *= energy[gin][i_gout];
     break;
-  case MG_GET_XS_SCATTER_MULT:
+  case MgxsType::SCATTER_MULT:
     if (gout != nullptr) {
       val *= energy[gin][i_gout] / mult[gin][i_gout];
     } else {
@@ -209,7 +209,7 @@ ScattData::get_xs(int xstype, int gin, const int* gout, const double* mu)
                                 energy[gin].begin(), 0.0);
     }
     break;
-  case MG_GET_XS_SCATTER_FMU_MULT:
+  case MgxsType::SCATTER_FMU_MULT:
     if ((gout != nullptr) && (mu != nullptr)) {
       val *= energy[gin][i_gout] * calc_f(gin, *gout, *mu);
     } else {
@@ -218,7 +218,7 @@ ScattData::get_xs(int xstype, int gin, const int* gout, const double* mu)
       fatal_error("Invalid call to get_xs");
     }
     break;
-  case MG_GET_XS_SCATTER_FMU:
+  case MgxsType::SCATTER_FMU:
     if ((gout != nullptr) && (mu != nullptr)) {
       val *= energy[gin][i_gout] * calc_f(gin, *gout, *mu) / mult[gin][i_gout];
     } else {
diff --git a/src/settings.cpp b/src/settings.cpp
index 59d2479be..ceee574f9 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -80,7 +80,7 @@ int32_t n_inactive {0};
 int32_t gen_per_batch {1};
 int64_t n_particles {-1};
 
-int electron_treatment {ELECTRON_TTB};
+ElectronTreatment electron_treatment {ElectronTreatment::TTB};
 std::array<double, 4> energy_cutoff {0.0, 1000.0, 0.0, 0.0};
 int legendre_to_tabular_points {C_NONE};
 int max_order {0};
@@ -90,10 +90,10 @@ ResScatMethod res_scat_method {ResScatMethod::rvs};
 double res_scat_energy_min {0.01};
 double res_scat_energy_max {1000.0};
 std::vector<std::string> res_scat_nuclides;
-int run_mode {-1};
+RunMode run_mode {RunMode::UNSET};
 std::unordered_set<int> sourcepoint_batch;
 std::unordered_set<int> statepoint_batch;
-int temperature_method {TEMPERATURE_NEAREST};
+TemperatureInterpolationType temperature_method {TemperatureInterpolationType::NEAREST};
 double temperature_tolerance {10.0};
 double temperature_default {293.6};
 std::array<double, 2> temperature_range {0.0, 0.0};
@@ -134,7 +134,7 @@ void get_run_parameters(pugi::xml_node node_base)
   if (!trigger_on) n_max_batches = n_batches;
 
   // Get number of inactive batches
-  if (run_mode == RUN_MODE_EIGENVALUE) {
+  if (run_mode == RunMode::EIGENVALUE) {
     if (check_for_node(node_base, "inactive")) {
       n_inactive = std::stoi(get_node_value(node_base, "inactive"));
     }
@@ -184,7 +184,7 @@ void read_settings_xml()
   // Check if settings.xml exists
   std::string filename = path_input + "settings.xml";
   if (!file_exists(filename)) {
-    if (run_mode != RUN_MODE_PLOTTING) {
+    if (run_mode != RunMode::PLOTTING) {
       std::stringstream msg;
       msg << "Settings XML file '" << filename << "' does not exist! In order "
         "to run OpenMC, you first need a set of input files; at a minimum, this "
@@ -292,19 +292,19 @@ void read_settings_xml()
 
   // Check run mode if it hasn't been set from the command line
   xml_node node_mode;
-  if (run_mode == C_NONE) {
+  if (run_mode == RunMode::UNSET) {
     if (check_for_node(root, "run_mode")) {
       std::string temp_str = get_node_value(root, "run_mode", true, true);
       if (temp_str == "eigenvalue") {
-        run_mode = RUN_MODE_EIGENVALUE;
+        run_mode = RunMode::EIGENVALUE;
       } else if (temp_str == "fixed source") {
-        run_mode = RUN_MODE_FIXEDSOURCE;
+        run_mode = RunMode::FIXEDSOURCE;
       } else if (temp_str == "plot") {
-        run_mode = RUN_MODE_PLOTTING;
+        run_mode = RunMode::PLOTTING;
       } else if (temp_str == "particle restart") {
-        run_mode = RUN_MODE_PARTICLE;
+        run_mode = RunMode::PARTICLE;
       } else if (temp_str == "volume") {
-        run_mode = RUN_MODE_VOLUME;
+        run_mode = RunMode::VOLUME;
       } else {
         fatal_error("Unrecognized run mode: " + temp_str);
       }
@@ -317,11 +317,11 @@ void read_settings_xml()
       // Make sure that either eigenvalue or fixed source was specified
       node_mode = root.child("eigenvalue");
       if (node_mode) {
-        run_mode = RUN_MODE_EIGENVALUE;
+        run_mode = RunMode::EIGENVALUE;
       } else {
         node_mode = root.child("fixed_source");
         if (node_mode) {
-          run_mode = RUN_MODE_FIXEDSOURCE;
+          run_mode = RunMode::FIXEDSOURCE;
         } else {
           fatal_error("<eigenvalue> or <fixed_source> not specified.");
         }
@@ -329,7 +329,7 @@ void read_settings_xml()
     }
   }
 
-  if (run_mode == RUN_MODE_EIGENVALUE || run_mode == RUN_MODE_FIXEDSOURCE) {
+  if (run_mode == RunMode::EIGENVALUE || run_mode == RunMode::FIXEDSOURCE) {
     // Read run parameters
     get_run_parameters(node_mode);
 
@@ -353,9 +353,9 @@ void read_settings_xml()
   if (check_for_node(root, "electron_treatment")) {
     auto temp_str = get_node_value(root, "electron_treatment", true, true);
     if (temp_str == "led") {
-      electron_treatment = ELECTRON_LED;
+      electron_treatment = ElectronTreatment::LED;
     } else if (temp_str == "ttb") {
-      electron_treatment = ELECTRON_TTB;
+      electron_treatment = ElectronTreatment::TTB;
     } else {
       fatal_error("Unrecognized electron treatment: " + temp_str + ".");
     }
@@ -731,9 +731,9 @@ void read_settings_xml()
   if (check_for_node(root, "temperature_method")) {
     auto temp = get_node_value(root, "temperature_method", true, true);
     if (temp == "nearest") {
-      temperature_method = TEMPERATURE_NEAREST;
+      temperature_method = TemperatureInterpolationType::NEAREST;
     } else if (temp == "interpolation") {
-      temperature_method = TEMPERATURE_INTERPOLATION;
+      temperature_method = TemperatureInterpolationType::INTERPOLATION;
     } else {
       fatal_error("Unknown temperature method: " + temp);
     }
@@ -772,7 +772,7 @@ void read_settings_xml()
   }
 
   // Check whether create fission sites
-  if (run_mode == RUN_MODE_FIXEDSOURCE) {
+  if (run_mode == RunMode::FIXEDSOURCE) {
     if (check_for_node(root, "create_fission_neutrons")) {
       create_fission_neutrons = get_node_value_bool(root, "create_fission_neutrons");
     }
diff --git a/src/simulation.cpp b/src/simulation.cpp
index eaed2badf..c7c096b58 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -100,9 +100,9 @@ int openmc_simulation_init()
 
   // Display header
   if (mpi::master) {
-    if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
+    if (settings::run_mode == RunMode::FIXEDSOURCE) {
       header("FIXED SOURCE TRANSPORT SIMULATION", 3);
-    } else if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+    } else if (settings::run_mode == RunMode::EIGENVALUE) {
       header("K EIGENVALUE SIMULATION", 3);
       if (settings::verbosity >= 7) print_columns();
     }
@@ -278,7 +278,7 @@ void allocate_banks()
   // Allocate source bank
   simulation::source_bank.resize(simulation::work_per_rank);
 
-  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+  if (settings::run_mode == RunMode::EIGENVALUE) {
 #ifdef _OPENMP
     // If OpenMP is being used, each thread needs its own private fission
     // bank. Since the private fission banks need to be combined at the end of
@@ -306,7 +306,7 @@ void initialize_batch()
   // Increment current batch
   ++simulation::current_batch;
 
-  if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
+  if (settings::run_mode == RunMode::FIXEDSOURCE) {
     int b = simulation::current_batch;
     write_message("Simulating batch " + std::to_string(b), 6);
   }
@@ -353,7 +353,7 @@ void finalize_batch()
     simulation::n_realizations = 0;
   }
 
-  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+  if (settings::run_mode == RunMode::EIGENVALUE) {
     // Write batch output
     if (mpi::master && settings::verbosity >= 7) print_batch_keff();
   }
@@ -397,7 +397,7 @@ void finalize_batch()
 
 void initialize_generation()
 {
-  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+  if (settings::run_mode == RunMode::EIGENVALUE) {
     // Clear out the fission bank
     simulation::fission_bank.clear();
 
@@ -406,7 +406,7 @@ void initialize_generation()
 
     // Store current value of tracklength k
     simulation::keff_generation = simulation::global_tallies(
-      K_TRACKLENGTH, RESULT_VALUE);
+      GlobalTally::K_TRACKLENGTH, TallyResult::VALUE);
   }
 }
 
@@ -419,16 +419,16 @@ void finalize_generation()
   {
     #pragma omp critical(increment_global_tallies)
     {
-      if (settings::run_mode == RUN_MODE_EIGENVALUE) {
-        gt(K_COLLISION, RESULT_VALUE) += global_tally_collision;
-        gt(K_ABSORPTION, RESULT_VALUE) += global_tally_absorption;
-        gt(K_TRACKLENGTH, RESULT_VALUE) += global_tally_tracklength;
+      if (settings::run_mode == RunMode::EIGENVALUE) {
+        gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += global_tally_collision;
+        gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) += global_tally_absorption;
+        gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) += global_tally_tracklength;
       }
-      gt(LEAKAGE, RESULT_VALUE) += global_tally_leakage;
+      gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage;
     }
 
     // reset threadprivate tallies
-    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+    if (settings::run_mode == RunMode::EIGENVALUE) {
       global_tally_collision = 0.0;
       global_tally_absorption = 0.0;
       global_tally_tracklength = 0.0;
@@ -437,7 +437,7 @@ void finalize_generation()
   }
 
 
-  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+  if (settings::run_mode == RunMode::EIGENVALUE) {
 #ifdef _OPENMP
     // Join the fission bank from each thread into one global fission bank
     join_bank_from_threads();
@@ -460,7 +460,7 @@ void finalize_generation()
       }
     }
 
-  } else if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
+  } else if (settings::run_mode == RunMode::FIXEDSOURCE) {
     // For fixed-source mode, we need to sample the external source
     fill_source_bank_fixedsource();
   }
diff --git a/src/state_point.cpp b/src/state_point.cpp
index 4c159c167..5c99c227c 100644
--- a/src/state_point.cpp
+++ b/src/state_point.cpp
@@ -83,10 +83,10 @@ openmc_statepoint_write(const char* filename, bool* write_source)
     write_dataset(file_id, "energy_mode", settings::run_CE ?
       "continuous-energy" : "multi-group");
     switch (settings::run_mode) {
-    case RUN_MODE_FIXEDSOURCE:
+    case RunMode::FIXEDSOURCE:
       write_dataset(file_id, "run_mode", "fixed source");
       break;
-    case RUN_MODE_EIGENVALUE:
+    case RunMode::EIGENVALUE:
       write_dataset(file_id, "run_mode", "eigenvalue");
       break;
     }
@@ -101,7 +101,7 @@ openmc_statepoint_write(const char* filename, bool* write_source)
     write_attribute(file_id, "source_present", write_source_);
 
     // Write out information for eigenvalue run
-    if (settings::run_mode == RUN_MODE_EIGENVALUE)
+    if (settings::run_mode == RunMode::EIGENVALUE)
       write_eigenvalue_hdf5(file_id);
 
     hid_t tallies_group = create_group(file_id, "tallies");
@@ -116,13 +116,13 @@ openmc_statepoint_write(const char* filename, bool* write_source)
         hid_t deriv_group = create_group(derivs_group,
           "derivative " + std::to_string(deriv.id));
         write_dataset(deriv_group, "material", deriv.diff_material);
-        if (deriv.variable == DIFF_DENSITY) {
+        if (deriv.variable == WithRespectTo::DENSITY) {
           write_dataset(deriv_group, "independent variable", "density");
-        } else if (deriv.variable == DIFF_NUCLIDE_DENSITY) {
+        } else if (deriv.variable == WithRespectTo::NUCLIDE_DENSITY) {
           write_dataset(deriv_group, "independent variable", "nuclide_density");
           write_dataset(deriv_group, "nuclide",
             data::nuclides[deriv.diff_nuclide]->name_);
-        } else if (deriv.variable == DIFF_TEMPERATURE) {
+        } else if (deriv.variable == WithRespectTo::TEMPERATURE) {
           write_dataset(deriv_group, "independent variable", "temperature");
         } else {
           fatal_error("Independent variable for derivative "
@@ -179,11 +179,11 @@ openmc_statepoint_write(const char* filename, bool* write_source)
           continue;
         }
 
-        if (tally->estimator_ == ESTIMATOR_ANALOG) {
+        if (tally->estimator_ == TallyEstimator::ANALOG) {
           write_dataset(tally_group, "estimator", "analog");
-        } else if (tally->estimator_ == ESTIMATOR_TRACKLENGTH) {
+        } else if (tally->estimator_ == TallyEstimator::TRACKLENGTH) {
           write_dataset(tally_group, "estimator", "tracklength");
-        } else if (tally->estimator_ == ESTIMATOR_COLLISION) {
+        } else if (tally->estimator_ == TallyEstimator::COLLISION) {
           write_dataset(tally_group, "estimator", "collision");
         }
 
@@ -275,11 +275,11 @@ openmc_statepoint_write(const char* filename, bool* write_source)
     write_dataset(runtime_group, "simulation", time_inactive.elapsed()
       + time_active.elapsed());
     write_dataset(runtime_group, "transport", time_transport.elapsed());
-    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+    if (settings::run_mode == RunMode::EIGENVALUE) {
       write_dataset(runtime_group, "inactive batches", time_inactive.elapsed());
     }
     write_dataset(runtime_group, "active batches", time_active.elapsed());
-    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+    if (settings::run_mode == RunMode::EIGENVALUE) {
       write_dataset(runtime_group, "synchronizing fission bank", time_bank.elapsed());
       write_dataset(runtime_group, "sampling source sites", time_bank_sample.elapsed());
       write_dataset(runtime_group, "SEND-RECV source sites", time_bank_sendrecv.elapsed());
@@ -363,9 +363,9 @@ void load_state_point()
   // Read and overwrite run information except number of batches
   read_dataset(file_id, "run_mode", word);
   if (word == "fixed source") {
-    settings::run_mode = RUN_MODE_FIXEDSOURCE;
+    settings::run_mode = RunMode::FIXEDSOURCE;
   } else if (word == "eigenvalue") {
-    settings::run_mode = RUN_MODE_EIGENVALUE;
+    settings::run_mode = RunMode::EIGENVALUE;
   }
   read_attribute(file_id, "photon_transport", settings::photon_transport);
   read_dataset(file_id, "n_particles", settings::n_particles);
@@ -388,7 +388,7 @@ void load_state_point()
   }
 
   // Read information specific to eigenvalue run
-  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+  if (settings::run_mode == RunMode::EIGENVALUE) {
     read_dataset(file_id, "n_inactive", temp);
     read_eigenvalue_hdf5(file_id);
 
@@ -455,7 +455,7 @@ void load_state_point()
   }
 
   // Read source if in eigenvalue mode
-  if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+  if (settings::run_mode == RunMode::EIGENVALUE) {
 
     // Check if source was written out separately
     if (!source_present) {
@@ -723,7 +723,7 @@ void write_tally_results_nr(hid_t file_id)
 
     // Get view of accumulated tally values
     auto values_view = xt::view(t->results_, xt::all(), xt::all(),
-      xt::range(RESULT_SUM, RESULT_SUM_SQ + 1));
+      xt::range(static_cast<int>(TallyResult::SUM), static_cast<int>(TallyResult::SUM_SQ) + 1));
 
     // Make copy of tally values in contiguous array
     xt::xtensor<double, 2> values = values_view;
@@ -750,7 +750,7 @@ void write_tally_results_nr(hid_t file_id)
       // Put in temporary tally result
       xt::xtensor<double, 3> results_copy = xt::zeros_like(t->results_);
       auto copy_view = xt::view(results_copy, xt::all(), xt::all(),
-        xt::range(RESULT_SUM, RESULT_SUM_SQ + 1));
+        xt::range(static_cast<int>(TallyResult::SUM), static_cast<int>(TallyResult::SUM_SQ) + 1));
       copy_view = values;
 
       // Write reduced tally results to file
diff --git a/src/surface.cpp b/src/surface.cpp
index abdaacbed..0e25e4ff9 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -16,15 +16,6 @@
 
 namespace openmc {
 
-//==============================================================================
-// Module constant definitions
-//==============================================================================
-
-extern "C" const int BC_TRANSMIT {0};
-extern "C" const int BC_VACUUM {1};
-extern "C" const int BC_REFLECT {2};
-extern "C" const int BC_PERIODIC {3};
-extern "C" const int BC_WHITE {4};
 //==============================================================================
 // Global variables
 //==============================================================================
@@ -141,18 +132,18 @@ Surface::Surface(pugi::xml_node surf_node)
     std::string surf_bc = get_node_value(surf_node, "boundary", true, true);
 
     if (surf_bc == "transmission" || surf_bc == "transmit" ||surf_bc.empty()) {
-      bc_ = BC_TRANSMIT;
+      bc_ = Bc::TRANSMIT;
 
     } else if (surf_bc == "vacuum") {
-      bc_ = BC_VACUUM;
+      bc_ = Bc::VACUUM;
 
     } else if (surf_bc == "reflective" || surf_bc == "reflect"
                || surf_bc == "reflecting") {
-      bc_ = BC_REFLECT;
+      bc_ = Bc::REFLECT;
     } else if (surf_bc == "white") {
-      bc_ = BC_WHITE;
+      bc_ = Bc::WHITE;
     } else if (surf_bc == "periodic") {
-      bc_ = BC_PERIODIC;
+      bc_ = Bc::PERIODIC;
     } else {
       std::stringstream err_msg;
       err_msg << "Unknown boundary condition \"" << surf_bc
@@ -161,7 +152,7 @@ Surface::Surface(pugi::xml_node surf_node)
     }
 
   } else {
-    bc_ = BC_TRANSMIT;
+    bc_ = Bc::TRANSMIT;
   }
 
 }
@@ -229,19 +220,19 @@ CSGSurface::to_hdf5(hid_t group_id) const
   hid_t surf_group = create_group(group_id, group_name);
 
   switch(bc_) {
-    case BC_TRANSMIT :
+    case Bc::TRANSMIT :
       write_string(surf_group, "boundary_type", "transmission", false);
       break;
-    case BC_VACUUM :
+    case Bc::VACUUM :
       write_string(surf_group, "boundary_type", "vacuum", false);
       break;
-    case BC_REFLECT :
+    case Bc::REFLECT :
       write_string(surf_group, "boundary_type", "reflective", false);
       break;
-    case BC_WHITE :
+    case Bc::WHITE :
       write_string(surf_group, "boundary_type", "white", false);
       break;
-    case BC_PERIODIC :
+    case Bc::PERIODIC :
       write_string(surf_group, "boundary_type", "periodic", false);
       break;
   }
@@ -1204,7 +1195,7 @@ void read_surfaces(pugi::xml_node node)
          zmin {INFTY}, zmax {-INFTY};
   int i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax;
   for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) {
-    if (model::surfaces[i_surf]->bc_ == BC_PERIODIC) {
+    if (model::surfaces[i_surf]->bc_ == Surface::Bc::PERIODIC) {
       // Downcast to the PeriodicSurface type.
       Surface* surf_base = model::surfaces[i_surf].get();
       auto surf = dynamic_cast<PeriodicSurface*>(surf_base);
@@ -1249,7 +1240,7 @@ void read_surfaces(pugi::xml_node node)
 
   // Set i_periodic for periodic BC surfaces.
   for (int i_surf = 0; i_surf < model::surfaces.size(); i_surf++) {
-    if (model::surfaces[i_surf]->bc_ == BC_PERIODIC) {
+    if (model::surfaces[i_surf]->bc_ == Surface::Bc::PERIODIC) {
       // Downcast to the PeriodicSurface type.
       Surface* surf_base = model::surfaces[i_surf].get();
       auto surf = dynamic_cast<PeriodicSurface*>(surf_base);
@@ -1298,7 +1289,7 @@ void read_surfaces(pugi::xml_node node)
       }
 
       // Make sure the opposite surface is also periodic.
-      if (model::surfaces[surf->i_periodic_]->bc_ != BC_PERIODIC) {
+      if (model::surfaces[surf->i_periodic_]->bc_ != Surface::Bc::PERIODIC) {
         std::stringstream err_msg;
         err_msg << "Could not find matching surface for periodic boundary "
                    "condition on surface " << surf->id_;
@@ -1311,12 +1302,12 @@ void read_surfaces(pugi::xml_node node)
   // surface
   bool boundary_exists = false;
   for (const auto& surf : model::surfaces) {
-    if (surf->bc_ != BC_TRANSMIT) {
+    if (surf->bc_ != Surface::Bc::TRANSMIT) {
       boundary_exists = true;
       break;
     }
   }
-  if (settings::run_mode != RUN_MODE_PLOTTING && !boundary_exists) {
+  if (settings::run_mode != RunMode::PLOTTING && !boundary_exists) {
     fatal_error("No boundary conditions were applied to any surfaces!");
   }
 }
diff --git a/src/tallies/derivative.cpp b/src/tallies/derivative.cpp
index f6f10cf7a..c538a3921 100644
--- a/src/tallies/derivative.cpp
+++ b/src/tallies/derivative.cpp
@@ -41,10 +41,10 @@ TallyDerivative::TallyDerivative(pugi::xml_node node)
   std::string variable_str = get_node_value(node, "variable");
 
   if (variable_str == "density") {
-    variable = DIFF_DENSITY;
+    variable = WithRespectTo::DENSITY;
 
   } else if (variable_str == "nuclide_density") {
-    variable = DIFF_NUCLIDE_DENSITY;
+    variable = WithRespectTo::NUCLIDE_DENSITY;
 
     std::string nuclide_name = get_node_value(node, "nuclide");
     bool found = false;
@@ -62,7 +62,7 @@ TallyDerivative::TallyDerivative(pugi::xml_node node)
     }
 
   } else if (variable_str == "temperature") {
-    variable = DIFF_TEMPERATURE;
+    variable = WithRespectTo::TEMPERATURE;
 
   } else {
     std::stringstream out;
@@ -146,11 +146,11 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
   // d_c / d_rho = sigma_MT * const
   // (1 / c) * (d_c / d_rho) = 1 / rho
 
-  case DIFF_DENSITY:
+  case WithRespectTo::DENSITY:
     switch (tally.estimator_) {
 
-    case ESTIMATOR_ANALOG:
-    case ESTIMATOR_COLLISION:
+    case TallyEstimator::ANALOG:
+    case TallyEstimator::COLLISION:
       switch (score_bin) {
 
       case SCORE_TOTAL:
@@ -186,10 +186,10 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
   // (1 / c) * (d_c / d_N) = sigma_MT_i / Sigma_MT
   // where i is the perturbed nuclide.
 
-  case DIFF_NUCLIDE_DENSITY:
+  case WithRespectTo::NUCLIDE_DENSITY:
     switch (tally.estimator_) {
 
-    case ESTIMATOR_ANALOG:
+    case TallyEstimator::ANALOG:
       if (p->event_nuclide_ != deriv.diff_nuclide) {
         score *= flux_deriv;
         return;
@@ -217,7 +217,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
       }
       break;
 
-    case ESTIMATOR_COLLISION:
+    case TallyEstimator::COLLISION:
       switch (score_bin) {
 
       case SCORE_TOTAL:
@@ -312,10 +312,10 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
   // computed by multipole_deriv_eval.  It only works for the resolved
   // resonance range and requires multipole data.
 
-  case DIFF_TEMPERATURE:
+  case WithRespectTo::TEMPERATURE:
     switch (tally.estimator_) {
 
-    case ESTIMATOR_ANALOG:
+    case TallyEstimator::ANALOG:
       {
         // Find the index of the event nuclide.
         int i;
@@ -401,7 +401,7 @@ apply_derivative_to_score(const Particle* p, int i_tally, int i_nuclide,
       }
       break;
 
-    case ESTIMATOR_COLLISION:
+    case TallyEstimator::COLLISION:
       if (i_nuclide != -1) {
         const auto& nuc {data::nuclides[i_nuclide]};
         if (!multipole_in_range(nuc.get(), p->E_last_)) {
@@ -578,7 +578,7 @@ score_track_derivative(const Particle* p, double distance)
 
     switch (deriv.variable) {
 
-    case DIFF_DENSITY:
+    case WithRespectTo::DENSITY:
       // phi is proportional to e^(-Sigma_tot * dist)
       // (1 / phi) * (d_phi / d_rho) = - (d_Sigma_tot / d_rho) * dist
       // (1 / phi) * (d_phi / d_rho) = - Sigma_tot / rho * dist
@@ -586,7 +586,7 @@ score_track_derivative(const Particle* p, double distance)
         / material.density_gpcc_;
       break;
 
-    case DIFF_NUCLIDE_DENSITY:
+    case WithRespectTo::NUCLIDE_DENSITY:
       // 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
@@ -594,7 +594,7 @@ score_track_derivative(const Particle* p, double distance)
         * p->neutron_xs_[deriv.diff_nuclide].total;
       break;
 
-    case DIFF_TEMPERATURE:
+    case WithRespectTo::TEMPERATURE:
       for (auto i = 0; i < material.nuclide_.size(); ++i) {
         const auto& nuc {*data::nuclides[material.nuclide_[i]]};
         if (multipole_in_range(&nuc, p->E_last_)) {
@@ -625,14 +625,14 @@ void score_collision_derivative(const Particle* p)
 
     switch (deriv.variable) {
 
-    case DIFF_DENSITY:
+    case WithRespectTo::DENSITY:
       // phi is proportional to Sigma_s
       // (1 / phi) * (d_phi / d_rho) = (d_Sigma_s / d_rho) / Sigma_s
       // (1 / phi) * (d_phi / d_rho) = 1 / rho
       deriv.flux_deriv += 1. / material.density_gpcc_;
       break;
 
-    case DIFF_NUCLIDE_DENSITY:
+    case WithRespectTo::NUCLIDE_DENSITY:
       if (p->event_nuclide_ != deriv.diff_nuclide) continue;
       // Find the index in this material for the diff_nuclide.
       int i;
@@ -653,7 +653,7 @@ void score_collision_derivative(const Particle* p)
       deriv.flux_deriv += 1. / material.atom_density_(i);
       break;
 
-    case DIFF_TEMPERATURE:
+    case WithRespectTo::TEMPERATURE:
       // Loop over the material's nuclides until we find the event nuclide.
       for (auto i_nuc : material.nuclide_) {
         const auto& nuc {*data::nuclides[i_nuc]};
diff --git a/src/tallies/filter_azimuthal.cpp b/src/tallies/filter_azimuthal.cpp
index eb02319c2..2cf03a127 100644
--- a/src/tallies/filter_azimuthal.cpp
+++ b/src/tallies/filter_azimuthal.cpp
@@ -50,11 +50,11 @@ void AzimuthalFilter::set_bins(gsl::span<double> bins)
 }
 
 void
-AzimuthalFilter::get_all_bins(const Particle* p, int estimator,
+AzimuthalFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                               FilterMatch& match) const
 {
   double phi;
-  if (estimator == ESTIMATOR_TRACKLENGTH) {
+  if (estimator == TallyEstimator::TRACKLENGTH) {
     phi = std::atan2(p->u().y, p->u().x);
   } else {
     phi = std::atan2(p->u_last_.y, p->u_last_.x);
diff --git a/src/tallies/filter_cell.cpp b/src/tallies/filter_cell.cpp
index 9e08c0c21..3077c42ea 100644
--- a/src/tallies/filter_cell.cpp
+++ b/src/tallies/filter_cell.cpp
@@ -48,7 +48,7 @@ CellFilter::set_cells(gsl::span<int32_t> cells)
 }
 
 void
-CellFilter::get_all_bins(const Particle* p, int estimator,
+CellFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                          FilterMatch& match) const
 {
   for (int i = 0; i < p->n_coord_; i++) {
diff --git a/src/tallies/filter_cell_instance.cpp b/src/tallies/filter_cell_instance.cpp
index a79e3017c..3417ca077 100644
--- a/src/tallies/filter_cell_instance.cpp
+++ b/src/tallies/filter_cell_instance.cpp
@@ -54,7 +54,7 @@ CellInstanceFilter::set_cell_instances(gsl::span<CellInstance> instances)
     Expects(x.index_cell >= 0);
     Expects(x.index_cell < model::cells.size());
     const auto& c {model::cells[x.index_cell]};
-    if (c->type_ != FILL_MATERIAL) {
+    if (c->type_ != Fill::MATERIAL) {
       throw std::invalid_argument{"Cell " + std::to_string(c->id_) + " is not "
         "filled with a material. Only material cells can be used in a cell "
         "instance filter."};
@@ -67,7 +67,7 @@ CellInstanceFilter::set_cell_instances(gsl::span<CellInstance> instances)
 }
 
 void
-CellInstanceFilter::get_all_bins(const Particle* p, int estimator,
+CellInstanceFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                          FilterMatch& match) const
 {
   gsl::index index_cell = p->coord_[p->n_coord_ - 1].cell;
diff --git a/src/tallies/filter_cellborn.cpp b/src/tallies/filter_cellborn.cpp
index 150962297..6743b7b90 100644
--- a/src/tallies/filter_cellborn.cpp
+++ b/src/tallies/filter_cellborn.cpp
@@ -5,7 +5,7 @@
 namespace openmc {
 
 void
-CellbornFilter::get_all_bins(const Particle* p, int estimator,
+CellbornFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                              FilterMatch& match) const
 {
   auto search = map_.find(p->cell_born_);
diff --git a/src/tallies/filter_cellfrom.cpp b/src/tallies/filter_cellfrom.cpp
index b8e5d0b07..49c0f4f79 100644
--- a/src/tallies/filter_cellfrom.cpp
+++ b/src/tallies/filter_cellfrom.cpp
@@ -5,7 +5,7 @@
 namespace openmc {
 
 void
-CellFromFilter::get_all_bins(const Particle* p, int estimator,
+CellFromFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                              FilterMatch& match) const
 {
   for (int i = 0; i < p->n_coord_last_; i++) {
diff --git a/src/tallies/filter_delayedgroup.cpp b/src/tallies/filter_delayedgroup.cpp
index 2b6978308..65076e2ac 100644
--- a/src/tallies/filter_delayedgroup.cpp
+++ b/src/tallies/filter_delayedgroup.cpp
@@ -37,7 +37,7 @@ DelayedGroupFilter::set_groups(gsl::span<int> groups)
 }
 
 void
-DelayedGroupFilter::get_all_bins(const Particle* p, int estimator,
+DelayedGroupFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                                  FilterMatch& match) const
 {
   match.bins_.push_back(0);
diff --git a/src/tallies/filter_distribcell.cpp b/src/tallies/filter_distribcell.cpp
index f9c7ab4d6..9efe31759 100644
--- a/src/tallies/filter_distribcell.cpp
+++ b/src/tallies/filter_distribcell.cpp
@@ -38,16 +38,16 @@ DistribcellFilter::set_cell(int32_t cell)
 }
 
 void
-DistribcellFilter::get_all_bins(const Particle* p, int estimator,
+DistribcellFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                                 FilterMatch& match) const
 {
   int offset = 0;
   auto distribcell_index = model::cells[cell_]->distribcell_index_;
   for (int i = 0; i < p->n_coord_; i++) {
     auto& c {*model::cells[p->coord_[i].cell]};
-    if (c.type_ == FILL_UNIVERSE) {
+    if (c.type_ == Fill::UNIVERSE) {
       offset += c.offset_[distribcell_index];
-    } else if (c.type_ == FILL_LATTICE) {
+    } else if (c.type_ == Fill::LATTICE) {
       auto& lat {*model::lattices[p->coord_[i+1].lattice]};
       int i_xyz[3] {p->coord_[i+1].lattice_x,
                     p->coord_[i+1].lattice_y,
diff --git a/src/tallies/filter_energy.cpp b/src/tallies/filter_energy.cpp
index 7c8b5fbee..a2c69c947 100644
--- a/src/tallies/filter_energy.cpp
+++ b/src/tallies/filter_energy.cpp
@@ -56,11 +56,11 @@ EnergyFilter::set_bins(gsl::span<const double> bins)
 }
 
 void
-EnergyFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+EnergyFilter::get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
 const
 {
   if (p->g_ != F90_NONE && matches_transport_groups_) {
-    if (estimator == ESTIMATOR_TRACKLENGTH) {
+    if (estimator == TallyEstimator::TRACKLENGTH) {
       match.bins_.push_back(data::mg.num_energy_groups_ - p->g_ - 1);
     } else {
       match.bins_.push_back(data::mg.num_energy_groups_ - p->g_last_ - 1);
@@ -100,7 +100,7 @@ EnergyFilter::text_label(int bin) const
 //==============================================================================
 
 void
-EnergyoutFilter::get_all_bins(const Particle* p, int estimator,
+EnergyoutFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                               FilterMatch& match) const
 {
   if (p->g_ != F90_NONE && matches_transport_groups_) {
diff --git a/src/tallies/filter_energyfunc.cpp b/src/tallies/filter_energyfunc.cpp
index 3e3fbd909..59cc4fcaa 100644
--- a/src/tallies/filter_energyfunc.cpp
+++ b/src/tallies/filter_energyfunc.cpp
@@ -55,7 +55,7 @@ EnergyFunctionFilter::set_data(gsl::span<const double> energy,
 }
 
 void
-EnergyFunctionFilter::get_all_bins(const Particle* p, int estimator,
+EnergyFunctionFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                                    FilterMatch& match) const
 {
   if (p->E_last_ >= energy_.front() && p->E_last_ <= energy_.back()) {
diff --git a/src/tallies/filter_legendre.cpp b/src/tallies/filter_legendre.cpp
index 3961e06d4..240054dde 100644
--- a/src/tallies/filter_legendre.cpp
+++ b/src/tallies/filter_legendre.cpp
@@ -24,7 +24,7 @@ LegendreFilter::set_order(int order)
 }
 
 void
-LegendreFilter::get_all_bins(const Particle* p, int estimator,
+LegendreFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                              FilterMatch& match) const
 {
   std::vector<double> wgt(n_bins_);
diff --git a/src/tallies/filter_material.cpp b/src/tallies/filter_material.cpp
index bfc216244..f0c2bbf8e 100644
--- a/src/tallies/filter_material.cpp
+++ b/src/tallies/filter_material.cpp
@@ -47,7 +47,7 @@ MaterialFilter::set_materials(gsl::span<const int32_t> materials)
 }
 
 void
-MaterialFilter::get_all_bins(const Particle* p, int estimator,
+MaterialFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                              FilterMatch& match) const
 {
   auto search = map_.find(p->material_);
diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp
index e6fdad31e..f80e31dc1 100644
--- a/src/tallies/filter_mesh.cpp
+++ b/src/tallies/filter_mesh.cpp
@@ -31,10 +31,10 @@ MeshFilter::from_xml(pugi::xml_node node)
 }
 
 void
-MeshFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+MeshFilter::get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
 const
 {
-  if (estimator != ESTIMATOR_TRACKLENGTH) {
+  if (estimator != TallyEstimator::TRACKLENGTH) {
     auto bin = model::meshes[mesh_]->get_bin(p->r());
     if (bin >= 0) {
       match.bins_.push_back(bin);
diff --git a/src/tallies/filter_meshsurface.cpp b/src/tallies/filter_meshsurface.cpp
index edb0492e9..80e398880 100644
--- a/src/tallies/filter_meshsurface.cpp
+++ b/src/tallies/filter_meshsurface.cpp
@@ -8,7 +8,7 @@
 namespace openmc {
 
 void
-MeshSurfaceFilter::get_all_bins(const Particle* p, int estimator,
+MeshSurfaceFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                                 FilterMatch& match) const
 {
   model::meshes[mesh_]->surface_bins_crossed(p, match.bins_);
@@ -23,47 +23,47 @@ MeshSurfaceFilter::text_label(int bin) const
 
   // Get flattend mesh index and surface index.
   int i_mesh = bin / (4 * n_dim);
-  int i_surf = (bin % (4 * n_dim)) + 1;
+  MeshDir surf_dir = static_cast<MeshDir>(bin % (4 * n_dim));
 
   // Get mesh index part of label.
   std::string out = MeshFilter::text_label(i_mesh);
 
   // Get surface part of label.
-  switch (i_surf) {
-    case OUT_LEFT:
+  switch (surf_dir) {
+    case MeshDir::OUT_LEFT:
       out += " Outgoing, x-min";
       break;
-    case IN_LEFT:
+    case MeshDir::IN_LEFT:
       out += " Incoming, x-min";
       break;
-    case OUT_RIGHT:
+    case MeshDir::OUT_RIGHT:
       out += " Outgoing, x-max";
       break;
-    case IN_RIGHT:
+    case MeshDir::IN_RIGHT:
       out += " Incoming, x-max";
       break;
-    case OUT_BACK:
+    case MeshDir::OUT_BACK:
       out += " Outgoing, y-min";
       break;
-    case IN_BACK:
+    case MeshDir::IN_BACK:
       out += " Incoming, y-min";
       break;
-    case OUT_FRONT:
+    case MeshDir::OUT_FRONT:
       out += " Outgoing, y-max";
       break;
-    case IN_FRONT:
+    case MeshDir::IN_FRONT:
       out += " Incoming, y-max";
       break;
-    case OUT_BOTTOM:
+    case MeshDir::OUT_BOTTOM:
       out += " Outgoing, z-min";
       break;
-    case IN_BOTTOM:
+    case MeshDir::IN_BOTTOM:
       out += " Incoming, z-min";
       break;
-    case OUT_TOP:
+    case MeshDir::OUT_TOP:
       out += " Outgoing, z-max";
       break;
-    case IN_TOP:
+    case MeshDir::IN_TOP:
       out += " Incoming, z-max";
       break;
   }
diff --git a/src/tallies/filter_mu.cpp b/src/tallies/filter_mu.cpp
index 9aaebea8c..a345fe2bd 100644
--- a/src/tallies/filter_mu.cpp
+++ b/src/tallies/filter_mu.cpp
@@ -49,7 +49,7 @@ MuFilter::set_bins(gsl::span<double> bins)
 }
 
 void
-MuFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+MuFilter::get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
 const
 {
   if (p->mu_ >= bins_.front() && p->mu_ <= bins_.back()) {
diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp
index eb419fec9..caccbc961 100644
--- a/src/tallies/filter_particle.cpp
+++ b/src/tallies/filter_particle.cpp
@@ -40,7 +40,7 @@ ParticleFilter::set_particles(gsl::span<Particle::Type> particles)
 }
 
 void
-ParticleFilter::get_all_bins(const Particle* p, int estimator,
+ParticleFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                              FilterMatch& match) const
 {
   for (auto i = 0; i < particles_.size(); i++) {
diff --git a/src/tallies/filter_polar.cpp b/src/tallies/filter_polar.cpp
index edbbc8922..f35c1d3cd 100644
--- a/src/tallies/filter_polar.cpp
+++ b/src/tallies/filter_polar.cpp
@@ -50,11 +50,11 @@ PolarFilter::set_bins(gsl::span<double> bins)
 }
 
 void
-PolarFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match)
+PolarFilter::get_all_bins(const Particle* p, TallyEstimator estimator, FilterMatch& match)
 const
 {
   double theta;
-  if (estimator == ESTIMATOR_TRACKLENGTH) {
+  if (estimator == TallyEstimator::TRACKLENGTH) {
     theta = std::acos(p->u().z);
   } else {
     theta = std::acos(p->u_last_.z);
diff --git a/src/tallies/filter_sph_harm.cpp b/src/tallies/filter_sph_harm.cpp
index 80d6878b7..bbeebf5b1 100644
--- a/src/tallies/filter_sph_harm.cpp
+++ b/src/tallies/filter_sph_harm.cpp
@@ -44,7 +44,7 @@ SphericalHarmonicsFilter::set_cosine(gsl::cstring_span cosine)
 }
 
 void
-SphericalHarmonicsFilter::get_all_bins(const Particle* p, int estimator,
+SphericalHarmonicsFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                                        FilterMatch& match) const
 {
   // Determine cosine term for scatter expansion if necessary
diff --git a/src/tallies/filter_sptl_legendre.cpp b/src/tallies/filter_sptl_legendre.cpp
index 40d82f4c1..71cc3f399 100644
--- a/src/tallies/filter_sptl_legendre.cpp
+++ b/src/tallies/filter_sptl_legendre.cpp
@@ -61,7 +61,7 @@ SpatialLegendreFilter::set_minmax(double min, double max)
 }
 
 void
-SpatialLegendreFilter::get_all_bins(const Particle* p, int estimator,
+SpatialLegendreFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                                     FilterMatch& match) const
 {
   // Get the coordinate along the axis of interest.
diff --git a/src/tallies/filter_surface.cpp b/src/tallies/filter_surface.cpp
index 72b356f21..5613f039a 100644
--- a/src/tallies/filter_surface.cpp
+++ b/src/tallies/filter_surface.cpp
@@ -49,7 +49,7 @@ SurfaceFilter::set_surfaces(gsl::span<int32_t> surfaces)
 }
 
 void
-SurfaceFilter::get_all_bins(const Particle* p, int estimator,
+SurfaceFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                             FilterMatch& match) const
 {
   auto search = map_.find(std::abs(p->surface_)-1);
diff --git a/src/tallies/filter_universe.cpp b/src/tallies/filter_universe.cpp
index dffdee621..289e96f56 100644
--- a/src/tallies/filter_universe.cpp
+++ b/src/tallies/filter_universe.cpp
@@ -47,7 +47,7 @@ UniverseFilter::set_universes(gsl::span<int32_t> universes)
 }
 
 void
-UniverseFilter::get_all_bins(const Particle* p, int estimator,
+UniverseFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                              FilterMatch& match) const
 {
   for (int i = 0; i < p->n_coord_; i++) {
diff --git a/src/tallies/filter_zernike.cpp b/src/tallies/filter_zernike.cpp
index b732ef793..826d9f508 100644
--- a/src/tallies/filter_zernike.cpp
+++ b/src/tallies/filter_zernike.cpp
@@ -25,7 +25,7 @@ ZernikeFilter::from_xml(pugi::xml_node node)
 }
 
 void
-ZernikeFilter::get_all_bins(const Particle* p, int estimator,
+ZernikeFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                             FilterMatch& match) const
 {
   // Determine the normalized (r,theta) coordinates.
@@ -86,7 +86,7 @@ ZernikeFilter::set_order(int order)
 //==============================================================================
 
 void
-ZernikeRadialFilter::get_all_bins(const Particle* p, int estimator,
+ZernikeRadialFilter::get_all_bins(const Particle* p, TallyEstimator estimator,
                                   FilterMatch& match) const
 {
   // Determine the normalized radius coordinate.
diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp
index d65f4d7ac..a9e757ada 100644
--- a/src/tallies/tally.cpp
+++ b/src/tallies/tally.cpp
@@ -310,15 +310,15 @@ Tally::Tally(pugi::xml_node node)
     // Change the tally estimator if a filter demands it
     std::string filt_type = f->type();
     if (filt_type == "energyout" || filt_type == "legendre") {
-      estimator_ = ESTIMATOR_ANALOG;
+      estimator_ = TallyEstimator::ANALOG;
     } else if (filt_type == "sphericalharmonics") {
       auto sf = dynamic_cast<SphericalHarmonicsFilter*>(f);
       if (sf->cosine() == SphericalHarmonicsCosine::scatter) {
-        estimator_ = ESTIMATOR_ANALOG;
+        estimator_ = TallyEstimator::ANALOG;
       }
     } else if (filt_type == "spatiallegendre" || filt_type == "zernike"
       || filt_type == "zernikeradial") {
-      estimator_ = ESTIMATOR_COLLISION;
+      estimator_ = TallyEstimator::COLLISION;
     }
   }
 
@@ -369,7 +369,7 @@ Tally::Tally(pugi::xml_node node)
       for (auto p : pf->particles()) {
         if (p == Particle::Type::electron ||
             p == Particle::Type::positron) {
-          estimator_ = ESTIMATOR_ANALOG;
+          estimator_ = TallyEstimator::ANALOG;
         }
       }
     }
@@ -402,13 +402,13 @@ Tally::Tally(pugi::xml_node node)
 
     // Only analog or collision estimators are supported for differential
     // tallies.
-    if (estimator_ == ESTIMATOR_TRACKLENGTH) {
-      estimator_ = ESTIMATOR_COLLISION;
+    if (estimator_ == TallyEstimator::TRACKLENGTH) {
+      estimator_ = TallyEstimator::COLLISION;
     }
 
     const auto& deriv = model::tally_derivs[deriv_];
-    if (deriv.variable == DIFF_NUCLIDE_DENSITY
-      || deriv.variable == DIFF_TEMPERATURE) {
+    if (deriv.variable == WithRespectTo::NUCLIDE_DENSITY
+      || deriv.variable == WithRespectTo::TEMPERATURE) {
       for (int i_nuc : nuclides_) {
         if (has_energyout && i_nuc == -1) {
           fatal_error("Error on tally " + std::to_string(id_)
@@ -437,29 +437,29 @@ Tally::Tally(pugi::xml_node node)
   if (check_for_node(node, "estimator")) {
     std::string est = get_node_value(node, "estimator");
     if (est == "analog") {
-      estimator_ = ESTIMATOR_ANALOG;
+      estimator_ = TallyEstimator::ANALOG;
     } else if (est == "tracklength" || est == "track-length"
       || est == "pathlength" || est == "path-length") {
       // If the estimator was set to an analog estimator, this means the
       // tally needs post-collision information
-      if (estimator_ == ESTIMATOR_ANALOG) {
+      if (estimator_ == TallyEstimator::ANALOG) {
         throw std::runtime_error{"Cannot use track-length estimator for tally "
           + std::to_string(id_)};
       }
 
       // Set estimator to track-length estimator
-      estimator_ = ESTIMATOR_TRACKLENGTH;
+      estimator_ = TallyEstimator::TRACKLENGTH;
 
     } else if (est == "collision") {
       // If the estimator was set to an analog estimator, this means the
       // tally needs post-collision information
-      if (estimator_ == ESTIMATOR_ANALOG) {
+      if (estimator_ == TallyEstimator::ANALOG) {
         throw std::runtime_error{"Cannot use collision estimator for tally " +
           std::to_string(id_)};
       }
 
       // Set estimator to collision estimator
-      estimator_ = ESTIMATOR_COLLISION;
+      estimator_ = TallyEstimator::COLLISION;
 
     } else {
       throw std::runtime_error{"Invalid estimator '" + est + "' on tally " +
@@ -615,20 +615,20 @@ Tally::set_scores(const std::vector<std::string>& scores)
 
     case SCORE_SCATTER:
       if (legendre_present)
-        estimator_ = ESTIMATOR_ANALOG;
+        estimator_ = TallyEstimator::ANALOG;
     case SCORE_NU_FISSION:
     case SCORE_DELAYED_NU_FISSION:
     case SCORE_PROMPT_NU_FISSION:
       if (energyout_present)
-        estimator_ = ESTIMATOR_ANALOG;
+        estimator_ = TallyEstimator::ANALOG;
       break;
 
     case SCORE_NU_SCATTER:
       if (settings::run_CE) {
-        estimator_ = ESTIMATOR_ANALOG;
+        estimator_ = TallyEstimator::ANALOG;
       } else {
         if (energyout_present || legendre_present)
-          estimator_ = ESTIMATOR_ANALOG;
+          estimator_ = TallyEstimator::ANALOG;
       }
       break;
 
@@ -647,9 +647,9 @@ Tally::set_scores(const std::vector<std::string>& scores)
         if (meshsurface_present)
           fatal_error("Cannot tally mesh surface currents in the same tally as "
             "normal surface currents");
-        type_ = TALLY_SURFACE;
+        type_ = TallyType::SURFACE;
       } else if (meshsurface_present) {
-        type_ = TALLY_MESH_SURFACE;
+        type_ = TallyType::MESH_SURFACE;
       } else {
         fatal_error("Cannot tally currents without surface type filters");
       }
@@ -677,7 +677,7 @@ Tally::set_scores(const std::vector<std::string>& scores)
   }
 
   // Make sure current scores are not mixed in with volumetric scores.
-  if (type_ == TALLY_SURFACE || type_ == TALLY_MESH_SURFACE) {
+  if (type_ == TallyType::SURFACE || type_ == TallyType::MESH_SURFACE) {
     if (scores_.size() != 1)
       fatal_error("Cannot tally other scores in the same tally as surface "
         "currents");
@@ -823,7 +823,7 @@ void Tally::accumulate()
   if (mpi::master || !settings::reduce_tallies) {
     // Calculate total source strength for normalization
     double total_source = 0.0;
-    if (settings::run_mode == RUN_MODE_FIXEDSOURCE) {
+    if (settings::run_mode == RunMode::FIXEDSOURCE) {
       for (const auto& s : model::external_sources) {
         total_source += s.strength();
       }
@@ -837,10 +837,10 @@ void Tally::accumulate()
     // Accumulate each result
     for (int i = 0; i < results_.shape()[0]; ++i) {
       for (int j = 0; j < results_.shape()[1]; ++j) {
-        double val = results_(i, j, RESULT_VALUE) * norm;
-        results_(i, j, RESULT_VALUE) = 0.0;
-        results_(i, j, RESULT_SUM) += val;
-        results_(i, j, RESULT_SUM_SQ) += val*val;
+        double val = results_(i, j, TallyResult::VALUE) * norm;
+        results_(i, j, TallyResult::VALUE) = 0.0;
+        results_(i, j, TallyResult::SUM) += val;
+        results_(i, j, TallyResult::SUM_SQ) += val*val;
       }
     }
   }
@@ -872,7 +872,7 @@ void read_tallies_xml()
   read_meshes(root);
 
   // We only need the mesh info for plotting
-  if (settings::run_mode == RUN_MODE_PLOTTING) return;
+  if (settings::run_mode == RunMode::PLOTTING) return;
 
   // Read data for tally derivatives
   read_tally_derivatives(root);
@@ -908,7 +908,7 @@ void reduce_tally_results()
     auto& tally {model::tallies[i_tally]};
 
     // Get view of accumulated tally values
-    auto values_view = xt::view(tally->results_, xt::all(), xt::all(), RESULT_VALUE);
+    auto values_view = xt::view(tally->results_, xt::all(), xt::all(), TallyResult::VALUE);
 
     // Make copy of tally values in contiguous array
     xt::xtensor<double, 2> values = values_view;
@@ -928,7 +928,7 @@ void reduce_tally_results()
 
   // Get view of global tally values
   auto& gt = simulation::global_tallies;
-  auto gt_values_view = xt::view(gt, xt::all(), RESULT_VALUE);
+  auto gt_values_view = xt::view(gt, xt::all(), TallyResult::VALUE);
 
   // Make copy of values in contiguous array
   xt::xtensor<double, 1> gt_values = gt_values_view;
@@ -969,12 +969,12 @@ accumulate_tallies()
   if (mpi::master || !settings::reduce_tallies) {
     auto& gt = simulation::global_tallies;
 
-    if (settings::run_mode == RUN_MODE_EIGENVALUE) {
+    if (settings::run_mode == RunMode::EIGENVALUE) {
       if (simulation::current_batch > settings::n_inactive) {
         // Accumulate products of different estimators of k
-        double k_col = gt(K_COLLISION, RESULT_VALUE) / simulation::total_weight;
-        double k_abs = gt(K_ABSORPTION, RESULT_VALUE) / simulation::total_weight;
-        double k_tra = gt(K_TRACKLENGTH, RESULT_VALUE) / simulation::total_weight;
+        double k_col = gt(GlobalTally::K_COLLISION, TallyResult::VALUE) / simulation::total_weight;
+        double k_abs = gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) / simulation::total_weight;
+        double k_tra = gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) / simulation::total_weight;
         simulation::k_col_abs += k_col * k_abs;
         simulation::k_col_tra += k_col * k_tra;
         simulation::k_abs_tra += k_abs * k_tra;
@@ -983,10 +983,10 @@ accumulate_tallies()
 
     // Accumulate results for global tallies
     for (int i = 0; i < N_GLOBAL_TALLIES; ++i) {
-      double val = gt(i, RESULT_VALUE)/simulation::total_weight;
-      gt(i, RESULT_VALUE) = 0.0;
-      gt(i, RESULT_SUM) += val;
-      gt(i, RESULT_SUM_SQ) += val*val;
+      double val = gt(i, TallyResult::VALUE)/simulation::total_weight;
+      gt(i, TallyResult::VALUE) = 0.0;
+      gt(i, TallyResult::SUM) += val;
+      gt(i, TallyResult::SUM_SQ) += val*val;
     }
   }
 
@@ -1014,24 +1014,24 @@ setup_active_tallies()
       model::active_tallies.push_back(i);
       switch (tally.type_) {
 
-      case TALLY_VOLUME:
+      case TallyType::VOLUME:
         switch (tally.estimator_) {
-          case ESTIMATOR_ANALOG:
+          case TallyEstimator::ANALOG:
             model::active_analog_tallies.push_back(i);
             break;
-          case ESTIMATOR_TRACKLENGTH:
+          case TallyEstimator::TRACKLENGTH:
             model::active_tracklength_tallies.push_back(i);
             break;
-          case ESTIMATOR_COLLISION:
+          case TallyEstimator::COLLISION:
             model::active_collision_tallies.push_back(i);
         }
         break;
 
-      case TALLY_MESH_SURFACE:
+      case TallyType::MESH_SURFACE:
         model::active_meshsurf_tallies.push_back(i);
         break;
 
-      case TALLY_SURFACE:
+      case TallyType::SURFACE:
         model::active_surface_tallies.push_back(i);
       }
 
@@ -1111,7 +1111,7 @@ openmc_tally_get_estimator(int32_t index, int* estimator)
     return OPENMC_E_OUT_OF_BOUNDS;
   }
 
-  *estimator = model::tallies[index]->estimator_;
+  *estimator = static_cast<int>(model::tallies[index]->estimator_);
   return 0;
 }
 
@@ -1127,11 +1127,11 @@ openmc_tally_set_estimator(int32_t index, const char* estimator)
 
   std::string est = estimator;
   if (est == "analog") {
-    t->estimator_ = ESTIMATOR_ANALOG;
+    t->estimator_ = TallyEstimator::ANALOG;
   } else if (est == "collision") {
-    t->estimator_ = ESTIMATOR_COLLISION;
+    t->estimator_ = TallyEstimator::COLLISION;
   } else if (est == "tracklength") {
-    t->estimator_ = ESTIMATOR_TRACKLENGTH;
+    t->estimator_ = TallyEstimator::TRACKLENGTH;
   } else {
     set_errmsg("Unknown tally estimator: " + est);
     return OPENMC_E_INVALID_ARGUMENT;
@@ -1170,7 +1170,7 @@ openmc_tally_get_type(int32_t index, int32_t* type)
     set_errmsg("Index in tallies array is out of bounds.");
     return OPENMC_E_OUT_OF_BOUNDS;
   }
-  *type = model::tallies[index]->type_;
+  *type = static_cast<int>(model::tallies[index]->type_);
 
   return 0;
 }
@@ -1183,11 +1183,11 @@ openmc_tally_set_type(int32_t index, const char* type)
     return OPENMC_E_OUT_OF_BOUNDS;
   }
   if (strcmp(type, "volume") == 0) {
-    model::tallies[index]->type_ = TALLY_VOLUME;
+    model::tallies[index]->type_ = TallyType::VOLUME;
   } else if (strcmp(type, "mesh-surface") == 0) {
-    model::tallies[index]->type_ = TALLY_MESH_SURFACE;
+    model::tallies[index]->type_ = TallyType::MESH_SURFACE;
   } else if (strcmp(type, "surface") == 0) {
-    model::tallies[index]->type_ = TALLY_SURFACE;
+    model::tallies[index]->type_ = TallyType::SURFACE;
   } else {
     std::stringstream errmsg;
     errmsg << "Unknown tally type: " << type;
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index 51eae2620..df6433a64 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -165,7 +165,7 @@ score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index)
 
   // Update the tally result
   #pragma omp atomic
-  tally.results_(filter_index, score_index, RESULT_VALUE) += score;
+  tally.results_(filter_index, score_index, TallyResult::VALUE) += score;
 
   // Reset the original delayed group bin
   dg_match.bins_[i_bin] = original_bin;
@@ -195,7 +195,7 @@ double get_nuc_fission_q(const Nuclide& nuc, const Particle* p, int score_bin)
 double score_fission_q(const Particle* p, int score_bin, const Tally& tally,
   double flux, int i_nuclide, double atom_density)
 {
-  if (tally.estimator_ == ESTIMATOR_ANALOG) {
+  if (tally.estimator_ == TallyEstimator::ANALOG) {
     const Nuclide& nuc {*data::nuclides[p->event_nuclide_]};
     if (settings::survival_biasing) {
       // No fission events occur if survival biasing is on -- need to
@@ -208,7 +208,7 @@ double score_fission_q(const Particle* p, int score_bin, const Tally& tally,
       }
     } else {
       // Skip any non-absorption events
-      if (p->event_ == EVENT_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
@@ -273,7 +273,7 @@ double score_neutron_heating(const Particle* p, const Tally& tally, double flux,
   if (i_nuclide >= 0) {
     const Nuclide& nuc {*data::nuclides[i_nuclide]};
     heating_xs = get_nuclide_neutron_heating(p, nuc, rxn_bin, i_nuclide);
-    if (tally.estimator_ == ESTIMATOR_ANALOG) {
+    if (tally.estimator_ == TallyEstimator::ANALOG) {
       heating_xs /= p->neutron_xs_[i_nuclide].total;
     } else {
       heating_xs *= atom_density;
@@ -288,13 +288,13 @@ double score_neutron_heating(const Particle* p, const Tally& tally, double flux,
         const Nuclide& nuc {*data::nuclides[j_nuclide]};
         heating_xs += atom_density * get_nuclide_neutron_heating(p, nuc, rxn_bin, j_nuclide);
       }
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         heating_xs /= p->macro_xs_.total;
       }
     }
   }
   score = heating_xs * flux;
-  if (tally.estimator_ == ESTIMATOR_ANALOG) {
+  if (tally.estimator_ == TallyEstimator::ANALOG) {
     // All events score to a heating tally bin. We actually use a
     // collision estimator in place of an analog one since there is no
     // reaction-wise heating cross section
@@ -394,7 +394,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
 
       // Update tally results
       #pragma omp atomic
-      tally.results_(filter_index, i_score, RESULT_VALUE) += score;
+      tally.results_(filter_index, i_score, TallyResult::VALUE) += score;
 
     } else if (score_bin == SCORE_DELAYED_NU_FISSION && g != 0) {
 
@@ -441,7 +441,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin)
 
         // Update tally results
         #pragma omp atomic
-        tally.results_(filter_index, i_score, RESULT_VALUE) += score*filter_weight;
+        tally.results_(filter_index, i_score, TallyResult::VALUE) += score*filter_weight;
       }
     }
   }
@@ -475,7 +475,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
 
 
     case SCORE_FLUX:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         // All events score to a flux bin. We actually use a collision estimator
         // in place of an analog one since there is no way to count 'events'
         // exactly for the flux
@@ -500,7 +500,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
 
 
     case SCORE_TOTAL:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         // All events will score to the total reaction rate. We can just use
         // use the weight of the particle entering the collision as the score
         if (settings::survival_biasing) {
@@ -522,7 +522,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
 
 
     case SCORE_INVERSE_VELOCITY:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         // All events score to an inverse velocity bin. We actually use a
         // collision estimator in place of an analog one since there is no way
         // to count 'events' exactly for the inverse velocity
@@ -543,9 +543,9 @@ score_general_ce(Particle* p, int i_tally, int start_index,
 
 
     case SCORE_SCATTER:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         // Skip any event where the particle didn't scatter
-        if (p->event_ != EVENT_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;
@@ -564,7 +564,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
     case SCORE_NU_SCATTER:
       // Only analog estimators are available.
       // Skip any event where the particle didn't scatter
-      if (p->event_ != EVENT_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
@@ -584,14 +584,14 @@ score_general_ce(Particle* p, int i_tally, int start_index,
 
 
     case SCORE_ABSORPTION:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         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;
         } else {
           // Skip any event where the particle wasn't absorbed
-          if (p->event_ == EVENT_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;
@@ -609,7 +609,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
 
     case SCORE_FISSION:
       if (p->macro_xs_.absorption == 0) continue;
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      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
@@ -623,7 +623,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
           }
         } else {
           // Skip any non-absorption events
-          if (p->event_ == EVENT_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
@@ -643,7 +643,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
 
     case SCORE_NU_FISSION:
       if (p->macro_xs_.absorption == 0) continue;
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         if (settings::survival_biasing || p->fission_) {
           if (tally.energyout_filter_ != C_NONE) {
             // Fission has multiple outgoing neutrons so this helper function
@@ -686,7 +686,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
 
     case SCORE_PROMPT_NU_FISSION:
       if (p->macro_xs_.absorption == 0) continue;
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         if (settings::survival_biasing || p->fission_) {
           if (tally.energyout_filter_ != C_NONE) {
             // Fission has multiple outgoing neutrons so this helper function
@@ -748,7 +748,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
 
     case SCORE_DELAYED_NU_FISSION:
       if (p->macro_xs_.absorption == 0) continue;
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         if (settings::survival_biasing || p->fission_) {
           if (tally.energyout_filter_ != C_NONE) {
             // Fission has multiple outgoing neutrons so this helper function
@@ -892,7 +892,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
 
     case SCORE_DECAY_RATE:
       if (p->macro_xs_.absorption == 0) continue;
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      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
@@ -1079,7 +1079,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
       score = 0.;
       // Kappa-fission values are determined from the Q-value listed for the
       // fission cross section.
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      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
@@ -1094,7 +1094,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
           }
         } else {
           // Skip any non-absorption events
-          if (p->event_ == EVENT_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
@@ -1141,7 +1141,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
 
 
     case ELASTIC:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         // Check if event MT matches
         if (p->event_mt_ != ELASTIC) continue;
         score = p->wgt_last_ * flux;
@@ -1180,7 +1180,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
     case N_GAMMA:
     case N_P:
     case N_A:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         // Check if the event MT matches
         if (p->event_mt_ != score_bin) continue;
         score = p->wgt_last_ * flux;
@@ -1219,7 +1219,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
         score = score_neutron_heating(p, tally, flux, HEATING,
             i_nuclide, atom_density);
       } else if (p->type_ == Particle::Type::photon) {
-        if (tally.estimator_ == ESTIMATOR_ANALOG) {
+        if (tally.estimator_ == TallyEstimator::ANALOG) {
           // Score direct energy deposition in the collision
           score = E - p->E_;
           // We need to substract the energy of the secondary particles since
@@ -1267,7 +1267,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
       break;
 
     default:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      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;
@@ -1326,7 +1326,7 @@ score_general_ce(Particle* p, int i_tally, int start_index,
 
     // Update tally results
     #pragma omp atomic
-    tally.results_(filter_index, score_index, RESULT_VALUE) += score;
+    tally.results_(filter_index, score_index, TallyResult::VALUE) += score;
   }
 }
 
@@ -1344,8 +1344,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
   // Set the direction and group to use with get_xs
   Direction p_u;
   int p_g;
-  if (tally.estimator_ == ESTIMATOR_ANALOG
-    || tally.estimator_ == ESTIMATOR_COLLISION) {
+  if (tally.estimator_ == TallyEstimator::ANALOG
+    || tally.estimator_ == TallyEstimator::COLLISION) {
 
     if (settings::survival_biasing) {
 
@@ -1358,7 +1358,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
         p_u = p->u_local();
         p_g = p->g_;
       }
-    } else if (p->event_ == EVENT_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
@@ -1398,7 +1398,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
 
 
     case SCORE_FLUX:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         // All events score to a flux bin. We actually use a collision estimator
         // in place of an analog one since there is no way to count 'events'
         // exactly for the flux
@@ -1417,7 +1417,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
 
 
     case SCORE_TOTAL:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         // All events will score to the total reaction rate. We can just use
         // use the weight of the particle entering the collision as the score
         if (settings::survival_biasing) {
@@ -1429,12 +1429,12 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
         }
         //TODO: should flux be multiplied in above instead of below?
         if (i_nuclide >= 0) {
-          score *= flux * atom_density * nuc_xs.get_xs(MG_GET_XS_TOTAL, p_g)
-            / macro_xs.get_xs(MG_GET_XS_TOTAL, p_g);
+          score *= flux * atom_density * nuc_xs.get_xs(MgxsType::TOTAL, p_g)
+            / macro_xs.get_xs(MgxsType::TOTAL, p_g);
         }
       } else {
         if (i_nuclide >= 0) {
-          score = atom_density * flux * nuc_xs.get_xs(MG_GET_XS_TOTAL, p_g);
+          score = atom_density * flux * nuc_xs.get_xs(MgxsType::TOTAL, p_g);
         } else {
           score = p->macro_xs_.total * flux;
         }
@@ -1443,8 +1443,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
 
 
     case SCORE_INVERSE_VELOCITY:
-      if (tally.estimator_ == ESTIMATOR_ANALOG
-        || tally.estimator_ == ESTIMATOR_COLLISION) {
+      if (tally.estimator_ == TallyEstimator::ANALOG
+        || tally.estimator_ == TallyEstimator::COLLISION) {
         // All events score to an inverse velocity bin. We actually use a
         // collision estimator in place of an analog one since there is no way
         // to count 'events' exactly for the inverse velocity
@@ -1456,51 +1456,51 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
           score = p->wgt_last_;
         }
         if (i_nuclide >= 0) {
-          score *= flux * nuc_xs.get_xs(MG_GET_XS_INVERSE_VELOCITY, p_g)
-            / macro_xs.get_xs(MG_GET_XS_TOTAL, p_g);
+          score *= flux * nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g)
+            / macro_xs.get_xs(MgxsType::TOTAL, p_g);
         } else {
-          score *= flux * macro_xs.get_xs(MG_GET_XS_INVERSE_VELOCITY, p_g)
-            / macro_xs.get_xs(MG_GET_XS_TOTAL, p_g);
+          score *= flux * macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g)
+            / macro_xs.get_xs(MgxsType::TOTAL, p_g);
         }
       } else {
         if (i_nuclide >= 0) {
-          score = flux * nuc_xs.get_xs(MG_GET_XS_INVERSE_VELOCITY, p_g);
+          score = flux * nuc_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g);
         } else {
-          score = flux * macro_xs.get_xs(MG_GET_XS_INVERSE_VELOCITY, p_g);
+          score = flux * macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, p_g);
         }
       }
       break;
 
 
     case SCORE_SCATTER:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         // Skip any event where the particle didn't scatter
-        if (p->event_ != EVENT_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;
         if (i_nuclide >= 0) {
           score *= atom_density * nuc_xs.get_xs(
-              MG_GET_XS_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(
-              MG_GET_XS_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(
-            MG_GET_XS_SCATTER_MULT, p_g, nullptr, &p->mu_, nullptr);
+            MgxsType::SCATTER_MULT, p_g, nullptr, &p->mu_, nullptr);
         } else {
           score = flux * macro_xs.get_xs(
-            MG_GET_XS_SCATTER_MULT, p_g, nullptr, &p->mu_, nullptr);
+            MgxsType::SCATTER_MULT, p_g, nullptr, &p->mu_, nullptr);
         }
       }
       break;
 
 
     case SCORE_NU_SCATTER:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         // Skip any event where the particle didn't scatter
-        if (p->event_ != EVENT_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
@@ -1510,42 +1510,42 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
         // adjust the score by the actual probability for that nuclide.
         if (i_nuclide >= 0) {
           score *= atom_density
-            * nuc_xs.get_xs(MG_GET_XS_SCATTER_FMU, p->g_last_, &p->g_,
+            * nuc_xs.get_xs(MgxsType::SCATTER_FMU, p->g_last_, &p->g_,
                             &p->mu_, nullptr)
-            / macro_xs.get_xs(MG_GET_XS_SCATTER_FMU, p->g_last_, &p->g_,
+            / macro_xs.get_xs(MgxsType::SCATTER_FMU, p->g_last_, &p->g_,
                               &p->mu_, nullptr);
         }
       } else {
         if (i_nuclide >= 0) {
-          score = atom_density * flux * nuc_xs.get_xs(MG_GET_XS_SCATTER, p_g);
+          score = atom_density * flux * nuc_xs.get_xs(MgxsType::SCATTER, p_g);
         } else {
-          score = flux * macro_xs.get_xs(MG_GET_XS_SCATTER, p_g);
+          score = flux * macro_xs.get_xs(MgxsType::SCATTER, p_g);
         }
       }
       break;
 
 
     case SCORE_ABSORPTION:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         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;
         } else {
           // Skip any event where the particle wasn't absorbed
-          if (p->event_ == EVENT_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;
         }
         if (i_nuclide >= 0) {
-          score *= atom_density * nuc_xs.get_xs(MG_GET_XS_ABSORPTION, p_g)
-            / macro_xs.get_xs(MG_GET_XS_ABSORPTION, p_g);
+          score *= atom_density * nuc_xs.get_xs(MgxsType::ABSORPTION, p_g)
+            / macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
         }
       } else {
         if (i_nuclide >= 0) {
           score = atom_density * flux
-            * nuc_xs.get_xs(MG_GET_XS_ABSORPTION, p_g);
+            * nuc_xs.get_xs(MgxsType::ABSORPTION, p_g);
         } else {
           score = p->macro_xs_.absorption * flux;
         }
@@ -1554,7 +1554,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
 
 
     case SCORE_FISSION:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      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
@@ -1562,31 +1562,31 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
           score = p->wgt_absorb_ * flux;
         } else {
           // Skip any non-absorption events
-          if (p->event_ == EVENT_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;
         }
         if (i_nuclide >= 0) {
-          score *= atom_density * nuc_xs.get_xs(MG_GET_XS_FISSION, p_g)
-            / macro_xs.get_xs(MG_GET_XS_ABSORPTION, p_g);
+          score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g)
+            / macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
         } else {
-          score *= macro_xs.get_xs(MG_GET_XS_FISSION, p_g)
-            / macro_xs.get_xs(MG_GET_XS_ABSORPTION, p_g);
+          score *= macro_xs.get_xs(MgxsType::FISSION, p_g)
+            / macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
         }
       } else {
         if (i_nuclide >= 0) {
-          score = atom_density * flux * nuc_xs.get_xs(MG_GET_XS_FISSION, p_g);
+          score = atom_density * flux * nuc_xs.get_xs(MgxsType::FISSION, p_g);
         } else {
-          score = flux * macro_xs.get_xs(MG_GET_XS_FISSION, p_g);
+          score = flux * macro_xs.get_xs(MgxsType::FISSION, p_g);
         }
       }
       break;
 
 
     case SCORE_NU_FISSION:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         if (settings::survival_biasing || p->fission_) {
           if (tally.energyout_filter_ != C_NONE) {
             // Fission has multiple outgoing neutrons so this helper function
@@ -1601,11 +1601,11 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
           // nu-fission
           score = p->wgt_absorb_ * flux;
           if (i_nuclide >= 0) {
-            score *= atom_density * nuc_xs.get_xs(MG_GET_XS_NU_FISSION, p_g)
-              / macro_xs.get_xs(MG_GET_XS_ABSORPTION, p_g);
+            score *= atom_density * nuc_xs.get_xs(MgxsType::NU_FISSION, p_g)
+              / macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
           } else {
-            score *= macro_xs.get_xs(MG_GET_XS_NU_FISSION, p_g)
-              / macro_xs.get_xs(MG_GET_XS_ABSORPTION, p_g);
+            score *= macro_xs.get_xs(MgxsType::NU_FISSION, p_g)
+              / macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
           }
         } else {
           // Skip any non-fission events
@@ -1617,23 +1617,23 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
           // score.
           score = simulation::keff * p->wgt_bank_ * flux;
           if (i_nuclide >= 0) {
-            score *= atom_density * nuc_xs.get_xs(MG_GET_XS_FISSION, p_g)
-              / macro_xs.get_xs(MG_GET_XS_FISSION, p_g);
+            score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g)
+              / macro_xs.get_xs(MgxsType::FISSION, p_g);
           }
         }
       } else {
         if (i_nuclide >= 0) {
           score = atom_density * flux
-            * nuc_xs.get_xs(MG_GET_XS_NU_FISSION, p_g);
+            * nuc_xs.get_xs(MgxsType::NU_FISSION, p_g);
         } else {
-          score = flux * macro_xs.get_xs(MG_GET_XS_NU_FISSION, p_g);
+          score = flux * macro_xs.get_xs(MgxsType::NU_FISSION, p_g);
         }
       }
       break;
 
 
     case SCORE_PROMPT_NU_FISSION:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         if (settings::survival_biasing || p->fission_) {
           if (tally.energyout_filter_ != C_NONE) {
             // Fission has multiple outgoing neutrons so this helper function
@@ -1649,11 +1649,11 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
           score = p->wgt_absorb_ * flux;
           if (i_nuclide >= 0) {
             score *= atom_density *
-              nuc_xs.get_xs(MG_GET_XS_PROMPT_NU_FISSION, p_g)
-              / macro_xs.get_xs(MG_GET_XS_ABSORPTION, p_g);
+              nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g)
+              / macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
           } else {
-            score *= macro_xs.get_xs(MG_GET_XS_PROMPT_NU_FISSION, p_g)
-              / macro_xs.get_xs(MG_GET_XS_ABSORPTION, p_g);
+            score *= macro_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g)
+              / macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
           }
         } else {
           // Skip any non-fission events
@@ -1668,23 +1668,23 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
           auto prompt_frac = 1. - n_delayed / static_cast<double>(p->n_bank_);
           score = simulation::keff * p->wgt_bank_ * prompt_frac * flux;
           if (i_nuclide >= 0) {
-            score *= atom_density * nuc_xs.get_xs(MG_GET_XS_FISSION, p_g)
-              / macro_xs.get_xs(MG_GET_XS_FISSION, p_g);
+            score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g)
+              / macro_xs.get_xs(MgxsType::FISSION, p_g);
           }
         }
       } else {
         if (i_nuclide >= 0) {
           score = atom_density * flux *
-            nuc_xs.get_xs(MG_GET_XS_PROMPT_NU_FISSION, p_g);
+            nuc_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g);
         } else {
-          score = flux * macro_xs.get_xs(MG_GET_XS_PROMPT_NU_FISSION, p_g);
+          score = flux * macro_xs.get_xs(MgxsType::PROMPT_NU_FISSION, p_g);
         }
       }
       break;
 
 
     case SCORE_DELAYED_NU_FISSION:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      if (tally.estimator_ == TallyEstimator::ANALOG) {
         if (settings::survival_biasing || p->fission_) {
           if (tally.energyout_filter_ != C_NONE) {
             // Fission has multiple outgoing neutrons so this helper function
@@ -1697,7 +1697,7 @@ score_general_mg(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
-          double abs_xs = macro_xs.get_xs(MG_GET_XS_ABSORPTION, p_g);
+          double abs_xs = macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
           if (abs_xs > 0.) {
             if (tally.delayedgroup_filter_ != C_NONE) {
               auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_];
@@ -1709,11 +1709,11 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
                 auto d = filt.groups()[d_bin] - 1;
                 score = p->wgt_absorb_ * flux;
                 if (i_nuclide >= 0) {
-                  score *= nuc_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g,
+                  score *= nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
                                          nullptr, nullptr, &d)
                     / abs_xs;
                 } else {
-                  score *= macro_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g,
+                  score *= macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
                                            nullptr, nullptr, &d)
                     / abs_xs;
                 }
@@ -1726,10 +1726,10 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
               // delayed-nu-fission xs to the absorption xs
               score = p->wgt_absorb_ * flux;
               if (i_nuclide >= 0) {
-                score *= nuc_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g)
+                score *= nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g)
                   / abs_xs;
               } else {
-                score *= macro_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g)
+                score *= macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g)
                   / abs_xs;
               }
             }
@@ -1755,8 +1755,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
               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(MG_GET_XS_FISSION, p_g)
-                  / macro_xs.get_xs(MG_GET_XS_FISSION, p_g);
+                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);
             }
@@ -1768,8 +1768,8 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
             score = simulation::keff * p->wgt_bank_ / p->n_bank_ * n_delayed
               * flux;
             if (i_nuclide >= 0) {
-              score *= atom_density * nuc_xs.get_xs(MG_GET_XS_FISSION, p_g)
-                / macro_xs.get_xs(MG_GET_XS_FISSION, p_g);
+              score *= atom_density * nuc_xs.get_xs(MgxsType::FISSION, p_g)
+                / macro_xs.get_xs(MgxsType::FISSION, p_g);
             }
           }
         }
@@ -1784,10 +1784,10 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
             auto d = filt.groups()[d_bin] - 1;
             if (i_nuclide >= 0) {
               score = flux * atom_density * nuc_xs.get_xs(
-                MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d);
+                MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d);
             } else {
               score = flux * macro_xs.get_xs(
-                MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d);
+                MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d);
             }
             score_fission_delayed_dg(i_tally, d_bin, score, score_index);
           }
@@ -1795,9 +1795,9 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
         } else {
           if (i_nuclide >= 0) {
             score = flux * atom_density *
-              nuc_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g);
+              nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g);
           } else {
-            score = flux * macro_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g);
+            score = flux * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g);
           }
         }
       }
@@ -1805,12 +1805,12 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
 
 
     case SCORE_DECAY_RATE:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      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
-          double abs_xs = macro_xs.get_xs(MG_GET_XS_ABSORPTION, p_g);
+          double abs_xs = macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
           if (abs_xs > 0) {
             if (tally.delayedgroup_filter_ != C_NONE) {
               auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_];
@@ -1822,14 +1822,14 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
                 auto d = filt.groups()[d_bin] - 1;
                 score = p->wgt_absorb_ * flux;
                 if (i_nuclide >= 0) {
-                  score *= nuc_xs.get_xs(MG_GET_XS_DECAY_RATE, p_g,
+                  score *= nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g,
                                          nullptr, nullptr, &d)
-                    * nuc_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g,
+                    * nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
                                     nullptr, nullptr, &d) / abs_xs;
                 } else {
-                  score *= macro_xs.get_xs(MG_GET_XS_DECAY_RATE, p_g, nullptr,
+                  score *= macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr,
                                            nullptr, &d)
-                    * macro_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g,
+                    * 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);
@@ -1844,15 +1844,15 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
               for (auto d = 0; d < data::mg.num_delayed_groups_; ++d) {
                 if (i_nuclide >= 0) {
                   score += p->wgt_absorb_ * flux
-                    * nuc_xs.get_xs(MG_GET_XS_DECAY_RATE, p_g, nullptr,
+                    * nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr,
                                     nullptr, &d)
-                    * nuc_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g,
+                    * nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
                                     nullptr, nullptr, &d) / abs_xs;
                 } else {
                   score += p->wgt_absorb_ * flux
-                    * macro_xs.get_xs(MG_GET_XS_DECAY_RATE, p_g, nullptr,
+                    * macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr,
                                       nullptr, &d)
-                    * macro_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g,
+                    * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g,
                                       nullptr, nullptr, &d) / abs_xs;
                 }
               }
@@ -1876,12 +1876,12 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
             if (d != -1) {
               if (i_nuclide >= 0) {
                 score += simulation::keff * atom_density * bank.wgt * flux
-                  * nuc_xs.get_xs(MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d)
-                  * nuc_xs.get_xs(MG_GET_XS_FISSION, p_g)
-                  / macro_xs.get_xs(MG_GET_XS_FISSION, p_g);
+                  * nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d)
+                  * nuc_xs.get_xs(MgxsType::FISSION, p_g)
+                  / macro_xs.get_xs(MgxsType::FISSION, p_g);
               } else {
                 score += simulation::keff * bank.wgt * flux
-                  * macro_xs.get_xs(MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d);
+                  * macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d);
               }
               if (tally.delayedgroup_filter_ != C_NONE) {
                 auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_];
@@ -1912,15 +1912,15 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
             auto d = filt.groups()[d_bin] - 1;
             if (i_nuclide >= 0) {
               score += atom_density * flux
-                * nuc_xs.get_xs(MG_GET_XS_DECAY_RATE, p_g, nullptr,
+                * nuc_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr,
                                 nullptr, &d)
-                * nuc_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr,
+                * nuc_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr,
                                 nullptr, &d);
             } else {
               score += flux
-                * macro_xs.get_xs(MG_GET_XS_DECAY_RATE, p_g, nullptr,
+                * macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr,
                                   nullptr, &d)
-                * macro_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr,
+                * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr,
                                   nullptr, &d);
             }
             score_fission_delayed_dg(i_tally, d_bin, score, score_index);
@@ -1931,12 +1931,12 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
           for (auto d = 0; d < data::mg.num_delayed_groups_; ++d) {
             if (i_nuclide >= 0) {
               score += atom_density * flux
-                * nuc_xs.get_xs(MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d)
-                * nuc_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d);
+                * 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);
             } else {
               score += flux
-                * macro_xs.get_xs(MG_GET_XS_DECAY_RATE, p_g, nullptr, nullptr, &d)
-                * macro_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d);
+                * macro_xs.get_xs(MgxsType::DECAY_RATE, p_g, nullptr, nullptr, &d)
+                * macro_xs.get_xs(MgxsType::DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d);
             }
           }
         }
@@ -1945,7 +1945,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
 
 
     case SCORE_KAPPA_FISSION:
-      if (tally.estimator_ == ESTIMATOR_ANALOG) {
+      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
@@ -1953,7 +1953,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
           score = p->wgt_absorb_ * flux;
         } else {
           // Skip any non-absorption events
-          if (p->event_ == EVENT_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
@@ -1961,18 +1961,18 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
         }
         if (i_nuclide >= 0) {
           score *= atom_density
-            * nuc_xs.get_xs(MG_GET_XS_KAPPA_FISSION, p_g)
-            / macro_xs.get_xs(MG_GET_XS_ABSORPTION, p_g);
+            * nuc_xs.get_xs(MgxsType::KAPPA_FISSION, p_g)
+            / macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
         } else {
           score *=
-            macro_xs.get_xs(MG_GET_XS_KAPPA_FISSION, p_g)
-            / macro_xs.get_xs(MG_GET_XS_ABSORPTION, p_g);
+            macro_xs.get_xs(MgxsType::KAPPA_FISSION, p_g)
+            / macro_xs.get_xs(MgxsType::ABSORPTION, p_g);
         }
       } else {
         if (i_nuclide >= 0) {
-          score = atom_density * flux * nuc_xs.get_xs(MG_GET_XS_KAPPA_FISSION, p_g);
+          score = atom_density * flux * nuc_xs.get_xs(MgxsType::KAPPA_FISSION, p_g);
         } else {
-          score = flux * macro_xs.get_xs(MG_GET_XS_KAPPA_FISSION, p_g);
+          score = flux * macro_xs.get_xs(MgxsType::KAPPA_FISSION, p_g);
         }
       }
       break;
@@ -1989,7 +1989,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index,
 
     // Update tally results
     #pragma omp atomic
-    tally.results_(filter_index, score_index, RESULT_VALUE) += score;
+    tally.results_(filter_index, score_index, TallyResult::VALUE) += score;
   }
 }
 
@@ -2289,7 +2289,7 @@ score_surface_tally(const Particle* p, const std::vector<int>& tallies)
       for (auto score_index = 0; score_index < tally.scores_.size();
            ++score_index) {
         #pragma omp atomic
-        tally.results_(filter_index, score_index, RESULT_VALUE) += score;
+        tally.results_(filter_index, score_index, TallyResult::VALUE) += score;
       }
     }
 
diff --git a/src/tallies/trigger.cpp b/src/tallies/trigger.cpp
index ebf3d83bb..7fb3475c9 100644
--- a/src/tallies/trigger.cpp
+++ b/src/tallies/trigger.cpp
@@ -31,8 +31,8 @@ get_tally_uncertainty(int i_tally, int score_index, int filter_index)
 {
   const auto& tally {model::tallies[i_tally]};
 
-  auto sum = tally->results_(filter_index, score_index, RESULT_SUM);
-  auto sum_sq = tally->results_(filter_index, score_index, RESULT_SUM_SQ);
+  auto sum = tally->results_(filter_index, score_index, TallyResult::SUM);
+  auto sum_sq = tally->results_(filter_index, score_index, TallyResult::SUM_SQ);
 
   int n = tally->n_realizations_;
   auto mean = sum / n;
@@ -111,7 +111,7 @@ check_tally_triggers(double& ratio, int& tally_id, int& score)
 double
 check_keff_trigger()
 {
-  if (settings::run_mode != RUN_MODE_EIGENVALUE) return 0.0;
+  if (settings::run_mode != RunMode::EIGENVALUE) return 0.0;
 
   double k_combined[2];
   openmc_get_keff(k_combined);
diff --git a/src/thermal.cpp b/src/thermal.cpp
index 5b0ee0d20..ab0dd1542 100644
--- a/src/thermal.cpp
+++ b/src/thermal.cpp
@@ -76,7 +76,7 @@ ThermalScattering::ThermalScattering(hid_t group, const std::vector<double>& tem
   }
 
   switch (settings::temperature_method) {
-  case TEMPERATURE_NEAREST:
+  case TemperatureInterpolationType::NEAREST:
     // Determine actual temperatures to read
     for (const auto& T : temperature) {
 
@@ -96,7 +96,7 @@ ThermalScattering::ThermalScattering(hid_t group, const std::vector<double>& tem
     }
     break;
 
-  case TEMPERATURE_INTERPOLATION:
+  case TemperatureInterpolationType::INTERPOLATION:
     // If temperature interpolation or multipole is selected, get a list of
     // bounding temperatures for each actual temperature present in the model
     for (const auto& T : temperature) {
@@ -156,7 +156,7 @@ ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp,
   // Determine temperature for S(a,b) table
   double kT = sqrtkT*sqrtkT;
   int i;
-  if (settings::temperature_method == TEMPERATURE_NEAREST) {
+  if (settings::temperature_method == TemperatureInterpolationType::NEAREST) {
     // If using nearest temperature, do linear search on temperature
     for (i = 0; i < kTs_.size(); ++i) {
       if (std::abs(kTs_[i] - kT) < K_BOLTZMANN*settings::temperature_tolerance) {
diff --git a/src/volume_calc.cpp b/src/volume_calc.cpp
index 903a111a2..10ea6bca4 100644
--- a/src/volume_calc.cpp
+++ b/src/volume_calc.cpp
@@ -45,11 +45,11 @@ VolumeCalculation::VolumeCalculation(pugi::xml_node node)
   // Read domain type (cell, material or universe)
   std::string domain_type = get_node_value(node, "domain_type");
   if (domain_type == "cell") {
-    domain_type_ = FILTER_CELL;
+    domain_type_ = TallyDomain::CELL;
   } else if (domain_type == "material") {
-    domain_type_ = FILTER_MATERIAL;
+    domain_type_ = TallyDomain::MATERIAL;
   } else if (domain_type == "universe") {
-    domain_type_ = FILTER_UNIVERSE;
+    domain_type_ = TallyDomain::UNIVERSE;
   } else {
     fatal_error(std::string("Unrecognized domain type for stochastic "
       "volume calculation: " + domain_type));
@@ -139,7 +139,7 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
         // If this location is not in the geometry at all, move on to next block
         if (!find_cell(&p, false)) continue;
 
-        if (domain_type_ == FILTER_MATERIAL) {
+        if (domain_type_ == TallyDomain::MATERIAL) {
           if (p.material_ != MATERIAL_VOID) {
             for (int i_domain = 0; i_domain < n; i_domain++) {
               if (model::materials[p.material_]->id_ == domain_ids_[i_domain]) {
@@ -148,7 +148,7 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
               }
             }
           }
-        } else if (domain_type_ == FILTER_CELL) {
+        } else if (domain_type_ == TallyDomain::CELL) {
           for (int level = 0; level < p.n_coord_; ++level) {
             for (int i_domain=0; i_domain < n; i_domain++) {
               if (model::cells[p.coord_[level].cell]->id_ == domain_ids_[i_domain]) {
@@ -157,7 +157,7 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
               }
             }
           }
-        } else if (domain_type_ == FILTER_UNIVERSE) {
+        } else if (domain_type_ == TallyDomain::UNIVERSE) {
           for (int level = 0; level < p.n_coord_; ++level) {
             for (int i_domain = 0; i_domain < n; ++i_domain) {
               if (model::universes[p.coord_[level].universe]->id_ == domain_ids_[i_domain]) {
@@ -382,13 +382,13 @@ void VolumeCalculation::to_hdf5(const std::string& filename,
     write_attribute(file_id, "iterations", 1);
   }
 
-  if (domain_type_ == FILTER_CELL) {
+  if (domain_type_ == TallyDomain::CELL) {
     write_attribute(file_id, "domain_type", "cell");
   }
-  else if (domain_type_ == FILTER_MATERIAL) {
+  else if (domain_type_ == TallyDomain::MATERIAL) {
     write_attribute(file_id, "domain_type", "material");
   }
-  else if (domain_type_ == FILTER_UNIVERSE) {
+  else if (domain_type_ == TallyDomain::UNIVERSE) {
     write_attribute(file_id, "domain_type", "universe");
   }
 
@@ -477,9 +477,9 @@ int openmc_calculate_volumes() {
 
     if (mpi::master) {
       std::string domain_type;
-      if (vol_calc.domain_type_ == FILTER_CELL) {
+      if (vol_calc.domain_type_ == VolumeCalculation::TallyDomain::CELL) {
         domain_type = "  Cell ";
-      } else if (vol_calc.domain_type_ == FILTER_MATERIAL) {
+      } else if (vol_calc.domain_type_ == VolumeCalculation::TallyDomain::MATERIAL) {
         domain_type = "  Material ";
       } else {
         domain_type = "  Universe ";
diff --git a/src/xsdata.cpp b/src/xsdata.cpp
index d6b92958a..c39d779ec 100644
--- a/src/xsdata.cpp
+++ b/src/xsdata.cpp
@@ -24,7 +24,7 @@ namespace openmc {
 // XsData class methods
 //==============================================================================
 
-XsData::XsData(bool fissionable, int scatter_format, int n_pol, int n_azi,
+XsData::XsData(bool fissionable, AngleDistributionType scatter_format, int n_pol, int n_azi,
                size_t n_groups, size_t n_d_groups) :
   n_g_(n_groups),
   n_dg_(n_d_groups)
@@ -32,8 +32,8 @@ XsData::XsData(bool fissionable, int scatter_format, int n_pol, int n_azi,
   size_t n_ang = n_pol * n_azi;
 
   // check to make sure scatter format is OK before we allocate
-  if (scatter_format != ANGLE_HISTOGRAM && scatter_format != ANGLE_TABULAR &&
-      scatter_format != ANGLE_LEGENDRE) {
+  if (scatter_format != AngleDistributionType::HISTOGRAM && scatter_format != AngleDistributionType::TABULAR &&
+      scatter_format != AngleDistributionType::LEGENDRE) {
     fatal_error("Invalid scatter_format!");
   }
   // allocate all [temperature][angle][in group] quantities
@@ -68,11 +68,11 @@ XsData::XsData(bool fissionable, int scatter_format, int n_pol, int n_azi,
 
 
   for (int a = 0; a < n_ang; a++) {
-    if (scatter_format == ANGLE_HISTOGRAM) {
+    if (scatter_format == AngleDistributionType::HISTOGRAM) {
       scatter.emplace_back(new ScattDataHistogram);
-    } else if (scatter_format == ANGLE_TABULAR) {
+    } else if (scatter_format == AngleDistributionType::TABULAR) {
       scatter.emplace_back(new ScattDataTabular);
-    } else if (scatter_format == ANGLE_LEGENDRE) {
+    } else if (scatter_format == AngleDistributionType::LEGENDRE) {
       scatter.emplace_back(new ScattDataLegendre);
     }
   }
@@ -81,8 +81,8 @@ XsData::XsData(bool fissionable, int scatter_format, int n_pol, int n_azi,
 //==============================================================================
 
 void
-XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format,
-  int final_scatter_format, int order_data, bool is_isotropic, int n_pol, int n_azi)
+XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, AngleDistributionType scatter_format,
+  AngleDistributionType final_scatter_format, int order_data, bool is_isotropic, int n_pol, int n_azi)
 {
   // Reconstruct the dimension information so it doesn't need to be passed
   size_t n_ang = n_pol * n_azi;
@@ -398,8 +398,8 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_isotropic)
 //==============================================================================
 
 void
-XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
-     int scatter_format, int final_scatter_format, int order_data)
+XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, AngleDistributionType scatter_format,
+    AngleDistributionType final_scatter_format, int order_data)
 {
   if (!object_exists(xsdata_grp, "scatter_data")) {
     fatal_error("Must provide scatter_data group!");
@@ -427,7 +427,7 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
   // Compare the number of orders given with the max order of the problem;
   // strip off the superfluous orders if needed
   int order_dim;
-  if (scatter_format == ANGLE_LEGENDRE) {
+  if (scatter_format == AngleDistributionType::LEGENDRE) {
     order_dim = std::min(order_data - 1, settings::max_order) + 1;
   } else {
     order_dim = order_data;
@@ -480,8 +480,8 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
   close_group(scatt_grp);
 
   // Finally, convert the Legendre data to tabular, if needed
-  if (scatter_format == ANGLE_LEGENDRE &&
-      final_scatter_format == ANGLE_TABULAR) {
+  if (scatter_format == AngleDistributionType::LEGENDRE &&
+      final_scatter_format == AngleDistributionType::TABULAR) {
     for (size_t a = 0; a < n_ang; a++) {
       ScattDataLegendre legendre_scatt;
       xt::xtensor<int, 1> in_gmin = xt::view(gmin, a, xt::all());
diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py
index ad198f255..264fcafdd 100644
--- a/tests/unit_tests/test_lib.py
+++ b/tests/unit_tests/test_lib.py
@@ -170,10 +170,6 @@ def test_settings(lib_init):
     assert settings.seed == 1
     settings.seed = 11
 
-    assert settings.run_mode == 'eigenvalue'
-    settings.run_mode = 'volume'
-    settings.run_mode = 'eigenvalue'
-
 
 def test_tally_mapping(lib_init):
     tallies = openmc.lib.tallies