diff --git a/tests/regression_tests/photon_production/inputs_true.dat b/tests/regression_tests/photon_production/inputs_true.dat
index 1821bea92f..31eab718b8 100644
--- a/tests/regression_tests/photon_production/inputs_true.dat
+++ b/tests/regression_tests/photon_production/inputs_true.dat
@@ -1,16 +1,16 @@
 <?xml version='1.0' encoding='utf-8'?>
 <geometry>
-  <cell id="13" material="void" region="-9 10 -11" universe="9" />
-  <cell id="14" material="13" region="-9 11 -12" universe="9" />
-  <cell id="15" material="void" region="~(-9 10 -12)" universe="9" />
-  <surface boundary="vacuum" coeffs="0.0 0.0 1.0" id="9" type="x-cylinder" />
-  <surface boundary="vacuum" coeffs="-1.0" id="10" type="x-plane" />
-  <surface coeffs="1.0" id="11" type="x-plane" />
-  <surface boundary="vacuum" coeffs="1000000000.0" id="12" type="x-plane" />
+  <cell id="1" material="void" region="-1 2 -3" universe="1" />
+  <cell id="2" material="1" region="-1 3 -4" universe="1" />
+  <cell id="3" material="void" region="~(-1 2 -4)" universe="1" />
+  <surface boundary="vacuum" coeffs="0.0 0.0 1.0" id="1" type="x-cylinder" />
+  <surface boundary="vacuum" coeffs="-1.0" id="2" type="x-plane" />
+  <surface coeffs="1.0" id="3" type="x-plane" />
+  <surface boundary="vacuum" coeffs="1000000000.0" id="4" type="x-plane" />
 </geometry>
 <?xml version='1.0' encoding='utf-8'?>
 <materials>
-  <material id="13">
+  <material id="1">
     <density units="g/cm3" value="2.6989" />
     <nuclide ao="1.0" name="Al27" />
   </material>
@@ -38,7 +38,7 @@
 <?xml version='1.0' encoding='utf-8'?>
 <tallies>
   <filter id="1" type="surface">
-    <bins>9</bins>
+    <bins>1</bins>
   </filter>
   <filter id="2" type="particle">
     <bins>photon</bins>
diff --git a/tests/regression_tests/photon_production/test.py b/tests/regression_tests/photon_production/test.py
index 1f645ebb94..4ed0a7ba72 100644
--- a/tests/regression_tests/photon_production/test.py
+++ b/tests/regression_tests/photon_production/test.py
@@ -1,90 +1,72 @@
-from math import pi
-
-import numpy as np
 import openmc
+import pytest
 
 from tests.testing_harness import PyAPITestHarness
 
 
-class SourceTestHarness(PyAPITestHarness):
-    def _build_inputs(self):
-        mat = openmc.Material()
-        mat.set_density('g/cm3', 2.6989)
-        mat.add_nuclide('Al27', 1.0)
-        materials = openmc.Materials([mat])
-        materials.export_to_xml()
+@pytest.fixture
+def model():
+    model = openmc.model.Model()
+    mat = openmc.Material()
+    mat.set_density('g/cm3', 2.6989)
+    mat.add_nuclide('Al27', 1.0)
+    model.materials.append(mat)
 
-        cyl = openmc.XCylinder(r=1.0, boundary_type='vacuum')
-        x_plane_left = openmc.XPlane(-1.0, 'vacuum')
-        x_plane_center = openmc.XPlane(1.0)
-        x_plane_right = openmc.XPlane(1.0e9, 'vacuum')
+    cyl = openmc.XCylinder(r=1.0, boundary_type='vacuum')
+    x_plane_left = openmc.XPlane(-1.0, 'vacuum')
+    x_plane_center = openmc.XPlane(1.0)
+    x_plane_right = openmc.XPlane(1.0e9, 'vacuum')
 
-        inner_cyl_left = openmc.Cell()
-        inner_cyl_right = openmc.Cell()
-        outer_cyl = openmc.Cell()
+    inner_cyl_left = openmc.Cell()
+    inner_cyl_right = openmc.Cell()
+    outer_cyl = openmc.Cell()
 
-        inner_cyl_left.region = -cyl & +x_plane_left & -x_plane_center
-        inner_cyl_right.region = -cyl & +x_plane_center & -x_plane_right
-        outer_cyl.region = ~(-cyl & +x_plane_left & -x_plane_right)
-        inner_cyl_right.fill = mat
-        geometry = openmc.Geometry([inner_cyl_left, inner_cyl_right, outer_cyl])
-        geometry.export_to_xml()
+    inner_cyl_left.region = -cyl & +x_plane_left & -x_plane_center
+    inner_cyl_right.region = -cyl & +x_plane_center & -x_plane_right
+    outer_cyl.region = ~(-cyl & +x_plane_left & -x_plane_right)
+    inner_cyl_right.fill = mat
+    model.geometry = openmc.Geometry([inner_cyl_left, inner_cyl_right, outer_cyl])
 
-        source = openmc.Source()
-        source.space = openmc.stats.Point((0, 0, 0))
-        source.angle = openmc.stats.Monodirectional()
-        source.energy = openmc.stats.Discrete([14.0e6], [1.0])
-        source.particle = 'neutron'
+    source = openmc.Source()
+    source.space = openmc.stats.Point((0, 0, 0))
+    source.angle = openmc.stats.Monodirectional()
+    source.energy = openmc.stats.Discrete([14.0e6], [1.0])
+    source.particle = 'neutron'
 
-        settings = openmc.Settings()
-        settings.particles = 10000
-        settings.run_mode = 'fixed source'
-        settings.batches = 1
-        settings.photon_transport = True
-        settings.electron_treatment = 'ttb'
-        settings.cutoff = {'energy_photon' : 1000.0}
-        settings.source = source
-        settings.export_to_xml()
+    model.settings.particles = 10000
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 1
+    model.settings.photon_transport = True
+    model.settings.electron_treatment = 'ttb'
+    model.settings.cutoff = {'energy_photon' : 1000.0}
+    model.settings.source = source
 
-        surface_filter = openmc.SurfaceFilter(cyl)
-        particle_filter = openmc.ParticleFilter('photon')
-        current_tally = openmc.Tally()
-        current_tally.filters = [surface_filter, particle_filter]
-        current_tally.scores = ['current']
-        tally_tracklength = openmc.Tally()
-        tally_tracklength.filters = [particle_filter]
-        tally_tracklength.scores = ['total', 'heating']
-        tally_tracklength.nuclides = ['Al27', 'total']
-        tally_tracklength.estimator = 'tracklength'
-        tally_collision = openmc.Tally()
-        tally_collision.filters = [particle_filter]
-        tally_collision.scores = ['total', 'heating']
-        tally_collision.nuclides = ['Al27', 'total']
-        tally_collision.estimator = 'collision'
-        tally_analog = openmc.Tally()
-        tally_analog.filters = [particle_filter]
-        tally_analog.scores = ['total', 'heating']
-        tally_analog.nuclides = ['Al27', 'total']
-        tally_analog.estimator = 'analog'
-        tallies = openmc.Tallies([current_tally, tally_tracklength,
-                                  tally_collision, tally_analog])
-        tallies.export_to_xml()
+    surface_filter = openmc.SurfaceFilter(cyl)
+    particle_filter = openmc.ParticleFilter('photon')
+    current_tally = openmc.Tally()
+    current_tally.filters = [surface_filter, particle_filter]
+    current_tally.scores = ['current']
+    tally_tracklength = openmc.Tally()
+    tally_tracklength.filters = [particle_filter]
+    tally_tracklength.scores = ['total', 'heating']
+    tally_tracklength.nuclides = ['Al27', 'total']
+    tally_tracklength.estimator = 'tracklength'
+    tally_collision = openmc.Tally()
+    tally_collision.filters = [particle_filter]
+    tally_collision.scores = ['total', 'heating']
+    tally_collision.nuclides = ['Al27', 'total']
+    tally_collision.estimator = 'collision'
+    tally_analog = openmc.Tally()
+    tally_analog.filters = [particle_filter]
+    tally_analog.scores = ['total', 'heating']
+    tally_analog.nuclides = ['Al27', 'total']
+    tally_analog.estimator = 'analog'
+    model.tallies.extend([current_tally, tally_tracklength,
+                          tally_collision, tally_analog])
 
-    def _get_results(self):
-        with openmc.StatePoint(self._sp_name) as sp:
-            outstr = ''
-            for i, tally_ind in enumerate(sp.tallies):
-                tally = sp.tallies[tally_ind]
-                results = np.zeros((tally.sum.size * 2, ))
-                results[0::2] = tally.sum.ravel()
-                results[1::2] = tally.sum_sq.ravel()
-                results = ['{0:12.6E}'.format(x) for x in results]
-
-                outstr += 'tally {}:\n'.format(i + 1)
-                outstr += '\n'.join(results) + '\n'
-            return outstr
+    return model
 
 
-def test_photon_production():
-    harness = SourceTestHarness('statepoint.1.h5')
+def test_photon_production(model):
+    harness = PyAPITestHarness('statepoint.1.h5', model)
     harness.main()