diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h
index e8504d48261..5f96751e86b 100644
--- a/include/openmc/geometry.h
+++ b/include/openmc/geometry.h
@@ -8,6 +8,7 @@
 
 #include "openmc/array.h"
 #include "openmc/constants.h"
+#include "openmc/position.h"
 #include "openmc/random_ray/source_region.h" // For hash_combine
 #include "openmc/vector.h"
 
@@ -86,6 +87,23 @@ int check_cell_overlap(GeometryState& p, bool error = true);
 
 int cell_instance_at_level(const GeometryState& p, int level);
 
+//==============================================================================
+//! Rotate a direction from a local coordinate frame into the root frame
+//!
+//! Surface and lattice normals are expressed in the local coordinate frame of
+//! the universe that contains them. Comparing such a normal against the
+//! particle's direction of travel (which is stored in the root frame at
+//! coordinate level zero) requires undoing the rotations that were applied
+//! while descending to \c level.
+//!
+//! \param p A particle whose coordinate levels give the chain of rotations
+//! \param level The level (zero indexed) that \c u is expressed in
+//! \param u A direction in the local frame of \c level
+//! \return The same direction expressed in the root coordinate frame
+//==============================================================================
+
+Direction rotate_to_root(const GeometryState& p, int level, Direction u);
+
 //==============================================================================
 //! Locate a particle in the geometry tree and set its geometry data fields.
 //!
diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h
index 24b7eb53c0a..65266fdb8a2 100644
--- a/include/openmc/particle_data.h
+++ b/include/openmc/particle_data.h
@@ -384,6 +384,14 @@ class GeometryState {
   int& surface() { return surface_; }
   const int& surface() const { return surface_; }
 
+  // Outward unit normal of the surface (or lattice boundary) currently being
+  // crossed, expressed in the root coordinate frame. Set by
+  // score_surface_tally() immediately before the tally filters are evaluated,
+  // so that filters can compare it against the root-frame direction u()
+  // without having to know which coordinate level the surface lives in.
+  Direction& surface_normal() { return surface_normal_; }
+  const Direction& surface_normal() const { return surface_normal_; }
+
   // Surface index based on the current value of the surface_ attribute
   int surface_index() const
   {
@@ -436,6 +444,9 @@ class GeometryState {
   int surface_ {
     SURFACE_NONE}; //!< surface token for surface the particle is currently on
 
+  //! Outward normal of the surface being crossed, in the root coordinate frame
+  Direction surface_normal_ {0.0, 0.0, 1.0};
+
   BoundaryInfo boundary_; //!< Info about the next intersection
 
   int material_ {-1};      //!< index for current material
diff --git a/include/openmc/surface.h b/include/openmc/surface.h
index 2d8580345a4..d6b9ce956fe 100644
--- a/include/openmc/surface.h
+++ b/include/openmc/surface.h
@@ -41,6 +41,21 @@ class Surface {
   unique_ptr<BoundaryCondition> bc_; //!< Boundary condition
   bool surf_source_ {false}; //!< Activate source banking for the surface?
 
+  //! Is this surface used only by cells in the root universe?
+  //!
+  //! evaluate() and normal() work in the local coordinate frame of the
+  //! universe holding the surface. That frame coincides with the root (lab)
+  //! frame only for surfaces in the root universe, so this flag marks the
+  //! surfaces for which a lab-frame position may be passed to them directly.
+  //! Set by finalize_geometry().
+  //!
+  //! Only meaningful for CSG surfaces. finalize_geometry() determines it by
+  //! walking the surfaces named in each cell's region, and DAGCell does not
+  //! report any (it does not override Cell::surfaces()), so a DAGMC surface
+  //! keeps the default here even when its universe is nested below the root
+  //! and transformed. Check geom_type() before relying on this flag.
+  bool root_frame_ {true};
+
   explicit Surface(pugi::xml_node surf_node);
   Surface();
 
diff --git a/src/geometry.cpp b/src/geometry.cpp
index ecacf0bffbf..34a4bc10b53 100644
--- a/src/geometry.cpp
+++ b/src/geometry.cpp
@@ -124,6 +124,23 @@ int cell_instance_at_level(const GeometryState& p, int level)
 
 //==============================================================================
 
+Direction rotate_to_root(const GeometryState& p, int level, Direction u)
+{
+  // Each coordinate level below the root was reached by applying the rotation
+  // matrix of the cell one level above it, so walk back up applying the
+  // inverse of each rotation in turn. Translations are irrelevant here since
+  // they do not affect directions.
+  for (int i = level; i > 0; --i) {
+    if (p.coord(i).rotated()) {
+      const auto& c {*model::cells[p.coord(i - 1).cell()]};
+      u = u.inverse_rotate(c.rotation_);
+    }
+  }
+  return u;
+}
+
+//==============================================================================
+
 bool find_cell_inner(
   GeometryState& p, const NeighborList* neighbor_list, bool verbose)
 {
diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp
index a740740c1e6..f8139326b03 100644
--- a/src/geometry_aux.cpp
+++ b/src/geometry_aux.cpp
@@ -282,6 +282,19 @@ void finalize_geometry()
 
   // Determine number of nested coordinate levels in the geometry
   model::n_coord_levels = maximum_levels(model::root_universe);
+
+  // Flag surfaces that are used outside of the root universe. Their local
+  // coordinate frame does not coincide with the root frame, so a lab-frame
+  // position cannot be handed to Surface::evaluate() or Surface::normal().
+  // Note that DAGCell reports no surfaces, so this leaves DAGMC surfaces at
+  // their default; see the comment on Surface::root_frame_.
+  for (const auto& c : model::cells) {
+    if (c->universe_ == model::root_universe)
+      continue;
+    for (auto token : c->surfaces()) {
+      model::surfaces[std::abs(token) - 1]->root_frame_ = false;
+    }
+  }
 }
 
 //==============================================================================
diff --git a/src/particle.cpp b/src/particle.cpp
index 98270af56fb..6a8e5cd4504 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -344,32 +344,60 @@ void Particle::event_cross_surface()
   surface() = boundary().surface();
   n_coord() = boundary().coord_level();
 
+  // The surface or lattice being crossed belongs to the universe at the
+  // coordinate level the boundary search found it on, so its normal is
+  // reported in that level's local frame while the particle direction used to
+  // score surface tallies lives in the root frame. The normal therefore has to
+  // be evaluated at the local position and rotated up into the root frame, and
+  // that has to happen before the crossing is carried out, since crossing
+  // invalidates the coordinate levels.
+  //
+  // Take the level from the boundary rather than from n_coord(). The two are
+  // equal here because of the assignment just above, but reading it from the
+  // boundary keeps this independent of that, and matches where the level came
+  // from originally.
+  int i_surf_level = boundary().coord_level() - 1;
+
   if (boundary().lattice_translation()[0] != 0 ||
       boundary().lattice_translation()[1] != 0 ||
       boundary().lattice_translation()[2] != 0) {
     // Particle crosses lattice boundary
 
-    int i_lattice = coord(boundary().coord_level() - 1).lattice();
-    bool verbose = settings::verbosity >= 10 || trace();
-    cross_lattice(*this, boundary(), verbose);
-    event() = TallyEvent::LATTICE;
+    int i_lattice = coord(i_surf_level).lattice();
 
-    // Score cell to cell partial currents
+    // Determine the lattice boundary normal in the root frame before crossing
+    bool normal_is_valid = false;
+    Direction normal;
     if (!model::active_surface_tallies.empty()) {
       auto& lat {*model::lattices[i_lattice]};
       bool is_valid;
-      Direction normal =
-        lat.get_normal(boundary().lattice_translation(), is_valid);
+      normal = lat.get_normal(boundary().lattice_translation(), is_valid);
       if (is_valid) {
-        normal /= normal.norm();
-        score_surface_tally(*this, model::active_surface_tallies, normal);
+        normal = rotate_to_root(*this, i_surf_level, normal / normal.norm());
+        normal_is_valid = true;
       }
     }
 
+    bool verbose = settings::verbosity >= 10 || trace();
+    cross_lattice(*this, boundary(), verbose);
+    event() = TallyEvent::LATTICE;
+
+    // Score cell to cell partial currents
+    if (normal_is_valid) {
+      score_surface_tally(*this, model::active_surface_tallies, normal);
+    }
+
   } else {
 
     const auto& surf {*model::surfaces[surface_index()].get()};
 
+    // Determine the surface normal in the root frame before crossing
+    Direction normal;
+    if (!model::active_surface_tallies.empty()) {
+      normal = surf.normal(coord(i_surf_level).r());
+      normal = rotate_to_root(*this, i_surf_level, normal / normal.norm());
+    }
+
     // Particle crosses surface
     // If BC, add particle to surface source before crossing surface
     if (surf.surf_source_ && surf.bc_) {
@@ -387,8 +415,6 @@ void Particle::event_cross_surface()
 
     // Score cell to cell partial currents
     if (!model::active_surface_tallies.empty()) {
-      Direction normal = surf.normal(r());
-      normal /= normal.norm();
       score_surface_tally(*this, model::active_surface_tallies, normal);
     }
   }
diff --git a/src/plot.cpp b/src/plot.cpp
index 68ce3c7d167..96711191b6e 100644
--- a/src/plot.cpp
+++ b/src/plot.cpp
@@ -1825,12 +1825,7 @@ void PhongRay::on_intersection()
 
     // Need to apply rotations to find the normal vector in
     // the base level universe's coordinate system.
-    for (int lev = surf_level - 1; lev >= 0; --lev) {
-      if (coord(lev + 1).rotated()) {
-        const Cell& c {*model::cells[coord(lev).cell()]};
-        normal = normal.inverse_rotate(c.rotation_);
-      }
-    }
+    normal = rotate_to_root(*this, surf_level, normal);
 
     // use the normal opposed to the ray direction
     if (normal.dot(u()) > 0.0) {
diff --git a/src/source.cpp b/src/source.cpp
index f951722cc18..d8528d56105 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -554,11 +554,19 @@ SourceSite FileSource::sample(uint64_t* seed) const
   // surface containing the source site, determine the signed half-space from
   // the particle direction. Otherwise, ignore the surface ID and allow the
   // normal cell search to locate the particle.
+  //
+  // The site position and direction are in the root coordinate frame, whereas
+  // evaluate() and normal() work in the local frame of the universe holding
+  // the surface. The half-space can therefore only be recovered here for
+  // surfaces in the root universe; for any other surface the frames differ by
+  // an unknown transform (a universe may be filled in several places with
+  // different rotations, so the transform cannot be recovered from the site
+  // alone) and the surface ID is dropped rather than signed incorrectly.
   if (site.surf_id != SURFACE_NONE) {
     auto it = model::surface_map.find(std::abs(site.surf_id));
     if (it != model::surface_map.end()) {
       const auto& surf = *model::surfaces[it->second];
-      if (surf.geom_type() == GeometryType::CSG &&
+      if (surf.geom_type() == GeometryType::CSG && surf.root_frame_ &&
           std::abs(surf.evaluate(site.r)) < FP_COINCIDENT) {
         int surf_id = std::abs(site.surf_id);
         site.surf_id =
diff --git a/src/tallies/filter_musurface.cpp b/src/tallies/filter_musurface.cpp
index 340149d4cff..2ab90b0aac3 100644
--- a/src/tallies/filter_musurface.cpp
+++ b/src/tallies/filter_musurface.cpp
@@ -3,7 +3,6 @@
 #include <cmath> // for abs, copysign
 
 #include "openmc/search.h"
-#include "openmc/surface.h"
 #include "openmc/tallies/tally_scoring.h"
 
 namespace openmc {
@@ -11,10 +10,11 @@ namespace openmc {
 void MuSurfaceFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
-  // Get surface normal (and make sure it is a unit vector)
-  const auto surf {model::surfaces[p.surface_index()].get()};
-  auto n = surf->normal(p.r());
-  n /= n.norm();
+  // Use the normal recorded for the crossing being scored. It is already a
+  // unit vector expressed in the root coordinate frame, which is the frame
+  // p.u() is in -- recomputing it here from the surface would give the normal
+  // in the local frame of whichever universe holds the surface.
+  Direction n = p.surface_normal();
 
   // Determine whether normal should be pointing in or out
   if (p.surface() < 0)
diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp
index 27e6671833a..f213ec51a65 100644
--- a/src/tallies/tally_scoring.cpp
+++ b/src/tallies/tally_scoring.cpp
@@ -2662,6 +2662,11 @@ void score_surface_tally(
 {
   double wgt = p.wgt_last();
 
+  // Make the normal available to filters (e.g. MuSurfaceFilter) that need it.
+  // The caller is responsible for supplying it in the root coordinate frame so
+  // that it can be compared directly against p.u().
+  p.surface_normal() = normal;
+
   double mu = std::clamp(p.u().dot(normal), -1.0, 1.0);
 
   // Sign for net current: +1 if crossing outward (in direction of normal),
diff --git a/tests/unit_tests/test_filter_musurface.py b/tests/unit_tests/test_filter_musurface.py
index ca0db71f0c6..745b4002615 100644
--- a/tests/unit_tests/test_filter_musurface.py
+++ b/tests/unit_tests/test_filter_musurface.py
@@ -1,3 +1,5 @@
+import math
+
 import openmc
 
 
@@ -36,3 +38,49 @@ def test_musurface(run_in_tmpdir):
         assert element == 0.0
 
 
+def test_musurface_rotated_universe(run_in_tmpdir):
+    """MuSurfaceFilter uses the surface normal in the root coordinate frame.
+
+    The plane lives in a universe filled into a cell rotated 45 degrees about
+    z, so its normal in the plane's own frame is (1, 0, 0) while the particle
+    direction is stored in the root frame. Binning must use the root-frame
+    normal, giving mu = cos(45 deg) rather than 1.
+    """
+    openmc.reset_auto_ids()
+
+    xplane = openmc.XPlane(0.0)
+    inner1 = openmc.Cell(region=-xplane)
+    inner2 = openmc.Cell(region=+xplane)
+    inner_univ = openmc.Universe(cells=[inner1, inner2])
+
+    sph = openmc.Sphere(r=10.0, boundary_type='vacuum')
+    root_cell = openmc.Cell(region=-sph, fill=inner_univ)
+    root_cell.rotation = (0.0, 0.0, 45.0)
+
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([root_cell])
+
+    src = openmc.IndependentSource()
+    src.space = openmc.stats.Point((-5.0, 0.0, 0.0))
+    src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0))
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 1
+    model.settings.particles = 100
+    model.settings.source = src
+
+    # 20 equal-width bins from -1 to 1; cos(45 deg) = 0.7071 falls in [0.7, 0.8)
+    tally = openmc.Tally()
+    tally.filters = [
+        openmc.MuSurfaceFilter(20),
+        openmc.SurfaceFilter([xplane]),
+    ]
+    tally.scores = ['current']
+    model.tallies = [tally]
+
+    model.run(apply_tally_results=True)
+    current_mu = tally.mean.ravel()
+
+    expected_bin = int((math.cos(math.radians(45.0)) + 1.0) / 0.1)
+    assert current_mu[expected_bin] == 1.0
+    assert current_mu.sum() == 1.0
diff --git a/tests/unit_tests/test_source_file.py b/tests/unit_tests/test_source_file.py
index f19ea74a67a..01debdd10a9 100644
--- a/tests/unit_tests/test_source_file.py
+++ b/tests/unit_tests/test_source_file.py
@@ -169,3 +169,108 @@ def test_source_file_photon_transport(run_in_tmpdir):
 
     # Running OpenMC should succeed
     model.run()
+
+
+def _rotated_universe_model(angle):
+    """Plane inside a universe filled into a cell rotated by `angle` about z.
+
+    Particles are born at (-5, 0, 0) travelling along +x and cross the plane
+    at the origin. The plane's normal in its own frame is +x, but in the root
+    frame it is rotated by `angle`, so for angle > 90 degrees the particle
+    enters the negative half-space even though its root-frame direction has a
+    positive dot product with the unrotated normal.
+    """
+    openmc.reset_auto_ids()
+
+    xplane = openmc.XPlane(0.0, surface_id=7)
+    inner1 = openmc.Cell(cell_id=1, region=-xplane)
+    inner2 = openmc.Cell(cell_id=2, region=+xplane)
+    inner_univ = openmc.Universe(cells=[inner1, inner2])
+
+    sph = openmc.Sphere(r=10.0, boundary_type='vacuum')
+    root_cell = openmc.Cell(cell_id=3, region=-sph, fill=inner_univ)
+    root_cell.rotation = (0.0, 0.0, angle)
+
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([root_cell])
+
+    src = openmc.IndependentSource()
+    src.space = openmc.stats.Point((-5.0, 0.0, 0.0))
+    src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0))
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 1
+    model.settings.particles = 20
+    model.settings.source = src
+    return model, inner1, inner2
+
+
+def test_surface_source_half_space_rotated_universe(run_in_tmpdir):
+    """A surface source on a surface in a rotated universe starts in the
+    correct half-space.
+
+    Surface source files store unsigned surface IDs, and the half-space is
+    recovered when the file is read back. That recovery compares the root-frame
+    particle direction against Surface::normal(), which is expressed in the
+    local frame of the universe holding the surface. With a 135 degree
+    rotation the two disagree in sign, so the particle used to be started in
+    the cell on the wrong side of the plane.
+    """
+    model, _, _ = _rotated_universe_model(135.0)
+    model.settings.surf_source_write = {
+        'surface_ids': [7], 'max_particles': 50}
+    model.run()
+
+    # Read the surface source back in and see which cell the particles start in
+    model, inner1, inner2 = _rotated_universe_model(135.0)
+    model.settings.source = openmc.FileSource('surface_source.h5')
+    tally = openmc.Tally()
+    tally.filters = [openmc.CellFilter([inner1, inner2])]
+    tally.scores = ['flux']
+    model.tallies = [tally]
+    model.run(apply_tally_results=True)
+
+    # The particle direction in the plane's own frame points into the negative
+    # half-space, so all of the track length belongs to the -x cell.
+    assert tally.mean.flat[0] == pytest.approx(10.0)
+    assert tally.mean.flat[1] == pytest.approx(0.0)
+
+
+def test_surface_source_half_space_root_universe(run_in_tmpdir):
+    """The half-space is still recovered for a surface in the root universe."""
+    openmc.reset_auto_ids()
+
+    xplane = openmc.XPlane(0.0, surface_id=7)
+    sph = openmc.Sphere(r=10.0, boundary_type='vacuum')
+
+    def build():
+        cell1 = openmc.Cell(cell_id=1, region=-sph & -xplane)
+        cell2 = openmc.Cell(cell_id=2, region=-sph & +xplane)
+        model = openmc.Model()
+        model.geometry = openmc.Geometry([cell1, cell2])
+        src = openmc.IndependentSource()
+        src.space = openmc.stats.Point((-5.0, 0.0, 0.0))
+        src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0))
+        model.settings.run_mode = 'fixed source'
+        model.settings.batches = 1
+        model.settings.particles = 20
+        model.settings.source = src
+        return model, cell1, cell2
+
+    model, _, _ = build()
+    model.settings.surf_source_write = {
+        'surface_ids': [7], 'max_particles': 50}
+    model.run()
+
+    openmc.reset_auto_ids()
+    model, cell1, cell2 = build()
+    model.settings.source = openmc.FileSource('surface_source.h5')
+    tally = openmc.Tally()
+    tally.filters = [openmc.CellFilter([cell1, cell2])]
+    tally.scores = ['flux']
+    model.tallies = [tally]
+    model.run(apply_tally_results=True)
+
+    # Travelling along +x through the plane, the particle enters the +x cell
+    assert tally.mean.flat[0] == pytest.approx(0.0)
+    assert tally.mean.flat[1] == pytest.approx(10.0)
diff --git a/tests/unit_tests/test_surface_flux.py b/tests/unit_tests/test_surface_flux.py
index b51c55fc3f4..4a9cff6d77f 100644
--- a/tests/unit_tests/test_surface_flux.py
+++ b/tests/unit_tests/test_surface_flux.py
@@ -204,3 +204,103 @@ def test_surface_filter_do_not_tally_virtual_surface_crossing(run_in_tmpdir):
 
     # Every particle crosses exit the cube with weight 1, so current = 1.0
     assert current_mean.sum() == pytest.approx(1.0, rel=1e-8)
+
+
+def test_surface_flux_rotated_universe(run_in_tmpdir):
+    """Surface flux uses the surface normal in the root coordinate frame.
+
+    The plane is defined inside a universe that is filled into a cell rotated
+    by 45 degrees about z, so the normal the surface reports in its own frame
+    is (1, 0, 0) while its normal in the root frame is 45 degrees away from
+    the particle direction. The surface-crossing estimator must use the
+    latter, giving w/|mu| = 1/cos(45 deg).
+    """
+    openmc.reset_auto_ids()
+
+    # Universe split by a plane whose local normal is +x
+    xplane = openmc.XPlane(0.0)
+    inner1 = openmc.Cell(region=-xplane)
+    inner2 = openmc.Cell(region=+xplane)
+    inner_univ = openmc.Universe(cells=[inner1, inner2])
+
+    sph = openmc.Sphere(r=10.0, boundary_type='vacuum')
+    root_cell = openmc.Cell(region=-sph, fill=inner_univ)
+    root_cell.rotation = (0.0, 0.0, 45.0)
+
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([root_cell])
+
+    src = openmc.IndependentSource()
+    src.space = openmc.stats.Point((-5.0, 0.0, 0.0))
+    src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0))
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 1
+    model.settings.particles = 100
+    model.settings.source = src
+
+    flux_tally = openmc.Tally()
+    flux_tally.filters = [openmc.SurfaceFilter([xplane])]
+    flux_tally.scores = ['flux']
+
+    current_tally = openmc.Tally()
+    current_tally.filters = [openmc.SurfaceFilter([xplane])]
+    current_tally.scores = ['current']
+
+    model.tallies = [flux_tally, current_tally]
+    model.run(apply_tally_results=True)
+
+    mu = math.cos(math.radians(45.0))
+    assert flux_tally.mean.flat[0] == pytest.approx(1.0 / mu)
+    # The net current is direction-signed, so it is unaffected by the rotation
+    assert current_tally.mean.flat[0] == pytest.approx(1.0)
+
+
+def test_lattice_surface_flux_rotated(run_in_tmpdir):
+    """Lattice boundary normals are also taken in the root coordinate frame.
+
+    Same setup as test_surface_flux_rotated_universe, except the crossing is a
+    lattice tile boundary rather than a CSG surface. RectLattice reports its
+    boundary normal in the lattice's own frame, which here is rotated by 45
+    degrees relative to the particle direction.
+    """
+    openmc.reset_auto_ids()
+
+    cell1 = openmc.Cell()
+    cell2 = openmc.Cell()
+    cell_outer = openmc.Cell()
+    univ1 = openmc.Universe(cells=[cell1])
+    univ2 = openmc.Universe(cells=[cell2])
+    univ_outer = openmc.Universe(cells=[cell_outer])
+
+    lattice = openmc.RectLattice()
+    lattice.lower_left = (-2.0, -10.0, -10.0)
+    lattice.pitch = (2.0, 20.0, 20.0)
+    lattice.universes = [[[univ1, univ2]]]
+    lattice.outer = univ_outer
+
+    sph = openmc.Sphere(r=10.0, boundary_type='vacuum')
+    root_cell = openmc.Cell(region=-sph, fill=lattice)
+    root_cell.rotation = (0.0, 0.0, 45.0)
+
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([root_cell])
+
+    src = openmc.IndependentSource()
+    src.space = openmc.stats.Point((-5.0, 0.0, 0.0))
+    src.angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0))
+
+    model.settings.run_mode = 'fixed source'
+    model.settings.batches = 1
+    model.settings.particles = 100
+    model.settings.source = src
+
+    tally = openmc.Tally()
+    tally.filters = [openmc.CellFromFilter([cell1]), openmc.CellFilter([cell2])]
+    tally.scores = ['flux']
+    model.tallies = [tally]
+
+    model.run(apply_tally_results=True)
+
+    mu = math.cos(math.radians(45.0))
+    assert tally.mean.flat[0] == pytest.approx(1.0 / mu)