diff --git a/openmc/aggregate.py b/openmc/aggregate.py
new file mode 100644
index 0000000000..67995a41f2
--- /dev/null
+++ b/openmc/aggregate.py
@@ -0,0 +1,409 @@
+import sys
+from numbers import Integral
+
+import numpy as np
+
+from openmc import Filter, Nuclide
+from openmc.cross import CrossScore, CrossNuclide, CrossFilter
+from openmc.filter import _FILTER_TYPES
+import openmc.checkvalue as cv
+
+
+if sys.version_info[0] >= 3:
+    basestring = str
+
+# Acceptable tally aggregation operations
+_TALLY_AGGREGATE_OPS = ['sum', 'mean']
+
+
+class AggregateScore(object):
+    """A special-purpose tally score used to encapsulate an aggregate of a
+    subset or all of tally's scores for tally aggregation.
+
+    Parameters
+    ----------
+    scores : Iterable of str or CrossScore
+        The scores included in the aggregation
+    aggregate_op : str
+        The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
+        to aggregate across a tally's scores with this AggregateScore
+
+    Attributes
+    ----------
+    scores : Iterable of str or CrossScore
+        The scores included in the aggregation
+    aggregate_op : str
+        The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
+        to aggregate across a tally's scores with this AggregateScore
+
+    """
+
+    def __init__(self, scores=None, aggregate_op=None):
+
+        self._scores = None
+        self._aggregate_op = None
+
+        if scores is not None:
+            self.scores = scores
+        if aggregate_op is not None:
+            self.aggregate_op = aggregate_op
+
+    def __hash__(self):
+        return hash(repr(self))
+
+    def __eq__(self, other):
+        return str(other) == str(self)
+
+    def __ne__(self, other):
+        return not self == other
+
+    def __deepcopy__(self, memo):
+        existing = memo.get(id(self))
+
+        # If this is the first time we have tried to copy this object, create a copy
+        if existing is None:
+            clone = type(self).__new__(type(self))
+            clone._scores = self.scores
+            clone._aggregate_op = self.aggregate_op
+
+            memo[id(self)] = clone
+
+            return clone
+
+        # If this object has been copied before, return the first copy made
+        else:
+            return existing
+
+    def __repr__(self):
+        string = ', '.join(map(str, self.scores))
+        string = '{0}({1})'.format(self.aggregate_op, string)
+        return string
+
+    @property
+    def scores(self):
+        return self._scores
+
+    @property
+    def aggregate_op(self):
+        return self._aggregate_op
+
+    @scores.setter
+    def scores(self, scores):
+        cv.check_iterable_type('scores', scores, basestring)
+        self._scores = scores
+
+    @aggregate_op.setter
+    def aggregate_op(self, aggregate_op):
+        cv.check_type('aggregate_op', aggregate_op, (basestring, CrossScore))
+        cv.check_value('aggregate_op', aggregate_op, _TALLY_AGGREGATE_OPS)
+        self._aggregate_op = aggregate_op
+
+
+class AggregateNuclide(object):
+    """A special-purpose tally nuclide used to encapsulate an aggregate of a
+    subset or all of tally's nuclides for tally aggregation.
+
+    Parameters
+    ----------
+    nuclides : Iterable of str or Nuclide or CrossNuclide
+        The nuclides included in the aggregation
+    aggregate_op : str
+        The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
+        to aggregate across a tally's nuclides with this AggregateNuclide
+
+    Attributes
+    ----------
+    nuclides : Iterable of str or Nuclide or CrossNuclide
+        The nuclides included in the aggregation
+    aggregate_op : str
+        The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
+        to aggregate across a tally's nuclides with this AggregateNuclide
+
+    """
+
+    def __init__(self, nuclides=None, aggregate_op=None):
+
+        self._nuclides = None
+        self._aggregate_op = None
+
+        if nuclides is not None:
+            self.nuclides = nuclides
+        if aggregate_op is not None:
+            self.aggregate_op = aggregate_op
+
+    def __hash__(self):
+        return hash(repr(self))
+
+    def __eq__(self, other):
+        return str(other) == str(self)
+
+    def __ne__(self, other):
+        return not self == other
+
+    def __deepcopy__(self, memo):
+        existing = memo.get(id(self))
+
+        # If this is the first time we have tried to copy this object, create a copy
+        if existing is None:
+            clone = type(self).__new__(type(self))
+            clone._nuclides = self.nuclides
+            clone._aggregate_op = self._aggregate_op
+
+            memo[id(self)] = clone
+
+            return clone
+
+        # If this object has been copied before, return the first copy made
+        else:
+            return existing
+
+    def __repr__(self):
+
+        # Append each nuclide in the aggregate to the string
+        string = '{0}('.format(self.aggregate_op)
+        names = [nuclide.name if isinstance(nuclide, Nuclide) else str(nuclide)
+                 for nuclide in self.nuclides]
+        string += ', '.join(map(str, names)) + ')'
+        return string
+
+    @property
+    def nuclides(self):
+        return self._nuclides
+
+    @property
+    def aggregate_op(self):
+        return self._aggregate_op
+
+    @nuclides.setter
+    def nuclides(self, nuclides):
+        cv.check_iterable_type('nuclides', nuclides,
+                               (basestring, Nuclide, CrossNuclide))
+        self._nuclides = nuclides
+
+    @aggregate_op.setter
+    def aggregate_op(self, aggregate_op):
+        cv.check_type('aggregate_op', aggregate_op, basestring)
+        cv.check_value('aggregate_op', aggregate_op, _TALLY_AGGREGATE_OPS)
+        self._aggregate_op = aggregate_op
+
+
+class AggregateFilter(object):
+    """A special-purpose tally filter used to encapsulate an aggregate of a
+    subset or all of a tally filter's bins for tally aggregation.
+
+    Parameters
+    ----------
+    aggregate_filter : Filter or CrossFilter
+        The filter included in the aggregation
+    bins : Iterable of tuple
+        The filter bins included in the aggregation
+    aggregate_op : str
+        The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
+        to aggregate across a tally filter's bins with this AggregateFilter
+
+    Attributes
+    ----------
+    type : str
+        The type of the aggregatefilter (e.g., 'sum(energy)', 'sum(cell)')
+    aggregate_filter : filter
+        The filter included in the aggregation
+    aggregate_op : str
+        The tally aggregation operator (e.g., 'sum', 'mean', etc.) used
+        to aggregate across a tally filter's bins with this AggregateFilter
+    bins : Iterable of tuple
+        The filter bins included in the aggregation
+    num_bins : Integral
+        The number of filter bins (always 1 if aggregate_filter is defined)
+    stride : Integral
+        The number of filter, nuclide and score bins within each of this
+        aggregatefilter's bins.
+
+    """
+
+    def __init__(self, aggregate_filter=None, bins=None, aggregate_op=None):
+
+        self._type = '{0}({1})'.format(aggregate_op, aggregate_filter.type)
+        self._bins = None
+        self._stride = None
+
+        self._aggregate_filter = None
+        self._aggregate_op = None
+
+        if aggregate_filter is not None:
+            self.aggregate_filter = aggregate_filter
+        if bins is not None:
+            self.bins = bins
+        if aggregate_op is not None:
+            self.aggregate_op = aggregate_op
+
+    def __hash__(self):
+        return hash((self.type, self.bins, self.aggregate_op))
+
+    def __eq__(self, other):
+        return str(other) == str(self)
+
+    def __ne__(self, other):
+        return not self == other
+
+    def __repr__(self):
+        string = 'AggregateFilter\n'
+        string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self.type)
+        string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', self.bins)
+        return string
+
+    def __deepcopy__(self, memo):
+        existing = memo.get(id(self))
+
+        # If this is the first time we have tried to copy this object, create a copy
+        if existing is None:
+            clone = type(self).__new__(type(self))
+            clone._type = self.type
+            clone._aggregate_filter = self.aggregate_filter
+            clone._aggregate_op = self.aggregate_op
+            clone._bins = self._bins
+            clone._stride = self.stride
+
+            memo[id(self)] = clone
+
+            return clone
+
+        # If this object has been copied before, return the first copy made
+        else:
+            return existing
+
+    @property
+    def aggregate_filter(self):
+        return self._aggregate_filter
+
+    @property
+    def aggregate_op(self):
+        return self._aggregate_op
+
+    @property
+    def type(self):
+        return self._type
+
+    @property
+    def bins(self):
+        return self._bins
+
+    @property
+    def num_bins(self):
+        return 1 if self.aggregate_filter else 0
+
+    @property
+    def stride(self):
+        return self._stride
+
+    @type.setter
+    def type(self, filter_type):
+        if filter_type not in _FILTER_TYPES.values():
+            msg = 'Unable to set AggregateFilter type to "{0}" since it ' \
+                  'is not one of the supported types'.format(filter_type)
+            raise ValueError(msg)
+
+        self._type = filter_type
+
+    @aggregate_filter.setter
+    def aggregate_filter(self, aggregate_filter):
+        cv.check_type('aggregate_filter', aggregate_filter, (Filter, CrossFilter))
+        self._aggregate_filter = aggregate_filter
+
+    @bins.setter
+    def bins(self, bins):
+        cv.check_iterable_type('bins', bins, (Integral, tuple))
+        self._bins = bins
+
+    @aggregate_op.setter
+    def aggregate_op(self, aggregate_op):
+        cv.check_type('aggregate_op', aggregate_op, basestring)
+        cv.check_value('aggregate_op', aggregate_op, _TALLY_AGGREGATE_OPS)
+        self._aggregate_op = aggregate_op
+
+    @stride.setter
+    def stride(self, stride):
+        self._stride = stride
+
+    def get_bin_index(self, filter_bin):
+        """Returns the index in the AggregateFilter for some bin.
+
+        Parameters
+        ----------
+        filter_bin : Integral or tuple of Real
+            A tuple of value(s) corresponding to the bin of interest in
+            the aggregated filter. The bin is the integer ID for 'material',
+            'surface', 'cell', 'cellborn', and 'universe' Filters. The bin
+            is the integer cell instance ID for 'distribcell' Filters. The
+            bin is a 2-tuple of floats for 'energy' and 'energyout' filters
+            corresponding to the energy boundaries of the bin of interest.
+            The bin is a (x,y,z) 3-tuple for 'mesh' filters corresponding to
+            the mesh cell of interest.
+
+        Returns
+        -------
+        filter_index : Integral
+             The index in the Tally data array for this filter bin. For an
+             AggregateTally the filter bin index is always unity.
+
+        Raises
+        ------
+        ValueError
+            When the filter_bin is not part of the aggregated filter's bins
+
+        """
+
+        if filter_bin not in self.bins:
+            msg = 'Unable to get the bin index for AggregateFilter since ' \
+                  '"{0}" is not one of the bins'.format(filter_bin)
+            raise ValueError(msg)
+        else:
+            return 0
+
+    def get_pandas_dataframe(self, datasize, summary=None):
+        """Builds a Pandas DataFrame for the AggregateFilter's bins.
+
+        This method constructs a Pandas DataFrame object for the AggregateFilter
+        with columns annotated by filter bin information. This is a helper
+        method for the Tally.get_pandas_dataframe(...) method.
+
+        Parameters
+        ----------
+        datasize : Integral
+            The total number of bins in the tally corresponding to this filter
+        summary : None or Summary
+            An optional Summary object to be used to construct columns for
+            distribcell tally filters (default is None). NOTE: This parameter
+            is not used by the AggregateFilter and simply mirrors the method
+            signature for the CrossFilter.
+
+        Returns
+        -------
+        pandas.DataFrame
+            A Pandas DataFrame with columns of strings that characterize the
+            aggregatefilter's bins. Each entry in the DataFrame will include
+            one or more aggregation operations used to construct the
+            aggregatefilter's bins. The number of rows in the DataFrame is the
+            same as the total number of bins in the corresponding tally, with
+            the filter bins appropriately tiled to map to the corresponding
+            tally bins.
+
+        See also
+        --------
+        Tally.get_pandas_dataframe(), Filter.get_pandas_dataframe(),
+        CrossFilter.get_pandas_dataframe()
+
+        """
+
+        import pandas as pd
+
+        # Construct a sring representing the filter aggregation
+        aggregate_bin = '{0}('.format(self.aggregate_op)
+        aggregate_bin += ', '.join(map(str, self.bins)) + ')'
+
+        # Construct NumPy array of bin repeated for each element in dataframe
+        aggregate_bin_array = np.array([aggregate_bin])
+        aggregate_bin_array = np.repeat(aggregate_bin_array, datasize)
+
+        # Construct Pandas DataFrame for the AggregateFilter
+        df = pd.DataFrame({self.aggregate_filter.type: aggregate_bin_array})
+        return df
\ No newline at end of file
diff --git a/openmc/cross.py b/openmc/cross.py
index c03ee51885..ee0fcb4c82 100644
--- a/openmc/cross.py
+++ b/openmc/cross.py
@@ -34,7 +34,7 @@ class CrossScore(object):
         The right score in the outer product
     binary_op : str
         The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
-        combine two tally's scores with this CrossNuclide
+        combine two tally's scores with this CrossScore
 
     """
 
@@ -245,6 +245,8 @@ class CrossFilter(object):
 
     Attributes
     ----------
+    type : str
+        The type of the crossfilter (e.g., 'energy / energy')
     left_filter : Filter or CrossFilter
         The left filter in the outer product
     right_filter : Filter or CrossFilter
@@ -252,6 +254,14 @@ class CrossFilter(object):
     binary_op : str
         The tally arithmetic binary operator (e.g., '+', '-', etc.) used to
         combine two tally's filter bins with this CrossFilter
+    bins : dict of Iterable
+        A dictionary of the bins from each filter keyed by the types of the
+        left / right filters
+    num_bins : Integral
+        The number of filter bins (always 1 if aggregate_filter is defined)
+    stride : Integral
+        The number of filter, nuclide and score bins within each of this
+        crossfilter's bins.
 
     """
 
@@ -310,7 +320,6 @@ class CrossFilter(object):
             clone._binary_op = self.binary_op
             clone._type = self.type
             clone._bins = self._bins
-            clone._num_bins = self.num_bins
             clone._stride = self.stride
 
             memo[id(self)] = clone
@@ -355,8 +364,8 @@ class CrossFilter(object):
     @type.setter
     def type(self, filter_type):
         if filter_type not in _FILTER_TYPES.values():
-            msg = 'Unable to set Filter type to "{0}" since it is not one ' \
-                  'of the supported types'.format(filter_type)
+            msg = 'Unable to set CrossFilter type to "{0}" since it ' \
+                  'is not one of the supported types'.format(filter_type)
             raise ValueError(msg)
 
         self._type = filter_type
diff --git a/openmc/filter.py b/openmc/filter.py
index d61fb648f6..a0c777fa07 100644
--- a/openmc/filter.py
+++ b/openmc/filter.py
@@ -35,8 +35,8 @@ class Filter(object):
     Attributes
     ----------
     type : str
-        The type of the tally filter.
-    bins : Integral or Iterable of Integral or Iterable of Real
+        The type of the tally filter
+    bins : Integral or Iterable of Real
         The bins for the filter
     num_bins : Integral
         The number of filter bins
@@ -227,12 +227,12 @@ class Filter(object):
 
         self._stride = stride
 
-    def can_merge(self, filter):
+    def can_merge(self, other):
         """Determine if filter can be merged with another.
 
         Parameters
         ----------
-        filter : Filter
+        other : Filter
             Filter to compare with
 
         Returns
@@ -242,11 +242,11 @@ class Filter(object):
 
         """
 
-        if not isinstance(filter, Filter):
+        if not isinstance(other, Filter):
             return False
 
         # Filters must be of the same type
-        elif self.type != filter.type:
+        elif self.type != other.type:
             return False
 
         # Distribcell filters cannot have more than one bin
@@ -264,12 +264,12 @@ class Filter(object):
         else:
             return True
 
-    def merge(self, filter):
+    def merge(self, other):
         """Merge this filter with another.
 
         Parameters
         ----------
-        filter : Filter
+        other : Filter
             Filter to merge with
 
         Returns
@@ -279,16 +279,16 @@ class Filter(object):
 
         """
 
-        if not self.can_merge(filter):
+        if not self.can_merge(other):
             msg = 'Unable to merge "{0}" with "{1}" ' \
-                  'filters'.format(self.type, filter.type)
+                  'filters'.format(self.type, other.type)
             raise ValueError(msg)
 
         # Create deep copy of filter to return as merged filter
         merged_filter = copy.deepcopy(self)
 
         # Merge unique filter bins
-        merged_bins = list(set(np.concatenate((self.bins, filter.bins))))
+        merged_bins = list(set(np.concatenate((self.bins, other.bins))))
         merged_filter.bins = merged_bins
         merged_filter.num_bins = len(merged_bins)
 
diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py
index 47154357ef..e20a6a320e 100644
--- a/openmc/mgxs/mgxs.py
+++ b/openmc/mgxs/mgxs.py
@@ -521,8 +521,8 @@ class MGXS(object):
                 self.tallies[key].add_trigger(trigger_clone)
 
             # Add all non-domain specific Filters (e.g., 'energy') to the Tally
-            for filter in filters:
-                self.tallies[key].add_filter(filter)
+            for add_filter in filters:
+                self.tallies[key].add_filter(add_filter)
 
             # If this is a by-nuclide cross-section, add all nuclides to Tally
             if self.by_nuclide and score != 'flux':
@@ -787,15 +787,15 @@ class MGXS(object):
             std_dev = tally.get_reshaped_data(value='std_dev')
 
             # Sum across all applicable fine energy group filters
-            for i, filter in enumerate(tally.filters):
-                if 'energy' not in filter.type:
+            for i, tally_filter in enumerate(tally.filters):
+                if 'energy' not in tally_filter.type:
                     continue
-                elif len(filter.bins) != len(fine_edges):
+                elif len(tally_filter.bins) != len(fine_edges):
                     continue
-                elif not np.allclose(filter.bins, fine_edges):
+                elif not np.allclose(tally_filter.bins, fine_edges):
                     continue
                 else:
-                    filter.bins = coarse_groups.group_edges
+                    tally_filter.bins = coarse_groups.group_edges
                     mean = np.add.reduceat(mean, energy_indices, axis=i)
                     std_dev = np.add.reduceat(std_dev**2, energy_indices, axis=i)
                     std_dev = np.sqrt(std_dev)
diff --git a/openmc/statepoint.py b/openmc/statepoint.py
index 3c0759f124..49ac43590e 100644
--- a/openmc/statepoint.py
+++ b/openmc/statepoint.py
@@ -3,6 +3,7 @@ import re
 import numpy as np
 
 import openmc
+import openmc.checkvalue as cv
 
 if sys.version > '3':
     long = int
@@ -377,18 +378,18 @@ class StatePoint(object):
                     bins = self._f['{0}{1}/bins'.format(subbase, j)].value
 
                     # Create Filter object
-                    filter = openmc.Filter(filter_type, bins)
-                    filter.num_bins = n_bins
+                    new_filter = openmc.Filter(filter_type, bins)
+                    new_filter.num_bins = n_bins
 
                     if filter_type == 'mesh':
                         mesh_ids = self._f['tallies/meshes/ids'].value
                         mesh_keys = self._f['tallies/meshes/keys'].value
 
                         key = mesh_keys[mesh_ids == bins][0]
-                        filter.mesh = self.meshes[key]
+                        new_filter.mesh = self.meshes[key]
 
                     # Add Filter to the Tally
-                    tally.add_filter(filter)
+                    tally.add_filter(new_filter)
 
                 # Read Nuclide bins
                 nuclide_names = \
@@ -406,11 +407,11 @@ class StatePoint(object):
 
                 # Compute and set the filter strides
                 for i in range(n_filters):
-                    filter = tally.filters[i]
-                    filter.stride = n_score_bins * len(nuclide_names)
+                    tally_filter = tally.filters[i]
+                    tally_filter.stride = n_score_bins * len(nuclide_names)
 
                     for j in range(i+1, n_filters):
-                        filter.stride *= tally.filters[j].num_bins
+                        tally_filter.stride *= tally.filters[j].num_bins
 
                 # Read scattering moment order strings (e.g., P3, Y1,2, etc.)
                 moments = self._f['{0}{1}/moment_orders'.format(
@@ -544,13 +545,13 @@ class StatePoint(object):
                 contains_filters = True
 
                 # Iterate over the Filters requested by the user
-                for filter in filters:
+                for outer_filter in filters:
                     contains_filters = False
 
                     # Test if requested filter is a subset of any of the test
                     # tally's filters and if so continue to next filter
-                    for test_filter in test_tally.filters:
-                        if test_filter.is_subset(filter):
+                    for inner_filter in test_tally.filters:
+                        if inner_filter.is_subset(outer_filter):
                             contains_filters = True
                             break
 
@@ -616,29 +617,29 @@ class StatePoint(object):
             tally.name = summary.tallies[tally_id].name
             tally.with_summary = True
 
-            for filter in tally.filters:
-                if filter.type == 'surface':
+            for tally_filter in tally.filters:
+                if tally_filter.type == 'surface':
                     surface_ids = []
-                    for bin in filter.bins:
+                    for bin in tally_filter.bins:
                         surface_ids.append(summary.surfaces[bin].id)
-                    filter.bins = surface_ids
+                    tally_filter.bins = surface_ids
 
-                if filter.type in ['cell', 'distribcell']:
+                if tally_filter.type in ['cell', 'distribcell']:
                     distribcell_ids = []
-                    for bin in filter.bins:
+                    for bin in tally_filter.bins:
                         distribcell_ids.append(summary.cells[bin].id)
-                    filter.bins = distribcell_ids
+                    tally_filter.bins = distribcell_ids
 
-                if filter.type == 'universe':
+                if tally_filter.type == 'universe':
                     universe_ids = []
-                    for bin in filter.bins:
+                    for bin in tally_filter.bins:
                         universe_ids.append(summary.universes[bin].id)
-                    filter.bins = universe_ids
+                    tally_filter.bins = universe_ids
 
-                if filter.type == 'material':
+                if tally_filter.type == 'material':
                     material_ids = []
-                    for bin in filter.bins:
+                    for bin in tally_filter.bins:
                         material_ids.append(summary.materials[bin].id)
-                    filter.bins = material_ids
+                    tally_filter.bins = material_ids
 
         self._summary = summary
diff --git a/openmc/summary.py b/openmc/summary.py
index ad1dd67cae..eb14d3bb81 100644
--- a/openmc/summary.py
+++ b/openmc/summary.py
@@ -556,11 +556,11 @@ class Summary(object):
                 bins = self._f['{0}/bins'.format(subsubbase)][...]
 
                 # Create Filter object
-                filter = openmc.Filter(filter_type, bins)
-                filter.num_bins = num_bins
+                new_filter = openmc.Filter(filter_type, bins)
+                new_filter.num_bins = num_bins
 
                 # Add Filter to the Tally
-                tally.add_filter(filter)
+                tally.add_filter(new_filter)
 
             # Add Tally to the global dictionary of all Tallies
             self.tallies[tally_id] = tally
diff --git a/openmc/tallies.py b/openmc/tallies.py
index 7ec1ca50cb..980c5fa2b5 100644
--- a/openmc/tallies.py
+++ b/openmc/tallies.py
@@ -13,6 +13,7 @@ import numpy as np
 
 from openmc import Mesh, Filter, Trigger, Nuclide
 from openmc.cross import CrossScore, CrossNuclide, CrossFilter
+from openmc.aggregate import AggregateScore, AggregateNuclide, AggregateFilter
 from openmc.filter import _FILTER_TYPES
 import openmc.checkvalue as cv
 from openmc.clean_xml import *
@@ -143,8 +144,8 @@ class Tally(object):
             clone._results_read = self._results_read
 
             clone._filters = []
-            for filter in self.filters:
-                clone.add_filter(copy.deepcopy(filter, memo))
+            for self_filter in self.filters:
+                clone.add_filter(copy.deepcopy(self_filter, memo))
 
             clone._nuclides = []
             for nuclide in self.nuclides:
@@ -174,8 +175,8 @@ class Tally(object):
         if len(self.filters) != len(other.filters):
             return False
 
-        for filter in self.filters:
-            if filter not in other.filters:
+        for self_filter in self.filters:
+            if self_filter not in other.filters:
                 return False
 
         # Check all nuclides
@@ -212,9 +213,9 @@ class Tally(object):
 
         string += '{0: <16}{1}\n'.format('\tFilters', '=\t')
 
-        for filter in self.filters:
-            string += '{0: <16}\t\t{1}\t{2}\n'.format('', filter.type,
-                                                          filter.bins)
+        for self_filter in self.filters:
+            string += '{0: <16}\t\t{1}\t{2}\n'.format('', self_filter.type,
+                                                          self_filter.bins)
 
         string += '{0: <16}{1}'.format('\tNuclides', '=\t')
 
@@ -259,12 +260,16 @@ class Tally(object):
     def num_scores(self):
         return len(self._scores)
 
+    @property
+    def num_filters(self):
+        return len(self.filters)
+
     @property
     def num_filter_bins(self):
         num_bins = 1
 
-        for filter in self.filters:
-            num_bins *= filter.num_bins
+        for self_filter in self.filters:
+            num_bins *= self_filter.num_bins
 
         return num_bins
 
@@ -455,33 +460,63 @@ class Tally(object):
         else:
             self._name = ''
 
-    def add_filter(self, filter):
+    def add_filter(self, new_filter):
         """Add a filter to the tally
 
         Parameters
         ----------
-        filter : openmc.filter.Filter
-            Filter to add
+        new_filter : Filter, CrossFilter or AggregateFilter
+            A filter to specify a discretization of the tally across some
+            dimension (e.g., 'energy', 'cell'). The filter should be a Filter
+            object when a user is adding filters to a Tally for input file
+            generation or when the Tally is created from a StatePoint. The
+            filter may be a CrossFilter or AggregateFilter for derived tallies
+            created by tally arithmetic.
 
         """
 
-        if not isinstance(filter, (Filter, CrossFilter)):
+        if not isinstance(new_filter, (Filter, CrossFilter, AggregateFilter)):
             msg = 'Unable to add Filter "{0}" to Tally ID="{1}" since it is ' \
-                  'not a Filter object'.format(filter, self.id)
+                  'not a Filter object'.format(new_filter, self.id)
             raise ValueError(msg)
 
-        self._filters.append(filter)
+        # If the filter is already in the Tally, raise an error
+        if new_filter in self.filters:
+            msg = 'Unable to add a duplicate filter "{0}" to Tally ID="{1}" ' \
+                  'since duplicate filters are not supported in the OpenMC ' \
+                  'Python API'.format(new_filter, self.id)
+            raise ValueError(msg)
+
+        self._filters.append(new_filter)
 
     def add_nuclide(self, nuclide):
         """Specify that scores for a particular nuclide should be accumulated
 
         Parameters
         ----------
-        nuclide : openmc.nuclide.Nuclide
-            Nuclide to add
+        nuclide : str, Nuclide, CrossNuclide or AggregateNuclide
+            Nuclide to add to the tally. The nuclide should be a Nuclide object
+            when a user is adding nuclides to a Tally for input file generation.
+            The nuclide is a str when a Tally is created from a StatePoint file
+            (e.g., 'H-1', 'U-235') unless a Summary has been linked with the
+            StatePoint. The nuclide may be a CrossNuclide or AggregateNuclide
+            for derived tallies created by tally arithmetic.
 
         """
 
+        if not isinstance(nuclide, (basestring, Nuclide,
+                                    CrossNuclide, AggregateNuclide)):
+            msg = 'Unable to add nuclide "{0}" to Tally ID="{1}" since it is ' \
+                  'not a Nuclide object'.format(nuclide)
+            raise ValueError(msg)
+
+        # If the nuclide is already in the Tally, raise an error
+        if nuclide in self.nuclides:
+            msg = 'Unable to add a duplicate nuclide "{0}" to Tally ID="{1}" ' \
+                  'since duplicate nuclides are not supported in the OpenMC ' \
+                  'Python API'.format(nuclide, self.id)
+            raise ValueError(msg)
+
         self._nuclides.append(nuclide)
 
     def add_score(self, score):
@@ -489,19 +524,23 @@ class Tally(object):
 
         Parameters
         ----------
-        score : str
-            Score to be accumulated, e.g. 'flux'
+        score : str, CrossScore or AggregateScore
+            A score to be accumulated (e.g., 'flux', 'nu-fission'). The score
+            should be a str when a user is adding scores to a Tally for input
+            file generation or when the Tally is created from a StatePoint. The
+            score may be a CrossScore or AggregateScore for derived tallies
+            created by tally arithmetic.
 
         """
 
-        if not isinstance(score, (basestring, CrossScore)):
+        if not isinstance(score, (basestring, CrossScore, AggregateScore)):
             msg = 'Unable to add score "{0}" to Tally ID="{1}" since it is ' \
                   'not a string'.format(score, self.id)
             raise ValueError(msg)
 
         # If the score is already in the Tally, raise an error
         if score in self.scores:
-            msg = 'Unable to add a duplicate score {0} to Tally ID="{1}" ' \
+            msg = 'Unable to add a duplicate score "{0}" to Tally ID="{1}" ' \
                   'since duplicate scores are not supported in the OpenMC ' \
                   'Python API'.format(score, self.id)
             raise ValueError(msg)
@@ -509,7 +548,7 @@ class Tally(object):
         # Normal score strings
         if isinstance(score, basestring):
             self._scores.append(score.strip())
-        # CrossScores
+        # CrossScores and AggrgateScore
         else:
             self._scores.append(score)
 
@@ -601,22 +640,22 @@ class Tally(object):
 
         self._scores.remove(score)
 
-    def remove_filter(self, filter):
+    def remove_filter(self, old_filter):
         """Remove a filter from the tally
 
         Parameters
         ----------
-        filter : openmc.filter.Filter
+        old_filter : openmc.filter.Filter
             Filter to remove
 
         """
 
-        if filter not in self.filters:
+        if old_filter not in self.filters:
             msg = 'Unable to remove filter "{0}" from Tally ID="{1}" since the ' \
-                  'Tally does not contain this filter'.format(filter, self.id)
+                  'Tally does not contain this filter'.format(old_filter, self.id)
             ValueError(msg)
 
-        self._filters.remove(filter)
+        self._filters.remove(old_filter)
 
     def remove_nuclide(self, nuclide):
         """Remove a nuclide from the tally
@@ -760,13 +799,13 @@ class Tally(object):
             element.set("name", self.name)
 
         # Optional Tally filters
-        for filter in self.filters:
+        for self_filter in self.filters:
             subelement = ET.SubElement(element, "filter")
-            subelement.set("type", str(filter.type))
+            subelement.set("type", str(self_filter.type))
 
-            if filter.bins is not None:
+            if self_filter.bins is not None:
                 bins = ''
-                for bin in filter.bins:
+                for bin in self_filter.bins:
                     bins += '{0} '.format(bin)
 
                 subelement.set("bins", bins.rstrip(' '))
@@ -818,7 +857,7 @@ class Tally(object):
 
         Returns
         -------
-        filter : openmc.filter.Filter
+        filter_found : openmc.filter.Filter
             Filter from this tally with matching type, or None if no matching
             Filter is found
 
@@ -829,21 +868,21 @@ class Tally(object):
 
         """
 
-        filter = None
+        filter_found = None
 
         # Look through all of this Tally's Filters for the type requested
         for test_filter in self.filters:
             if test_filter.type == filter_type:
-                filter = test_filter
+                filter_found = test_filter
                 break
 
         # If we did not find the Filter, throw an Exception
-        if filter is None:
+        if filter_found is None:
             msg = 'Unable to find filter type "{0}" in ' \
                   'Tally ID="{1}"'.format(filter_type, self.id)
             raise ValueError(msg)
 
-        return filter
+        return filter_found
 
     def get_filter_index(self, filter_type, filter_bin):
         """Returns the index in the Tally's results array for a Filter bin
@@ -868,10 +907,10 @@ class Tally(object):
         """
 
         # Find the equivalent Filter in this Tally's list of Filters
-        filter = self.find_filter(filter_type)
+        filter_found = self.find_filter(filter_type)
 
         # Get the index for the requested bin from the Filter and return it
-        filter_index = filter.get_bin_index(filter_bin)
+        filter_index = filter_found.get_bin_index(filter_bin)
         return filter_index
 
     def get_nuclide_index(self, nuclide):
@@ -991,12 +1030,12 @@ class Tally(object):
             filter_indices = []
 
             # Loop over all of the Tally's Filters
-            for i, filter in enumerate(self.filters):
+            for i, self_filter in enumerate(self.filters):
                 user_filter = False
 
                 # If a user-requested Filter, get the user-requested bins
                 for j, test_filter in enumerate(filters):
-                    if filter.type == test_filter:
+                    if self_filter.type == test_filter:
                         bins = filter_bins[j]
                         user_filter = True
                         break
@@ -1004,36 +1043,36 @@ class Tally(object):
                 # If not a user-requested Filter, get all bins
                 if not user_filter:
                     # Create list of 2- or 3-tuples tuples for mesh cell bins
-                    if filter.type == 'mesh':
-                        dimension = filter.mesh.dimension
+                    if self_filter.type == 'mesh':
+                        dimension = self_filter.mesh.dimension
                         xyz = map(lambda x: np.arange(1, x+1), dimension)
                         bins = list(itertools.product(*xyz))
 
                     # Create list of 2-tuples for energy boundary bins
-                    elif filter.type in ['energy', 'energyout']:
+                    elif self_filter.type in ['energy', 'energyout']:
                         bins = []
-                        for k in range(filter.num_bins):
-                            bins.append((filter.bins[k], filter.bins[k+1]))
+                        for k in range(self_filter.num_bins):
+                            bins.append((self_filter.bins[k], self_filter.bins[k+1]))
 
                     # Create list of cell instance IDs for distribcell Filters
-                    elif filter.type == 'distribcell':
-                        bins = np.arange(filter.num_bins)
+                    elif self_filter.type == 'distribcell':
+                        bins = np.arange(self_filter.num_bins)
 
                     # Create list of IDs for bins for all other filter types
                     else:
-                        bins = filter.bins
+                        bins = self_filter.bins
 
                 # Initialize a NumPy array for the Filter bin indices
                 filter_indices.append(np.zeros(len(bins), dtype=np.int))
 
                 # Add indices for each bin in this Filter to the list
                 for j, bin in enumerate(bins):
-                    filter_index = self.get_filter_index(filter.type, bin)
+                    filter_index = self.get_filter_index(self_filter.type, bin)
                     filter_indices[i][j] = filter_index
 
                 # Account for stride in each of the previous filters
                 for indices in filter_indices[:i]:
-                    indices *= filter.num_bins
+                    indices *= self_filter.num_bins
 
             # Apply outer product sum between all filter bin indices
             filter_indices = list(map(sum, itertools.product(*filter_indices)))
@@ -1275,8 +1314,8 @@ class Tally(object):
         if filters:
 
             # Append each Filter's DataFrame to the overall DataFrame
-            for filter in self.filters:
-                filter_df = filter.get_pandas_dataframe(data_size, summary)
+            for self_filter in self.filters:
+                filter_df = self_filter.get_pandas_dataframe(data_size, summary)
                 df = pd.concat([df, filter_df], axis=1)
 
         # Include DataFrame column for nuclides if user requested it
@@ -1365,8 +1404,8 @@ class Tally(object):
 
         # Build a new array shape with one dimension per filter
         new_shape = ()
-        for filter in self.filters:
-            new_shape += (filter.num_bins, )
+        for self_filter in self.filters:
+            new_shape += (self_filter.num_bins, )
         new_shape += (self.num_nuclides,)
         new_shape += (self.num_scores,)
 
@@ -1459,8 +1498,9 @@ class Tally(object):
             # Create an HDF5 sub-group for the Filters
             filter_group = tally_group.create_group('filters')
 
-            for filter in self.filters:
-                filter_group.create_dataset(filter.type, data=filter.bins)
+            for self_filter in self.filters:
+                filter_group.create_dataset(self_filter.type,
+                                            filter=self_filter.bins)
 
             # Add all results to the main HDF5 group for the Tally
             tally_group.create_dataset('sum', data=self.sum)
@@ -1503,8 +1543,8 @@ class Tally(object):
             tally_group['filters'] = {}
             filter_group = tally_group['filters']
 
-            for filter in self.filters:
-                filter_group[filter.type] = filter.bins
+            for self_filter in self.filters:
+                filter_group[self_filter.type] = self_filter.bins
 
             # Add all results to the main sub-dictionary for the Tally
             tally_group['sum'] = self.sum
@@ -1599,8 +1639,12 @@ class Tally(object):
             raise ValueError(msg)
 
         new_tally = Tally()
-        new_tally.with_batch_statistics = True
         new_tally._derived = True
+        new_tally.with_batch_statistics = True
+        new_tally._num_realizations = self.num_realizations
+        new_tally._estimator = self.estimator
+        new_tally._with_summary = self.with_summary
+        new_tally._sp_filename = self._sp_filename
 
         # Construct a combined derived name from the two tally operands
         if self.name != '' and other.name != '':
@@ -1698,14 +1742,21 @@ class Tally(object):
                 new_score = CrossScore(self_score, other_score, binary_op)
                 new_tally.add_score(new_score)
 
-        # Correct each Filter's stride
-        stride = new_tally.num_nuclides * new_tally.num_scores
-        for filter in reversed(new_tally.filters):
-            filter.stride = stride
-            stride *= filter.num_bins
+        # Update the new tally's filter strides
+        new_tally._update_filter_strides()
 
         return new_tally
 
+    def _update_filter_strides(self):
+        """Update each filter's stride based on the tally's nuclides and scores
+        for derived tallies created by tally arithmetic.
+        """
+
+        stride = self.num_nuclides * self.num_scores
+        for self_filter in reversed(self.filters):
+            self_filter.stride = stride
+            stride *= self_filter.num_bins
+
     def _align_tally_data(self, other, filter_product, nuclide_product,
                           score_product):
         """Aligns data from two tallies for tally arithmetic.
@@ -1748,26 +1799,26 @@ class Tally(object):
             set(other.filters).difference(set(self.filters))
 
         # Add filters present in self but not in other to other
-        for filter in other_missing_filters:
-            filter = copy.deepcopy(filter)
-            other._mean = np.repeat(other.mean, filter.num_bins, axis=0)
-            other._std_dev = np.repeat(other.std_dev, filter.num_bins, axis=0)
-            other.add_filter(filter)
+        for other_filter in other_missing_filters:
+            filter_copy = copy.deepcopy(other_filter)
+            other._mean = np.repeat(other.mean, filter_copy.num_bins, axis=0)
+            other._std_dev = np.repeat(other.std_dev, filter_copy.num_bins, axis=0)
+            other.add_filter(filter_copy)
 
         # Add filters present in other but not in self to self
-        for filter in self_missing_filters:
-            filter = copy.deepcopy(filter)
-            self._mean = np.repeat(self.mean, filter.num_bins, axis=0)
-            self._std_dev = np.repeat(self.std_dev, filter.num_bins, axis=0)
-            self.add_filter(filter)
+        for self_filter in self_missing_filters:
+            filter_copy = copy.deepcopy(self_filter)
+            self._mean = np.repeat(self.mean, filter_copy.num_bins, axis=0)
+            self._std_dev = np.repeat(self.std_dev, filter_copy.num_bins, axis=0)
+            self.add_filter(filter_copy)
 
         # Align other filters with self filters
-        for i, filter in enumerate(self.filters):
-            other_index = other.filters.index(filter)
+        for i, self_filter in enumerate(self.filters):
+            other_index = other.filters.index(self_filter)
 
             # If necessary, swap other filter
             if other_index != i:
-                other._swap_filters(filter, other.filters[i])
+                other._swap_filters(self_filter, other.filters[i])
 
         # Repeat and tile the data by nuclide in preparation for performing
         # the tensor product across nuclides.
@@ -1855,17 +1906,9 @@ class Tally(object):
                 if other_index != i:
                     other._swap_scores(score, other.scores[i])
 
-        # Correct the stride for other filters
-        stride = other.num_nuclides * other.num_scores
-        for filter in reversed(other.filters):
-            filter.stride = stride
-            stride *= filter.num_bins
-
-        # Correct the stride for self filters
-        stride = self.num_nuclides * self.num_scores
-        for filter in reversed(self.filters):
-            filter.stride = stride
-            stride *= filter.num_bins
+        # Update the tallies' filter strides
+        other._update_filter_strides()
+        self._update_filter_strides()
 
         data = {}
         data['self'] = {}
@@ -1927,16 +1970,13 @@ class Tally(object):
         self.filters[filter1_index] = filter2
         self.filters[filter2_index] = filter1
 
-        # Update the strides for each of the filters
-        stride = self.num_nuclides * self.num_scores
-        for filter in reversed(self.filters):
-            filter.stride = stride
-            stride *= filter.num_bins
+        # Update the tally's filter strides
+        self._update_filter_strides()
 
         # Construct lists of tuples for the bins in each of the two filters
         filters = [filter1.type, filter2.type]
         if filter1.type == 'distribcell':
-            filter1_bins = np.arange(filter.num_bins)
+            filter1_bins = np.arange(filter1.num_bins)
         else:
             filter1_bins = [(filter1.get_bin(i)) for i in range(filter1.num_bins)]
 
@@ -2161,8 +2201,8 @@ class Tally(object):
             new_tally.with_summary = self.with_summary
             new_tally.num_realization = self.num_realizations
 
-            for filter in self.filters:
-                new_tally.add_filter(filter)
+            for self_filter in self.filters:
+                new_tally.add_filter(self_filter)
             for nuclide in self.nuclides:
                 new_tally.add_nuclide(nuclide)
             for score in self.scores:
@@ -2235,8 +2275,8 @@ class Tally(object):
             new_tally.with_summary = self.with_summary
             new_tally.num_realization = self.num_realizations
 
-            for filter in self.filters:
-                new_tally.add_filter(filter)
+            for self_filter in self.filters:
+                new_tally.add_filter(self_filter)
             for nuclide in self.nuclides:
                 new_tally.add_nuclide(nuclide)
             for score in self.scores:
@@ -2310,8 +2350,8 @@ class Tally(object):
             new_tally.with_summary = self.with_summary
             new_tally.num_realization = self.num_realizations
 
-            for filter in self.filters:
-                new_tally.add_filter(filter)
+            for self_filter in self.filters:
+                new_tally.add_filter(self_filter)
             for nuclide in self.nuclides:
                 new_tally.add_nuclide(nuclide)
             for score in self.scores:
@@ -2385,8 +2425,8 @@ class Tally(object):
             new_tally.with_summary = self.with_summary
             new_tally.num_realization = self.num_realizations
 
-            for filter in self.filters:
-                new_tally.add_filter(filter)
+            for self_filter in self.filters:
+                new_tally.add_filter(self_filter)
             for nuclide in self.nuclides:
                 new_tally.add_nuclide(nuclide)
             for score in self.scores:
@@ -2449,8 +2489,8 @@ class Tally(object):
         if isinstance(power, Tally):
             new_tally = self.hybrid_product(power, binary_op='^')
 
-            # If original tally operands were sparse, sparsify the new tally
-            if self.sparse and other.sparse:
+            # If original tally operand was sparse, sparsify the new tally
+            if self.sparse:
                 new_tally.sparse = True
 
         elif isinstance(power, Real):
@@ -2464,8 +2504,8 @@ class Tally(object):
             new_tally.with_summary = self.with_summary
             new_tally.num_realization = self.num_realizations
 
-            for filter in self.filters:
-                new_tally.add_filter(filter)
+            for self_filter in self.filters:
+                new_tally.add_filter(self_filter)
             for nuclide in self.nuclides:
                 new_tally.add_nuclide(nuclide)
             for score in self.scores:
@@ -2605,7 +2645,7 @@ class Tally(object):
             parameter (e.g., [(1,), (0., 0.625e-6)]; default is []). Each bin
             in the list is the integer ID for 'material', 'surface', 'cell',
             'cellborn', and 'universe' Filters. Each bin is an integer for the
-            cell instance ID for 'distribcell Filters. Each bin is a 2-tuple of
+            cell instance ID for 'distribcell' Filters. Each bin is a 2-tuple of
             floats for 'energy' and 'energyout' filters corresponding to the
             energy boundaries of the bin of interest.  The bin is a (x,y,z)
             3-tuple for 'mesh' filters corresponding to the mesh cell of
@@ -2687,29 +2727,29 @@ class Tally(object):
 
             # Determine the filter indices from any of the requested filters
             for i, filter_type in enumerate(filters):
-                filter = new_tally.find_filter(filter_type)
+                find_filter = new_tally.find_filter(filter_type)
 
                 # Remove and/or reorder filter bins to user specifications
                 bin_indices = []
+                num_bins = 0
 
                 for filter_bin in filter_bins[i]:
-                    bin_index = filter.get_bin_index(filter_bin)
+                    bin_index = find_filter.get_bin_index(filter_bin)
                     if filter_type in ['energy', 'energyout']:
                         bin_indices.extend([bin_index, bin_index+1])
+                        num_bins += 1
                     elif filter_type == 'distribcell':
-                        indices = [(bin,) for bin in range(filter.num_bins)]
-                        bin_indices.extend(indices)
+                        bin_indices = [0]
+                        num_bins = find_filter.num_bins
                     else:
                         bin_indices.append(bin_index)
+                        num_bins += 1
 
-                filter.bins = filter.bins[bin_indices]
-                filter.num_bins = len(filter_bins[i])
+                find_filter.bins = find_filter.bins[bin_indices]
+                find_filter.num_bins = num_bins
 
-        # Correct each Filter's stride
-        stride = new_tally.num_nuclides * new_tally.num_scores
-        for filter in reversed(new_tally.filters):
-            filter.stride = stride
-            stride *= filter.num_bins
+        # Update the new tally's filter strides
+        new_tally._update_filter_strides()
 
         # If original tally was sparse, sparsify the sliced tally
         new_tally.sparse = self.sparse
@@ -2737,7 +2777,7 @@ class Tally(object):
             A list of the filter bins corresponding to the filter_type parameter
             Each bin in the list is the integer ID for 'material', 'surface',
             'cell', 'cellborn', and 'universe' Filters. Each bin is an integer
-            for the cell instance ID for 'distribcell Filters. Each bin is a
+            for the cell instance ID for 'distribcell' Filters. Each bin is a
             2-tuple of floats for 'energy' and 'energyout' filters corresponding
             to the energy boundaries of the bin of interest. Each bin is an
             (x,y,z) 3-tuple for 'mesh' filters corresponding to the mesh cell of
@@ -2755,65 +2795,109 @@ class Tally(object):
             A new tally which encapsulates the sum of data requested.
         """
 
-        # If user did not specify any scores, do not sum across scores
-        if len(scores) == 0:
-            scores = [[]]
-        # Sum across any scores specified by the user
-        else:
-            scores = [[score] for score in scores]
+        # Create new derived Tally for summation
+        tally_sum = Tally()
+        tally_sum._derived = True
+        tally_sum._estimator = self.estimator
+        tally_sum._num_realizations = self.num_realizations
+        tally_sum.with_batch_statistics = self.with_batch_statistics
+        tally_sum._with_summary = self.with_summary
+        tally_sum._sp_filename = self._sp_filename
+        tally_sum._results_read = self._results_read
 
-        # If user did not specify any nuclides, do not sum across nuclides
-        if len(nuclides) == 0:
-            nuclides = [[]]
-        # Sum across any nuclides specified by the user
-        else:
-            nuclides = [[nuclide] for nuclide in nuclides]
+        # Get tally data arrays reshaped with one dimension per filter
+        mean = self.get_reshaped_data(value='mean')
+        std_dev = self.get_reshaped_data(value='std_dev')
 
         # Sum across any filter bins specified by the user
         if filter_type in _FILTER_TYPES:
+            find_filter = self.find_filter(filter_type)
 
             # If user did not specify filter bins, sum across all bins
             if len(filter_bins) == 0:
-                filter = self.find_filter(filter_type)
-                filter_bins = [[(filter.get_bin(i),)] for i in range(filter.num_bins)]
+                bin_indices = np.arange(find_filter.num_bins)
+
+                if filter_type == 'distribcell':
+                    filter_bins = np.arange(find_filter.num_bins)
+                else:
+                    num_bins = find_filter.num_bins
+                    filter_bins = \
+                        [(find_filter.get_bin(i)) for i in range(num_bins)]
+
+            # Only sum across bins specified by the user
             else:
-                filter_bins = [[(filter_bin,)] for filter_bin in filter_bins]
+                bin_indices = \
+                    [find_filter.get_bin_index(bin) for bin in filter_bins]
 
-            filters = [[filter_type]]
-        # If user did not specify a filter type, do not sum across filter bins
+            # Sum across the bins in the user-specified filter
+            for i, self_filter in enumerate(self.filters):
+                if self_filter.type == filter_type:
+                    mean = np.take(mean, indices=bin_indices, axis=i)
+                    std_dev = np.take(std_dev, indices=bin_indices, axis=i)
+                    mean = np.sum(mean, axis=i, keepdims=True)
+                    std_dev = np.sum(std_dev**2, axis=i, keepdims=True)
+                    std_dev = np.sqrt(std_dev)
+
+                    # Add AggregateFilter to the tally sum
+                    if not remove_filter:
+                        filter_sum = \
+                            AggregateFilter(self_filter, filter_bins, 'sum')
+                        tally_sum.add_filter(filter_sum)
+
+                # Add a copy of each filter not summed across to the tally sum
+                else:
+                    tally_sum.add_filter(copy.deepcopy(self_filter))
+
+        # Add a copy of this tally's filters to the tally sum
         else:
-            filter_bins = [[]]
-            filters = [[]]
+            tally_sum._filters = copy.deepcopy(self.filters)
 
-        # Initialize Tally sum
-        tally_sum = 0
+        # Sum across any nuclides specified by the user
+        if len(nuclides) != 0:
+            nuclide_bins = [self.get_nuclide_index(nuclide) for nuclide in nuclides]
+            axis_index = self.num_filters
+            mean = np.take(mean, indices=nuclide_bins, axis=axis_index)
+            std_dev = np.take(std_dev, indices=nuclide_bins, axis=axis_index)
+            mean = np.sum(mean, axis=axis_index, keepdims=True)
+            std_dev = np.sum(std_dev**2, axis=axis_index, keepdims=True)
+            std_dev = np.sqrt(std_dev)
 
-        # Iterate over all Tally slice operands in summation
-        prod = [scores, filters, filter_bins, nuclides]
-        summed_filters = defaultdict(list)
-        for scores, filters, filter_bins, nuclides in itertools.product(*prod):
-            tally_slice = self.get_slice(scores, filters, filter_bins, nuclides)
+            # Add AggregateNuclide to the tally sum
+            nuclide_sum = AggregateNuclide(nuclides, 'sum')
+            tally_sum.add_nuclide(nuclide_sum)
 
-            # Remove filters summed across to avoid bulky CrossFilters
-            if filter_type:
-                filter = tally_slice.find_filter(filter_type)
-                tally_slice.remove_filter(filter)
-                summed_filters[filter_type].append(filter)
-
-            # Accumulate this Tally slice into the Tally sum
-            tally_sum += tally_slice
-
-        # Add back the filter(s) which were summed across to derived tally,
-        # if filter bins were input; otherwise, leave out summed filter(s)
-        if remove_filter and filter_type is not None:
-            # Rename tally sum indicating a summation over a particular filter
-            tally_sum.name = 'sum({0}, {1})'.format(self.name, filter_type)
+        # Add a copy of this tally's nuclides to the tally sum
         else:
-            for summed_filter_type in summed_filters:
-                filters = summed_filters[summed_filter_type]
-                for i in range(1, len(filters)):
-                    filters[i] = CrossFilter(filters[i-1], filters[i], '+')
-                tally_sum.add_filter(filters[-1])
+            tally_sum._nuclides = copy.deepcopy(self.nuclides)
+
+        # Sum across any scores specified by the user
+        if len(scores) != 0:
+            score_bins = [self.get_score_index(score) for score in scores]
+            axis_index = self.num_filters + 1
+            mean = np.take(mean, indices=score_bins, axis=axis_index)
+            std_dev = np.take(std_dev, indices=score_bins, axis=axis_index)
+            mean = np.sum(mean, axis=axis_index, keepdims=True)
+            std_dev = np.sum(std_dev**2, axis=axis_index, keepdims=True)
+            std_dev = np.sqrt(std_dev)
+
+            # Add AggregateScore to the tally sum
+            score_sum = AggregateScore(scores, 'sum')
+            tally_sum.add_score(score_sum)
+
+        # Add a copy of this tally's scores to the tally sum
+        else:
+            tally_sum._scores = copy.deepcopy(self.scores)
+
+        # Update the tally sum's filter strides
+        tally_sum._update_filter_strides()
+
+        # Reshape condensed data arrays with one dimension for all filters
+        mean = np.reshape(mean, tally_sum.shape)
+        std_dev = np.reshape(std_dev, tally_sum.shape)
+
+        # Assign tally sum's data with the new arrays
+        tally_sum._mean = mean
+        tally_sum._std_dev = std_dev
 
         # If original tally was sparse, sparsify the tally summation
         tally_sum.sparse = self.sparse
@@ -2880,11 +2964,8 @@ class Tally(object):
             new_tally._std_dev = np.zeros(new_tally.shape, dtype=np.float64)
             new_tally._std_dev[diag_indices, :, :] = self.std_dev
 
-        # Correct each Filter's stride
-        stride = new_tally.num_nuclides * new_tally.num_scores
-        for filter in reversed(new_tally.filters):
-            filter.stride = stride
-            stride *= filter.num_bins
+        # Update the new tally's filter strides
+        new_tally._update_filter_strides()
 
         # If original tally was sparse, sparsify the diagonalized tally
         new_tally.sparse = self.sparse
diff --git a/tests/test_tally_aggregation/inputs_true.dat b/tests/test_tally_aggregation/inputs_true.dat
new file mode 100644
index 0000000000..7b4276f59e
--- /dev/null
+++ b/tests/test_tally_aggregation/inputs_true.dat
@@ -0,0 +1 @@
+530a5e969901e153531f74aed46246b1e8783a0e2f347e472f7554c9970152f45d85499f17d7df9c35c74fed6f78d449aa70bf0c1f8947cd34d3a829483a0055
\ No newline at end of file
diff --git a/tests/test_tally_aggregation/results_true.dat b/tests/test_tally_aggregation/results_true.dat
new file mode 100644
index 0000000000..cde3e281c2
--- /dev/null
+++ b/tests/test_tally_aggregation/results_true.dat
@@ -0,0 +1 @@
+ba8bfe764fcc0484a4fdab8fdc4ff8ad0e4a98b1ff33e8687899c8cc6bf80cb28b3a59aeaec84bd74681b8b5f19f714292ccaa9c9d4ba852b2cc29872f612e10
\ No newline at end of file
diff --git a/tests/test_tally_aggregation/test_tally_aggregation.py b/tests/test_tally_aggregation/test_tally_aggregation.py
new file mode 100644
index 0000000000..a5c8d94141
--- /dev/null
+++ b/tests/test_tally_aggregation/test_tally_aggregation.py
@@ -0,0 +1,99 @@
+#!/usr/bin/env python
+
+import os
+import sys
+import glob
+import hashlib
+sys.path.insert(0, os.pardir)
+from testing_harness import PyAPITestHarness
+import openmc
+
+
+class TallyAggregationTestHarness(PyAPITestHarness):
+    def _build_inputs(self):
+
+        # The summary.h5 file needs to be created to read in the tallies
+        self._input_set.settings.output = {'summary': True}
+
+        # Initialize the tallies file
+        tallies_file = openmc.TalliesFile()
+
+        # Initialize the nuclides
+        u235 = openmc.Nuclide('U-235')
+        u238 = openmc.Nuclide('U-238')
+        pu239 = openmc.Nuclide('Pu-239')
+
+        # Initialize the filters
+        energy_filter = openmc.Filter(type='energy', bins=[0.0, 0.253e-6,
+                                                           1.0e-3, 1.0, 20.0])
+        distrib_filter   = openmc.Filter(type='distribcell', bins=[60])
+
+        # Initialized the tallies
+        tally = openmc.Tally(name='distribcell tally')
+        tally.add_filter(energy_filter)
+        tally.add_filter(distrib_filter)
+        tally.add_score('nu-fission')
+        tally.add_score('total')
+        tally.add_nuclide(u235)
+        tally.add_nuclide(u238)
+        tally.add_nuclide(pu239)
+        tallies_file.add_tally(tally)
+
+        # Export tallies to file
+        self._input_set.tallies = tallies_file
+        super(TallyAggregationTestHarness, self)._build_inputs()
+
+    def _get_results(self, hash_output=True):
+        """Digest info in the statepoint and return as a string."""
+
+        # Read the statepoint file.
+        statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0]
+        sp = openmc.StatePoint(statepoint)
+
+        # Read the summary file.
+        summary = glob.glob(os.path.join(os.getcwd(), 'summary.h5'))[0]
+        su = openmc.Summary(summary)
+        sp.link_with_summary(su)
+
+        # Extract the tally of interest
+        tally = sp.get_tally(name='distribcell tally')
+
+        # Perform tally aggregations across filter bins, nuclides and scores
+        outstr = ''
+
+        # Sum across all energy filter bins
+        tally_sum = tally.summation(filter_type='energy')
+        outstr += ', '.join(map(str, tally_sum.mean))
+        outstr += ', '.join(map(str, tally_sum.std_dev))
+
+        # Sum across all distribcell filter bins
+        tally_sum = tally.summation(filter_type='distribcell')
+        outstr += ', '.join(map(str, tally_sum.mean))
+        outstr += ', '.join(map(str, tally_sum.std_dev))
+
+        # Sum across all nuclides
+        tally_sum = tally.summation(nuclides=['U-235', 'U-238', 'Pu-239'])
+        outstr += ', '.join(map(str, tally_sum.mean))
+        outstr += ', '.join(map(str, tally_sum.std_dev))
+
+        # Sum across all scores
+        tally_sum = tally.summation(scores=['nu-fission', 'total'])
+        outstr += ', '.join(map(str, tally_sum.mean))
+        outstr += ', '.join(map(str, tally_sum.std_dev))
+
+        # Hash the results if necessary
+        if hash_output:
+            sha512 = hashlib.sha512()
+            sha512.update(outstr.encode('utf-8'))
+            outstr = sha512.hexdigest()
+
+        return outstr
+
+    def _cleanup(self):
+        super(TallyAggregationTestHarness, self)._cleanup()
+        f = os.path.join(os.getcwd(), 'tallies.xml')
+        if os.path.exists(f): os.remove(f)
+
+if __name__ == '__main__':
+    harness = TallyAggregationTestHarness('statepoint.10.h5', True)
+    harness.main()
diff --git a/tests/test_tally_arithmetic/test_tally_arithmetic.py b/tests/test_tally_arithmetic/test_tally_arithmetic.py
index 6643fd5a22..1954334b7f 100644
--- a/tests/test_tally_arithmetic/test_tally_arithmetic.py
+++ b/tests/test_tally_arithmetic/test_tally_arithmetic.py
@@ -100,8 +100,6 @@ class TallyArithmeticTestHarness(PyAPITestHarness):
                                          'entrywise')
         outstr += str(tally_3.mean)
 
-        print(outstr)
-
         # Hash the results if necessary
         if hash_output:
             sha512 = hashlib.sha512()