diff --git a/docs/source/devguide/statepoint.rst b/docs/source/devguide/statepoint.rst
index cd2785fc5..ffcc3d968 100644
--- a/docs/source/devguide/statepoint.rst
+++ b/docs/source/devguide/statepoint.rst
@@ -4,6 +4,621 @@
 State Point Binary File Specifications
 ======================================
 
+----------
+Revision 5
+----------
+
+**integer(4) REVISION_STATEPOINT**
+
+    Revision of the binary state point file. Any time a change is made in the
+    format of the state-point file, this integer is incremented.
+
+**integer(4) VERSION_MAJOR**
+
+    Major version number for OpenMC
+
+**integer(4) VERSION_MINOR**
+
+    Minor version number for OpenMC
+
+**integer(4) VERSION_RELEASE**
+
+    Release version number for OpenMC
+
+**character(19) time_stamp**
+
+    Date and time the state point was written.
+
+**integer(8) seed**
+
+    Pseudo-random number generator seed.
+
+**integer(4) run_mode**
+
+    run mode used. The modes are described in constants.F90.
+
+**integer(8) n_particles**
+
+    Number of particles used per generation.
+
+**integer(4) n_batches**
+
+    Total number of batches (active + inactive).
+
+**integer(4) current_batch**
+
+    The number of batches already simulated.
+
+if (run_mode == MODE_CRITICALITY)
+
+    **integer(4) n_inactive**
+
+        Number of inactive batches
+
+    **integer(4) gen_per_batch**
+
+        Number of generations per batch for criticality calculations
+
+    *do i = 1, current_batch*
+
+        **real(8) k_batch(i)**
+
+             k-effective for the i-th batch
+
+    *do i = 1, current_batch*
+
+        **real(8) entropy(i)**
+
+            Shannon entropy for the i-th batch
+
+**integer(4) n_meshes**
+
+    Number of meshes in tallies.xml file
+
+*do i = 1, n_meshes*
+
+    **integer(4) meshes(i) % type**
+
+        Type of mesh.
+
+    **integer(4) meshes(i) % n_dimension**
+
+        Number of dimensions for mesh (2 or 3).
+
+    **integer(4) meshes(i) % dimension(:)**
+
+        Number of mesh cells in each dimension.
+
+    **real(8) meshes(i) % lower_left(:)**
+
+        Coordinates of lower-left corner of mesh.
+
+    **real(8) meshes(i) % upper_right(:)**
+
+        Coordinates of upper-right corner of mesh.
+
+    **real(8) meshes(i) % width(:)**
+
+        Width of each mesh cell in each dimension.
+
+**integer(4) n_tallies**
+
+*do i = 1, n_tallies*
+
+    **integer(4) tallies(i) % n_realizations**
+
+        Number of realizations for the i-th tally.
+
+    **integer(4) size(tallies(i) % scores, 1)**
+
+        Total number of score bins for the i-th tally
+
+    **integer(4) size(tallies(i) % scores, 2)**
+
+        Total number of filter bins for the i-th tally
+
+    **integer(4) tallies(i) % n_filters**
+
+    *do j = 1, tallies(i) % n_filters*
+
+        **integer(4) tallies(i) % filter(j) % type**
+
+            Type of tally filter.
+
+        **integer(4) tallies(i) % filter(j) % n_bins**
+
+            Number of bins for filter.
+
+        **integer(4)/real(8) tallies(i) % filter(j) % bins(:)**
+
+            Value for each filter bin of this type.
+
+    **integer(4) tallies(i) % n_nuclide_bins**
+
+        Number of nuclide bins. If none are specified, this is just one.
+
+    *do j = 1, tallies(i) % n_nuclide_bins*
+
+        **integer(4) tallies(i) % nuclide_bins(j)**
+
+            Values of specified nuclide bins
+
+    **integer(4) tallies(i) % n_score_bins**
+
+        Number of scoring bins.
+
+    *do j = 1, tallies(i) % n_score_bins*
+
+        **integer(4) tallies(i) % score_bins(j)**
+
+            Values of specified scoring bins (e.g. SCORE_FLUX).
+
+**integer(4) n_realizations**
+
+    Number of realizations for global tallies.
+
+**integer(4) N_GLOBAL_TALLIES**
+
+    Number of global tally scores
+
+*do i = 1, N_GLOBAL_TALLIES*
+
+    **real(8) global_tallies(i) % sum**
+
+        Accumulated sum for the i-th global tally
+
+    **real(8) global_tallies(i) % sum_sq**
+
+        Accumulated sum of squares for the i-th global tally
+
+**integer(4) tallies_on**
+
+    Flag indicated if tallies are present in the file.
+
+if (tallies_on > 0)
+
+    *do i = 1, n_tallies*
+
+        *do k = 1, size(tallies(i) % scores, 2)*
+
+            *do j = 1, size(tallies(i) % scores, 1)*
+
+                **real(8) tallies(i) % scores(j,k) % sum**
+            
+                    Accumulated sum for the j-th score and k-th filter of the
+                    i-th tally
+
+                **real(8) tallies(i) % scores(j,k) % sum_sq**
+
+                    Accumulated sum of squares for the j-th score and k-th
+                    filter of the i-th tally
+
+if (run_mode == MODE_CRITICALITY)
+
+    *do i = 1, n_particles*
+
+        **real(8) source_bank(i) % wgt**
+
+            Weight of the i-th source particle
+
+        **real(8) source_bank(i) % xyz(1:3)**
+
+            Coordinates of the i-th source particle.
+
+        **real(8) source_bank(i) % uvw(1:3)**
+
+            Direction of the i-th source particle
+
+        **real(8) source_bank(i) % E**
+
+            Energy of the i-th source particle.
+
+----------
+Revision 4
+----------
+
+**integer(4) REVISION_STATEPOINT**
+
+    Revision of the binary state point file. Any time a change is made in the
+    format of the state-point file, this integer is incremented.
+
+**integer(4) VERSION_MAJOR**
+
+    Major version number for OpenMC
+
+**integer(4) VERSION_MINOR**
+
+    Minor version number for OpenMC
+
+**integer(4) VERSION_RELEASE**
+
+    Release version number for OpenMC
+
+**character(19) time_stamp**
+
+    Date and time the state point was written.
+
+**integer(8) seed**
+
+    Pseudo-random number generator seed.
+
+**integer(4) run_mode**
+
+    run mode used. The modes are described in constants.F90.
+
+**integer(8) n_particles**
+
+    Number of particles used per generation.
+
+**integer(4) n_batches**
+
+    Total number of batches (active + inactive).
+
+**integer(4) current_batch**
+
+    The number of batches already simulated.
+
+if (run_mode == MODE_CRITICALITY)
+
+    **integer(4) n_inactive**
+
+        Number of inactive batches
+
+    **integer(4) gen_per_batch**
+
+        Number of generations per batch for criticality calculations
+
+    *do i = 1, current_batch*
+
+        **real(8) k_batch(i)**
+
+             k-effective for the i-th batch
+
+    *do i = 1, current_batch*
+
+        **real(8) entropy(i)**
+
+            Shannon entropy for the i-th batch
+
+**integer(4) n_meshes**
+
+    Number of meshes in tallies.xml file
+
+*do i = 1, n_meshes*
+
+    **integer(4) meshes(i) % type**
+
+        Type of mesh.
+
+    **integer(4) meshes(i) % n_dimension**
+
+        Number of dimensions for mesh (2 or 3).
+
+    **integer(4) meshes(i) % dimension(:)**
+
+        Number of mesh cells in each dimension.
+
+    **real(8) meshes(i) % lower_left(:)**
+
+        Coordinates of lower-left corner of mesh.
+
+    **real(8) meshes(i) % upper_right(:)**
+
+        Coordinates of upper-right corner of mesh.
+
+    **real(8) meshes(i) % width(:)**
+
+        Width of each mesh cell in each dimension.
+
+**integer(4) n_tallies**
+
+*do i = 1, n_tallies*
+
+    **integer(4) size(tallies(i) % scores, 1)**
+
+        Total number of score bins for the i-th tally
+
+    **integer(4) size(tallies(i) % scores, 2)**
+
+        Total number of filter bins for the i-th tally
+
+    **integer(4) tallies(i) % n_filters**
+
+    *do j = 1, tallies(i) % n_filters*
+
+        **integer(4) tallies(i) % filter(j) % type**
+
+            Type of tally filter.
+
+        **integer(4) tallies(i) % filter(j) % n_bins**
+
+            Number of bins for filter.
+
+        **integer(4)/real(8) tallies(i) % filter(j) % bins(:)**
+
+            Value for each filter bin of this type.
+
+    **integer(4) tallies(i) % n_nuclide_bins**
+
+        Number of nuclide bins. If none are specified, this is just one.
+
+    *do j = 1, tallies(i) % n_nuclide_bins*
+
+        **integer(4) tallies(i) % nuclide_bins(j)**
+
+            Values of specified nuclide bins
+
+    **integer(4) tallies(i) % n_score_bins**
+
+        Number of scoring bins.
+
+    *do j = 1, tallies(i) % n_score_bins*
+
+        **integer(4) tallies(i) % score_bins(j)**
+
+            Values of specified scoring bins (e.g. SCORE_FLUX).
+
+**integer(4) N_GLOBAL_TALLIES**
+
+    Number of global tally scores
+
+*do i = 1, N_GLOBAL_TALLIES*
+
+    **real(8) global_tallies(i) % sum**
+
+        Accumulated sum for the i-th global tally
+
+    **real(8) global_tallies(i) % sum_sq**
+
+        Accumulated sum of squares for the i-th global tally
+
+**integer(4) tallies_on**
+
+    Flag indicated if tallies are present in the file.
+
+if (tallies_on > 0)
+
+    **integer(4) n_realizations**
+
+        Number of realizations for tally random variables.
+
+    *do i = 1, n_tallies*
+
+        *do k = 1, size(tallies(i) % scores, 2)*
+
+            *do j = 1, size(tallies(i) % scores, 1)*
+
+                **real(8) tallies(i) % scores(j,k) % sum**
+            
+                    Accumulated sum for the j-th score and k-th filter of the
+                    i-th tally
+
+                **real(8) tallies(i) % scores(j,k) % sum_sq**
+
+                    Accumulated sum of squares for the j-th score and k-th
+                    filter of the i-th tally
+
+if (run_mode == MODE_CRITICALITY)
+
+    *do i = 1, n_particles*
+
+        **real(8) source_bank(i) % wgt**
+
+            Weight of the i-th source particle
+
+        **real(8) source_bank(i) % xyz(1:3)**
+
+            Coordinates of the i-th source particle.
+
+        **real(8) source_bank(i) % uvw(1:3)**
+
+            Direction of the i-th source particle
+
+        **real(8) source_bank(i) % E**
+
+            Energy of the i-th source particle.
+
+----------
+Revision 3
+----------
+
+**integer(4) REVISION_STATEPOINT**
+
+    Revision of the binary state point file. Any time a change is made in the
+    format of the state-point file, this integer is incremented.
+
+**integer(4) VERSION_MAJOR**
+
+    Major version number for OpenMC
+
+**integer(4) VERSION_MINOR**
+
+    Minor version number for OpenMC
+
+**integer(4) VERSION_RELEASE**
+
+    Release version number for OpenMC
+
+**character(19) time_stamp**
+
+    Date and time the state point was written.
+
+**integer(8) seed**
+
+    Pseudo-random number generator seed.
+
+**integer(4) run_mode**
+
+    run mode used. The modes are described in constants.F90.
+
+**integer(8) n_particles**
+
+    Number of particles used per generation.
+
+**integer(4) n_batches**
+
+    Total number of batches (active + inactive).
+
+**integer(4) current_batch**
+
+    The number of batches already simulated.
+
+if (run_mode == MODE_CRITICALITY)
+
+    **integer(4) n_inactive**
+
+        Number of inactive batches
+
+    **integer(4) gen_per_batch**
+
+        Number of generations per batch for criticality calculations
+
+    *do i = 1, current_batch*
+
+        **real(8) k_batch(i)**
+
+             k-effective for the i-th batch
+
+    *do i = 1, current_batch*
+
+        **real(8) entropy(i)**
+
+            Shannon entropy for the i-th batch
+
+**integer(4) N_GLOBAL_TALLIES**
+
+    Number of global tally scores
+
+*do i = 1, N_GLOBAL_TALLIES*
+
+    **real(8) global_tallies(i) % sum**
+
+        Accumulated sum for the i-th global tally
+
+    **real(8) global_tallies(i) % sum_sq**
+
+        Accumulated sum of squares for the i-th global tally
+
+**integer(4) n_meshes**
+
+    Number of meshes in tallies.xml file
+
+*do i = 1, n_meshes*
+
+    **integer(4) meshes(i) % type**
+
+        Type of mesh.
+
+    **integer(4) meshes(i) % n_dimension**
+
+        Number of dimensions for mesh (2 or 3).
+
+    **integer(4) meshes(i) % dimension(:)**
+
+        Number of mesh cells in each dimension.
+
+    **real(8) meshes(i) % lower_left(:)**
+
+        Coordinates of lower-left corner of mesh.
+
+    **real(8) meshes(i) % upper_right(:)**
+
+        Coordinates of upper-right corner of mesh.
+
+    **real(8) meshes(i) % width(:)**
+
+        Width of each mesh cell in each dimension.
+
+**integer(4) n_tallies**
+
+*do i = 1, n_tallies*
+
+    **integer(4) size(tallies(i) % scores, 1)**
+
+        Total number of score bins for the i-th tally
+
+    **integer(4) size(tallies(i) % scores, 2)**
+
+        Total number of filter bins for the i-th tally
+
+    **integer(4) tallies(i) % n_filters**
+
+    *do j = 1, tallies(i) % n_filters*
+
+        **integer(4) tallies(i) % filter(j) % type**
+
+            Type of tally filter.
+
+        **integer(4) tallies(i) % filter(j) % n_bins**
+
+            Number of bins for filter.
+
+        **integer(4)/real(8) tallies(i) % filter(j) % bins(:)**
+
+            Value for each filter bin of this type.
+
+    **integer(4) tallies(i) % n_nuclide_bins**
+
+        Number of nuclide bins. If none are specified, this is just one.
+
+    *do j = 1, tallies(i) % n_nuclide_bins*
+
+        **integer(4) tallies(i) % nuclide_bins(j)**
+
+            Values of specified nuclide bins
+
+    **integer(4) tallies(i) % n_score_bins**
+
+        Number of scoring bins.
+
+    *do j = 1, tallies(i) % n_score_bins*
+
+        **integer(4) tallies(i) % score_bins(j)**
+
+            Values of specified scoring bins (e.g. SCORE_FLUX).
+
+**integer(4) tallies_on**
+
+    Flag indicated if tallies are present in the file.
+
+if (tallies_on > 0)
+
+    *do i = 1, n_tallies*
+
+        *do k = 1, size(tallies(i) % scores, 2)*
+
+            *do j = 1, size(tallies(i) % scores, 1)*
+
+                **real(8) tallies(i) % scores(j,k) % sum**
+            
+                    Accumulated sum for the j-th score and k-th filter of the
+                    i-th tally
+
+                **real(8) tallies(i) % scores(j,k) % sum_sq**
+
+                    Accumulated sum of squares for the j-th score and k-th
+                    filter of the i-th tally
+
+if (run_mode == MODE_CRITICALITY)
+
+    *do i = 1, n_particles*
+
+        **real(8) source_bank(i) % wgt**
+
+            Weight of the i-th source particle
+
+        **real(8) source_bank(i) % xyz(1:3)**
+
+            Coordinates of the i-th source particle.
+
+        **real(8) source_bank(i) % uvw(1:3)**
+
+            Direction of the i-th source particle
+
+        **real(8) source_bank(i) % E**
+
+            Energy of the i-th source particle.
+
 ----------
 Revision 2
 ----------
@@ -57,7 +672,9 @@ Revision 2
 
     if (entropy_on)
 
-        **real(8) entropy for the i-th batch**
+        **real(8) entropy(i)**
+
+            Shannon entropy for the i-th batch
 
 **integer(4) N_GLOBAL_TALLIES**
 
@@ -156,7 +773,9 @@ Revision 1
 
     if (entropy_on)
 
-        **real(8) entropy for the i-th batch**
+        **real(8) entropy(i)**
+
+            Shannon entropy for the i-th batch
 
 **integer(4) N_GLOBAL_TALLIES**
 
diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90
index 4942a025c..56b7c4dc8 100644
--- a/src/hdf5_interface.F90
+++ b/src/hdf5_interface.F90
@@ -929,6 +929,10 @@ contains
        call h5gcreate_f(tallies_group, "tally " // to_str(t % id), &
             temp_group, hdf5_err)
 
+       ! Write number of realizations
+       call hdf5_make_integer(temp_group, "n_realizations", &
+            t % n_realizations)
+
        ! Write size of each tally
        call hdf5_make_integer(temp_group, "n_filter_bins", size(t % scores, 1))
        call hdf5_make_integer(temp_group, "n_total_score_bins", &
@@ -1036,6 +1040,9 @@ contains
        call h5gclose_f(temp_group, hdf5_err)
     end do TALLY_METADATA
 
+    ! Write number of realizations for global tallies
+    call hdf5_make_integer(hdf5_state_point, "n_realizations", n_realizations)
+
     ! Write out global tallies sum and sum_sq
     call hdf5_make_integer(hdf5_state_point, "n_global_tallies", &
          N_GLOBAL_TALLIES)
@@ -1054,9 +1061,6 @@ contains
        ! Indicate that tallies are on
        call hdf5_make_integer(tallies_group, "tallies_present", 1)
 
-       ! Write number of realizations
-       call hdf5_make_integer(tallies_group, "n_realizations", n_realizations)
-
        ! Write tally sum and sum_sq
        TALLY_SCORES: do i = 1, n_tallies
           ! Get pointer to tally
diff --git a/src/state_point.F90 b/src/state_point.F90
index c7283e929..6e1b03b32 100644
--- a/src/state_point.F90
+++ b/src/state_point.F90
@@ -67,6 +67,10 @@ contains
        call write_tally_scores_nr(fh)
 
     elseif (master) then
+       ! Write number of realizations
+       call MPI_FILE_WRITE(fh, n_realizations, 1, MPI_INTEGER, &
+            MPI_STATUS_IGNORE, mpi_err)
+
        ! Write global tallies
        call MPI_FILE_WRITE(fh, N_GLOBAL_TALLIES, 1, MPI_INTEGER, &
             MPI_STATUS_IGNORE, mpi_err)
@@ -79,10 +83,6 @@ contains
           call MPI_FILE_WRITE(fh, temp, 1, MPI_INTEGER, &
                MPI_STATUS_IGNORE, mpi_err)
 
-          ! Write number of realizations
-          call MPI_FILE_WRITE(fh, n_realizations, 1, MPI_INTEGER, &
-               MPI_STATUS_IGNORE, mpi_err)
-
           ! Write all tally scores
           TALLY_SCORES: do i = 1, n_tallies
              t => tallies(i)
@@ -188,6 +188,9 @@ contains
        ! Get pointer to tally
        t => tallies(i)
 
+       ! Number of realizations
+       write(UNIT_STATE) t % n_realizations
+
        ! Write size of each tally
        write(UNIT_STATE) size(t % scores, 1)
        write(UNIT_STATE) size(t % scores, 2)
@@ -240,6 +243,9 @@ contains
        write(UNIT_STATE) t % score_bins
     end do TALLY_METADATA
 
+    ! Number of realizations for global tallies
+    write(UNIT_STATE) n_realizations
+
     ! Write out global tallies sum and sum_sq
     write(UNIT_STATE) N_GLOBAL_TALLIES
     GLOBAL_TALLIES_LOOP: do i = 1, N_GLOBAL_TALLIES
@@ -251,9 +257,6 @@ contains
        ! Indicate that tallies are on
        write(UNIT_STATE) 1
 
-       ! Number of realizations
-       write(UNIT_STATE) n_realizations
-
        TALLY_SCORES: do i = 1, n_tallies
           ! Get pointer to tally
           t => tallies(i)
@@ -382,6 +385,10 @@ contains
        ! Get pointer to tally
        t => tallies(i)
 
+       ! Write number of realizations
+       call MPI_FILE_WRITE(fh, t % n_realizations, 1, MPI_INTEGER, &
+            MPI_STATUS_IGNORE, mpi_err)
+
        ! Write size of each tally
        n = t % n_score_bins * t % n_nuclide_bins
        call MPI_FILE_WRITE(fh, n, 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err)
@@ -477,8 +484,12 @@ contains
     ! ==========================================================================
     ! COLLECT AND WRITE GLOBAL TALLIES
 
-    ! Write number of global tallies
     if (master) then
+       ! Write number of realizations
+       call MPI_FILE_WRITE(fh, n_realizations, 1, MPI_INTEGER, &
+            MPI_STATUS_IGNORE, mpi_err)
+
+       ! Write number of global tallies
        call MPI_FILE_WRITE(fh, N_GLOBAL_TALLIES, 1, MPI_INTEGER, &
             MPI_STATUS_IGNORE, mpi_err)
     end if
@@ -515,10 +526,6 @@ contains
           temp = 1
           call MPI_FILE_WRITE(fh, temp, 1, MPI_INTEGER, &
                MPI_STATUS_IGNORE, mpi_err)
-
-          ! Write number of realizations
-          call MPI_FILE_WRITE(fh, n_realizations, 1, MPI_INTEGER, &
-               MPI_STATUS_IGNORE, mpi_err)
        end if
 
        ! Write all tally scores
@@ -686,6 +693,10 @@ contains
        end if
 
        TALLY_METADATA: do i = 1, n_tallies
+          ! Read number of realizations for global tallies
+          call MPI_FILE_READ(fh, tallies(i) % n_realizations, 1, &
+               MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err)
+
           ! Read dimensions of tally filters and scores and make sure they
           ! match
           call MPI_FILE_READ(fh, temp, 2, MPI_INTEGER, &
@@ -746,6 +757,10 @@ contains
           deallocate(int_array)
        end do TALLY_METADATA
 
+       ! Read number of realizations for global tallies
+       call MPI_FILE_READ(fh, n_realizations, 1, MPI_INTEGER, &
+            MPI_STATUS_IGNORE, mpi_err)
+
        ! Read number of global tallies and make sure it matches
        call MPI_FILE_READ(fh, temp, 1, MPI_INTEGER, MPI_STATUS_IGNORE, mpi_err)
        if (temp(1) /= N_GLOBAL_TALLIES) then
@@ -765,10 +780,6 @@ contains
 
        ! Read sum and sum squared
        if (temp(1) == 1) then
-          ! Read number of realizations
-          call MPI_FILE_READ(fh, n_realizations, 1, MPI_INTEGER, &
-               MPI_STATUS_IGNORE, mpi_err)
-
           TALLY_SCORES: do i = 1, n_tallies
              n = size(tallies(i) % scores, 1) * size(tallies(i) % scores, 2)
              call MPI_FILE_READ(fh, tallies(i) % scores, n, MPI_TALLYSCORE, &
@@ -882,6 +893,9 @@ contains
        end if
 
        TALLY_METADATA: do i = 1, n_tallies
+          ! Read number of realizations
+          read(UNIT_STATE) tallies(i) % n_realizations
+
           ! Read dimensions of tally filters and scores and make sure they
           ! match
           read(UNIT_STATE) temp(1:2)
@@ -932,6 +946,9 @@ contains
           deallocate(int_array)
        end do TALLY_METADATA
 
+       ! Read number of realizations for global tallies
+       read(UNIT_STATE) n_realizations
+
        ! Read number of global tallies and make sure it matches
        read(UNIT_STATE) temp(1)
        if (temp(1) /= N_GLOBAL_TALLIES) then
@@ -948,9 +965,6 @@ contains
        ! Read sum and sum squared
        read(UNIT_STATE) temp(1)
        if (temp(1) == 1) then
-          ! Read number of realizations
-          read(UNIT_STATE) n_realizations
-
           TALLY_SCORES: do i = 1, n_tallies
              do k = 1, size(tallies(i) % scores, 2)
                 do j = 1, size(tallies(i) % scores, 1)
diff --git a/src/utils/statepoint.py b/src/utils/statepoint.py
index a1c12a880..87ac92705 100644
--- a/src/utils/statepoint.py
+++ b/src/utils/statepoint.py
@@ -144,6 +144,9 @@ class StatePoint(BinaryFile):
             t = Tally()
             self.tallies.append(t)
 
+            # Read number of realizations
+            t.n_realizations = self._get_int()[0]
+
             # Read sizes of tallies
             t.n_score_bins = self._get_int()[0]
             t.n_filter_bins = self._get_int()[0]
@@ -186,6 +189,9 @@ class StatePoint(BinaryFile):
         if not self._metadata:
             self._read_metadata()
 
+        # Number of realizations for global tallies
+        self.n_realizations = self._get_int()[0]
+
         # Read global tallies
         n_global_tallies = self._get_int()[0]
         self.global_tallies = np.array(self._get_double(2*n_global_tallies))
@@ -196,7 +202,6 @@ class StatePoint(BinaryFile):
 
         # Read tally results
         if tallies_present:
-            self.n_realizations = self._get_int()[0]
 
             for t in self.tallies:
                 n = t.n_score_bins * t.n_filter_bins
diff --git a/src/utils/statepoint_histogram.py b/src/utils/statepoint_histogram.py
index e3dc012b5..4cde3b521 100755
--- a/src/utils/statepoint_histogram.py
+++ b/src/utils/statepoint_histogram.py
@@ -19,6 +19,9 @@ score = int(argv[3]) if len(argv) > 3 else 1
 # Create StatePoint object
 sp = StatePoint(filename)
 
+# Read number of realizations for global tallies
+sp.n_realizations = sp._get_int()[0]
+
 # Read global tallies
 n_global_tallies = sp._get_int()[0]
 sp.global_tallies = np.array(sp._get_double(2*n_global_tallies))
@@ -31,12 +34,12 @@ tallies_present = sp._get_int()[0]
 if not tallies_present:
     raise Exception("No tally data in state point!")
 
-# Calculate t-value for 95% two-sided CI
-n = sp._get_int()[0]
-t_value = scipy.stats.t.ppf(0.975, n - 1)
-
 # Loop over all tallies
 for i, t in enumerate(sp.tallies):
+    # Calculate t-value for 95% two-sided CI
+    n = t.n_realizations
+    t_value = scipy.stats.t.ppf(0.975, n - 1)
+
     # Create lists for tallies
     mean = []
     uncertainties = []
diff --git a/src/utils/statepoint_meshplot.py b/src/utils/statepoint_meshplot.py
index 1307c2b84..78e04a57e 100755
--- a/src/utils/statepoint_meshplot.py
+++ b/src/utils/statepoint_meshplot.py
@@ -19,6 +19,9 @@ score = int(argv[3]) if len(argv) > 3 else 1
 # Create StatePoint object
 sp = StatePoint(filename)
 
+# Read number of realizations for global tallies
+sp.n_realizations = sp._get_int()[0]
+
 # Read global tallies
 n_global_tallies = sp._get_int()[0]
 sp.global_tallies = np.array(sp._get_double(2*n_global_tallies))
@@ -31,13 +34,13 @@ tallies_present = sp._get_int()[0]
 if not tallies_present:
     raise Exception("No tally data in state point!")
 
-# Calculate t-value for 95% two-sided CI
-n = sp._get_int()[0]
-t_value = scipy.stats.t.ppf(0.975, n - 1)
-
 # Loop over all tallies
 print("Reading data...")
 for t in sp.tallies:
+    # Calculate t-value for 95% two-sided CI
+    n = t.n_realizations
+    t_value = scipy.stats.t.ppf(0.975, n - 1)
+
     n_bins = t.n_score_bins * t.n_filter_bins
     i_mesh = [f.type for f in t.filters].index('mesh')
 
diff --git a/src/utils/tally_conv.py b/src/utils/tally_conv.py
index 8079a6c7a..9741a5db3 100755
--- a/src/utils/tally_conv.py
+++ b/src/utils/tally_conv.py
@@ -77,6 +77,9 @@ for i_batch in range(len(files)):
     
     # Create StatePoint object
     sp = StatePoint(batch_filename)
+
+    # Read number of realizations for global tallies
+    sp.n_realizations = sp._get_int()[0]
     
     # Read global tallies
     n_global_tallies = sp._get_int()[0]
@@ -95,15 +98,15 @@ for i_batch in range(len(files)):
     uncert[i_batch] = [None for x in range(len(sp.tallies))]
     scoreType[i_batch] = [None for x in range(len(sp.tallies))]     
     
-    # Calculate t-value for 95% two-sided CI
-    n = sp._get_int()[0]
-    t_value = scipy.stats.t.ppf(0.975, n - 1)
-    
-    # Store the batch count    
-    active_batches[i_batch] = n
-    
     # Loop over all tallies
     for i_tally, t in enumerate(sp.tallies):
+        # Calculate t-value for 95% two-sided CI
+        n = t.n_realizations
+        t_value = scipy.stats.t.ppf(0.975, n - 1)
+    
+        # Store the batch count    
+        active_batches[i_batch] = n
+    
         # Resize the 2nd dimension      
         mean[i_batch][i_tally] = [None for x in range(t.n_filter_bins)]
         uncert[i_batch][i_tally] = [None for x in range(t.n_filter_bins)]