diff --git a/src/particle.cpp b/src/particle.cpp
index 7a003feb9..a0f1d66bd 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -44,34 +44,30 @@ namespace openmc {
 
 double Particle::speed() const
 {
-  // Determine mass in eV/c^2
-  double mass;
-  switch (this->type()) {
-  case ParticleType::neutron:
-    mass = MASS_NEUTRON_EV;
-    break;
-  case ParticleType::photon:
-    mass = 0.0;
-    break;
-  case ParticleType::electron:
-  case ParticleType::positron:
-    mass = MASS_ELECTRON_EV;
-    break;
-  }
-
-  if (this->E() < 1.0e-9 * mass) {
-    // If the energy is much smaller than the mass, revert to non-relativistic
-    // formula. The 1e-9 criterion is specifically chosen as the point below
-    // which the error from using the non-relativistic formula is less than the
-    // round-off eror when using the relativistic formula (see analysis at
-    // https://gist.github.com/paulromano/da3b473fe3df33de94b265bdff0c7817)
-    return C_LIGHT * std::sqrt(2 * this->E() / mass);
+  if (settings::run_CE) {
+    // Determine mass in eV/c^2
+    double mass;
+    switch (this->type()) {
+    case ParticleType::neutron:
+      mass = MASS_NEUTRON_EV;
+      break;
+    case ParticleType::photon:
+      mass = 0.0;
+      break;
+    case ParticleType::electron:
+    case ParticleType::positron:
+      mass = MASS_ELECTRON_EV;
+      break;
+    }
+    // Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<<m:
+    return C_LIGHT * std::sqrt(this->E() * (this->E() + 2 * mass)) /
+           (this->E() + mass);
   } else {
-    // Calculate inverse of Lorentz factor
-    const double inv_gamma = mass / (this->E() + mass);
-
-    // Calculate speed via v = c * sqrt(1 - γ^-2)
-    return C_LIGHT * std::sqrt(1 - inv_gamma * inv_gamma);
+    auto& macro_xs = data::mg.macro_xs_[this->material()];
+    int macro_t = this->mg_xs_cache().t;
+    int macro_a = macro_xs.get_angle_index(this->u());
+    return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, this->g(), nullptr,
+                   nullptr, nullptr, macro_t, macro_a);
   }
 }
 
@@ -241,8 +237,14 @@ void Particle::event_advance()
     collision_distance() = -std::log(prn(current_seed())) / macro_xs().total;
   }
 
-  // Select smaller of the two distances
-  double distance = std::min(boundary().distance(), collision_distance());
+  double speed = this->speed();
+  double time_cutoff = settings::time_cutoff[static_cast<int>(type())];
+  double distance_cutoff =
+    (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
+
+  // Select smaller of the three distances
+  double distance =
+    std::min({boundary().distance(), collision_distance(), distance_cutoff});
 
   // Advance particle in space and time
   // Short-term solution until the surface source is revised and we can use
@@ -250,23 +252,10 @@ void Particle::event_advance()
   for (int j = 0; j < n_coord(); ++j) {
     coord(j).r() += distance * coord(j).u();
   }
-  double dt = distance / this->speed();
+  double dt = distance / speed;
   this->time() += dt;
   this->lifetime() += dt;
 
-  // Kill particle if its time exceeds the cutoff
-  bool hit_time_boundary = false;
-  double time_cutoff = settings::time_cutoff[static_cast<int>(type())];
-  if (time() > time_cutoff) {
-    double dt = time() - time_cutoff;
-    time() = time_cutoff;
-    lifetime() = time_cutoff;
-
-    double push_back_distance = speed() * dt;
-    this->move_distance(-push_back_distance);
-    hit_time_boundary = true;
-  }
-
   // Score track-length tallies
   if (!model::active_tracklength_tallies.empty()) {
     score_tracklength_tally(*this, distance);
@@ -284,7 +273,7 @@ void Particle::event_advance()
   }
 
   // Set particle weight to zero if it hit the time boundary
-  if (hit_time_boundary) {
+  if (distance == distance_cutoff) {
     wgt() = 0.0;
   }
 }
diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp
index c97a3d6f3..5ebef9c14 100644
--- a/src/physics_mg.cpp
+++ b/src/physics_mg.cpp
@@ -137,6 +137,7 @@ void create_fission_sites(Particle& p)
     SourceSite site;
     site.r = p.r();
     site.particle = ParticleType::neutron;
+    site.time = p.time();
     site.wgt = 1. / weight;
     site.parent_id = p.id();
     site.progeny_id = p.n_progeny()++;
diff --git a/tests/regression_tests/time_cutoff/results_true.dat b/tests/regression_tests/time_cutoff/results_true.dat
index 1cafceb77..d3d5e1b3c 100644
--- a/tests/regression_tests/time_cutoff/results_true.dat
+++ b/tests/regression_tests/time_cutoff/results_true.dat
@@ -1,5 +1,5 @@
 tally 1:
-2.000000E+03
-4.000000E+05
+1.383148E+02
+1.913099E+03
 0.000000E+00
 0.000000E+00
diff --git a/tests/unit_tests/test_mesh.py b/tests/unit_tests/test_mesh.py
index 3c8e988c0..67ca4028e 100644
--- a/tests/unit_tests/test_mesh.py
+++ b/tests/unit_tests/test_mesh.py
@@ -620,7 +620,7 @@ def test_mesh_material_volumes_serialize():
     assert new_volumes.by_element(3) == [(2, 1.0)]
 
 
-def test_raytrace_mesh_infinite_loop():
+def test_raytrace_mesh_infinite_loop(run_in_tmpdir):
     # Create a model with one large spherical cell
     sphere = openmc.Sphere(r=100, boundary_type='vacuum')
     cell = openmc.Cell(region=-sphere)