diff --git a/src/photon.cpp b/src/photon.cpp
index 1be515a62..bb715cdea 100644
--- a/src/photon.cpp
+++ b/src/photon.cpp
@@ -128,7 +128,7 @@ PhotonInteraction::PhotonInteraction(hid_t group)
   }
 
   shells_.resize(n_shell);
-  cross_sections_.resize({energy_.size(), n_shell});
+  cross_sections_ = xt::zeros<double>({energy_.size(), n_shell});
 
   // Create mapping from designator to index
   std::unordered_map<int, int> shell_map;
@@ -166,8 +166,8 @@ PhotonInteraction::PhotonInteraction(hid_t group)
     read_dataset(tgroup, "xs", xs);
 
     auto cross_section = xt::view(cross_sections_,
-      xt::range(shell.threshold, shell.threshold + xs.size()), i);
-    cross_section = xt::where(xs > 0.0, xt::log(xs), -500.0);
+      xt::range(shell.threshold, xt::placeholders::_), i);
+    cross_section = xt::where(xs > 0, xt::log(xs), 0);
 
     if (object_exists(tgroup, "transitions")) {
       // Determine dimensions of transitions
@@ -571,10 +571,13 @@ void PhotonInteraction::calculate_xs(Particle& p) const
     incoherent_(i_grid) + f * (incoherent_(i_grid + 1) - incoherent_(i_grid)));
 
   // Calculate microscopic photoelectric cross section
-  const auto& xs_upper = xt::row(cross_sections_, i_grid);
-  const auto& xs_lower = xt::row(cross_sections_, i_grid + 1);
+  xs.photoelectric = 0.0;
+  const auto& xs_lower = xt::row(cross_sections_, i_grid);
+  const auto& xs_upper = xt::row(cross_sections_, i_grid + 1);
 
-  xs.photoelectric = xt::sum(xt::exp(xs_upper + f * (xs_upper - xs_lower)))[0];
+  for (int i = 0; i < xs_upper.size(); ++i)
+    if (xs_lower(i) != 0)
+      xs.photoelectric += std::exp(xs_lower(i) + f * (xs_upper(i) - xs_lower(i)));
 
   // Calculate microscopic pair production cross section
   xs.pair_production = std::exp(
diff --git a/src/physics.cpp b/src/physics.cpp
index 38d74dbed..26b3a3520 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -29,6 +29,7 @@
 
 #include <algorithm> // for max, min, max_element
 #include <cmath>     // for sqrt, exp, log, abs, copysign
+#include <xtensor/xview.hpp>
 
 namespace openmc {
 
@@ -344,25 +345,24 @@ void sample_photon_reaction(Particle& p)
 
   // Photoelectric effect
   double prob_after = prob + micro.photoelectric;
+  int i_grid = micro.index_grid;
+  double f = micro.interp_factor;
+  const auto& xs_lower = xt::row(element.cross_sections_, i_grid);
+  const auto& xs_upper = xt::row(element.cross_sections_, i_grid + 1);
+
   if (prob_after > cutoff) {
     for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) {
       const auto& shell {element.shells_[i_shell]};
 
-      // Get grid index and interpolation factor
-      int i_grid = micro.index_grid;
-      double f = micro.interp_factor;
-
       // Check threshold of reaction
-      int i_start = shell.threshold;
-      if (i_grid < i_start)
+      if (xs_lower(i_shell) == 0)
         continue;
 
-      // Evaluation subshell photoionization cross section
-      double xs = std::exp(element.cross_sections_(i_grid - i_start) +
-                           f * (element.cross_sections_(i_grid + 1 - i_start) -
-                                 element.cross_sections_(i_grid - i_start)));
+      //  Evaluation subshell photoionization cross section
+      prob += std::exp(xs_lower(i_shell) +
+                       f * (xs_upper(i_shell) -
+                             xs_lower(i_shell)));
 
-      prob += xs;
       if (prob > cutoff) {
         double E_electron = p.E() - shell.binding_energy;