diff --git a/openmc_dagmc_wrapper/Tally.py b/openmc_dagmc_wrapper/Tally.py
index 568cb60..901dc51 100644
--- a/openmc_dagmc_wrapper/Tally.py
+++ b/openmc_dagmc_wrapper/Tally.py
@@ -6,7 +6,7 @@
 from openmc.data import REACTION_MT, REACTION_NAME
 
 from openmc_dagmc_wrapper import Materials
-from .utils import create_material, silently_remove_file, find_bounding_box
+from .utils import find_bounding_box
 
 
 class Tally(openmc.Tally):
@@ -33,7 +33,10 @@ def tally_type(self, value):
         output_options = (
             [
                 "TBR",
+                "flux",
                 "heating",
+                "photon_heating",
+                "neutron_heating",
                 "neutron_flux",
                 "photon_flux",
                 "absorption",
@@ -58,6 +61,7 @@ def tally_type(self, value):
 
     def set_score(self):
         flux_scores = [
+            "flux",
             "neutron_flux",
             "photon_flux",
             "neutron_fast_flux",
@@ -454,7 +458,17 @@ def compute_filters(tally_type):
     neutron_particle_filter = openmc.ParticleFilter(["neutron"])
 
     additional_filters = []
-    if tally_type == "neutron_fast_flux":
+    if tally_type == "neutron_flux":
+        additional_filters = [neutron_particle_filter]
+    elif tally_type == "photon_flux":
+        additional_filters = [photon_particle_filter]
+
+    elif tally_type == "neutron_heating":
+        additional_filters = [neutron_particle_filter]
+    elif tally_type == "photon_heating":
+        additional_filters = [photon_particle_filter]
+
+    elif tally_type == "neutron_fast_flux":
         energy_bins = [1e6, 1000e6]
         energy_filter = openmc.EnergyFilter(energy_bins)
         additional_filters = [neutron_particle_filter, energy_filter]
@@ -462,6 +476,7 @@ def compute_filters(tally_type):
         energy_bins = [1e6, 1000e6]
         energy_filter = openmc.EnergyFilter(energy_bins)
         additional_filters = [photon_particle_filter, energy_filter]
+
     elif tally_type == "neutron_spectra":
         energy_bins = openmc.mgxs.GROUP_STRUCTURES["CCFE-709"]
         energy_filter = openmc.EnergyFilter(energy_bins)
@@ -470,6 +485,7 @@ def compute_filters(tally_type):
         energy_bins = openmc.mgxs.GROUP_STRUCTURES["CCFE-709"]
         energy_filter = openmc.EnergyFilter(energy_bins)
         additional_filters = [photon_particle_filter, energy_filter]
+
     elif tally_type == "neutron_effective_dose":
         energy_function_filter_n = openmc.EnergyFunctionFilter(
             energy_bins_n, dose_coeffs_n
diff --git a/tests/test_tallies/test_mesh_tally_2d.py b/tests/test_tallies/test_mesh_tally_2d.py
index c488091..5981523 100644
--- a/tests/test_tallies/test_mesh_tally_2d.py
+++ b/tests/test_tallies/test_mesh_tally_2d.py
@@ -53,19 +53,19 @@ def test_shape_of_resulting_png(self):
             },
         )
         tally1 = odw.MeshTally2D(
-            tally_type="flux",
+            tally_type="neutron_flux",
             plane="xy",
             bounding_box=self.h5m_filename_smaller,
             mesh_resolution=(10, 200)
         )
         tally2 = odw.MeshTally2D(
-            tally_type="flux",
+            tally_type="neutron_flux",
             plane="xz",
             bounding_box=self.h5m_filename_smaller,
             mesh_resolution=(20, 100)
         )
         tally3 = odw.MeshTally2D(
-            tally_type="flux",
+            tally_type="neutron_flux",
             plane="yz",
             bounding_box=self.h5m_filename_smaller,
             mesh_resolution=(30, 500)
@@ -88,6 +88,6 @@ def test_shape_of_resulting_png(self):
 
         odw.process_results(statepoint_file, fusion_power=1e9)
 
-        assert Path('flux_on_2D_mesh_xy.png').exists()
-        assert Path('flux_on_2D_mesh_xz.png').exists()
-        assert Path('flux_on_2D_mesh_yz.png').exists()
+        assert Path('neutron_flux_on_2D_mesh_xy.png').exists()
+        assert Path('neutron_flux_on_2D_mesh_xz.png').exists()
+        assert Path('neutron_flux_on_2D_mesh_yz.png').exists()