diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index b37a5880c..c837bf085 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -1059,6 +1059,7 @@ class ScatterMatrixXS(MultiGroupXS):
             scatter_p1 = scatter_p1.get_slice(scores=['scatter-P1'])
             energy_filter = openmc.Filter(type='energy')
             energy_filter.bins = self.energy_groups.group_edges
+            energy_filter.num_bins = self.num_groups
             scatter_p1 = scatter_p1.diagonalize_filter(energy_filter)
             rxn_tally = self.tallies['scatter'] - scatter_p1
         else:
@@ -1281,6 +1282,9 @@ class Chi(MultiGroupXS):
         sum_nu_fission_in = nu_fission_in.summation(filters=['energy'],
                                                     filter_bins=energy_bins)
 
+        # FIXME: Need ability to override energy groups with group numbers
+        # FIXME: Reverse from fast to thermal with energy groups
+
         # FIXME: CrossFilter for energy + energy messes up tally arithmetic
         sum_nu_fission_in.remove_filter(sum_nu_fission_in.filters[-1])
 
@@ -1291,7 +1295,12 @@ class Chi(MultiGroupXS):
                                        filter_bins=energy_bins)
 
         # FIXME: CrossFilter for energy + energy messes up tally arithmetic
-        norm.remove_filter(sum_nu_fission_in.filters[-1])
+        norm.remove_filter(norm.filters[-1])
+
+        energy_filter = openmc.Filter(type='energyout')
+        energy_filter.bins = self.energy_groups.group_edges
+        energy_filter.num_bins = self.num_groups
+        norm = norm.tile_filter(energy_filter)
 
         self._xs_tally /= norm
         self._xs_tally._mean = np.nan_to_num(self._xs_tally.mean)
diff --git a/openmc/tallies.py b/openmc/tallies.py
index c0fb1615f..28a89af91 100644
--- a/openmc/tallies.py
+++ b/openmc/tallies.py
@@ -2387,6 +2387,74 @@ class Tally(object):
 
         return tally_sum
 
+    def tile_filter(self, new_filter):
+        """Combines filters, scores and nuclides with another tally.
+
+        This is a helper method for the tally arithmetic methods. The filters,
+        scores and nuclides from both tallies are enumerated into all possible
+        combinations and expressed as CrossFilter, CrossScore and
+        CrossNuclide objects in the new derived tally.
+
+        Parameters
+        ----------
+        other : Tally
+            The tally on the right hand side of the outer product
+        binary_op : {'+', '-', '*', '/', '^'}
+            The binary operation in the outer product
+
+        Returns
+        -------
+        Tally
+            A new Tally outer that is the outer product with this one.
+
+        """
+
+        cv.check_type('new_filter', new_filter, Filter)
+
+        if new_filter in self.filters:
+            msg = 'Unable to tile Tally ID="{0}" which already ' \
+                  'contains a "{1}" filter'.format(self.id, new_filter.type)
+            raise ValueError(msg)
+
+        new_tally = copy.deepcopy(self)
+        new_tally.add_filter(new_filter)
+
+        num_filter_bins = new_tally.num_filter_bins
+        num_nuclides = new_tally.num_nuclides
+        num_score_bins = new_tally.num_score_bins
+        new_shape = (num_filter_bins, num_nuclides, num_score_bins)
+
+        repeat_indices = np.arange(0, new_tally.num_bins, new_filter.num_bins)
+        repeat_factor = new_filter.num_bins
+
+        if self.sum is not None:
+            new_tally._sum = np.zeros(new_shape, dtype=np.float64)
+        if self.sum_sq is not None:
+            new_tally._sum_sq = np.zeros(new_shape, dtype=np.float64)
+        if self.mean is not None:
+            new_tally._mean = np.zeros(new_shape, dtype=np.float64)
+        if self.std_dev is not None:
+            new_tally._std_dev = np.zeros(new_shape, dtype=np.float64)
+
+        for i in range(repeat_factor):
+            if self.sum is not None:
+                new_tally._sum[repeat_indices+i, :, :] = self.sum
+            if self.sum_sq is not None:
+                new_tally._sum_sq[repeat_indices+i, :, :] = self.sum_sq
+            if self.mean is not None:
+                new_tally._mean[repeat_indices+i, :, :] = self.mean
+            if self.std_dev is not None:
+                new_tally._std_dev[repeat_indices+i, :, :] = self.std_dev
+
+        # Correct each Filter's stride
+        stride = new_tally.num_nuclides * new_tally.num_score_bins
+        for filter in reversed(new_tally.filters):
+            filter.stride = stride
+            stride *= filter.num_bins
+
+        return new_tally
+
+
     def diagonalize_filter(self, new_filter):
         """Combines filters, scores and nuclides with another tally.
 
@@ -2424,7 +2492,14 @@ class Tally(object):
         num_score_bins = new_tally.num_score_bins
         new_shape = (num_filter_bins, num_nuclides, num_score_bins)
 
-        diag_indices = np.arange(0, new_tally.num_filter_bins, new_filter.num_bins+1)
+        indices = np.arange(0, new_filter.num_bins**2, new_filter.num_bins+1)
+        diag_indices = np.zeros(self.num_bins, dtype=np.int)
+        diag_factor = self.num_bins / new_filter.num_bins
+
+        for i in range(diag_factor):
+            start = i * new_filter.num_bins
+            end = (i+1) * new_filter.num_bins
+            diag_indices[start:end] = indices + (i * new_filter.num_bins**2)
 
         if self.sum is not None:
             new_tally._sum = np.zeros(new_shape, dtype=np.float64)