diff --git a/‎include/openmc/photon.h‎
Lines changed: 1 addition & 1 deletion b/‎include/openmc/photon.h‎
Lines changed: 1 addition & 1 deletion
diff --git a/‎src/photon.cpp‎
Lines changed: 11 additions & 15 deletions b/‎src/photon.cpp‎
Lines changed: 11 additions & 15 deletions
diff --git a/‎src/physics.cpp‎
Lines changed: 3 additions & 3 deletions b/‎src/physics.cpp‎
Lines changed: 3 additions & 3 deletions
@@ -36,7 +36,6 @@ class ElectronSubshell {
   int threshold;
   double n_electrons;
   double binding_energy;
-  xt::xtensor<double, 1> cross_section;
   vector<Transition> transitions;
 };
 
@@ -82,6 +81,7 @@ class PhotonInteraction {
 
   // Photoionization and atomic relaxation data
   vector<ElectronSubshell> shells_;
+  xt::xtensor<double, 2> cross_sections_;
 
   // Compton profile data
   xt::xtensor<double, 2> profile_pdf_;
 
@@ -16,6 +16,9 @@
 #include "xtensor/xbuilder.hpp"
 #include "xtensor/xoperation.hpp"
 #include "xtensor/xview.hpp"
+#include "xtensor/xmath.hpp"
+#include "xtensor/xslice.hpp"
+#include "xtensor/xtensor_forward.hpp"
 
 #include <cmath>
 #include <fmt/core.h>
@@ -125,6 +128,7 @@ PhotonInteraction::PhotonInteraction(hid_t group)
   }
 
   shells_.resize(n_shell);
+  cross_sections_.resize({energy_.size(), n_shell});
 
   // Create mapping from designator to index
   std::unordered_map<int, int> shell_map;
@@ -155,13 +159,14 @@ PhotonInteraction::PhotonInteraction(hid_t group)
     read_attribute(tgroup, "num_electrons", shell.n_electrons);
 
     // Read subshell cross section
+    xt::xtensor<double, 1> xs;
     dset = open_dataset(tgroup, "xs");
     read_attribute(dset, "threshold_idx", shell.threshold);
     close_dataset(dset);
-    read_dataset(tgroup, "xs", shell.cross_section);
+    read_dataset(tgroup, "xs", xs);
 
-    auto& xs = shell.cross_section;
-    xs = xt::where(xs > 0.0, xt::log(xs), -500.0);
+    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);
 
     if (object_exists(tgroup, "transitions")) {
       // Determine dimensions of transitions
@@ -565,19 +570,10 @@ void PhotonInteraction::calculate_xs(Particle& p) const
     incoherent_(i_grid) + f * (incoherent_(i_grid + 1) - incoherent_(i_grid)));
 
   // Calculate microscopic photoelectric cross section
-  xs.photoelectric = 0.0;
-  for (const auto& shell : shells_) {
-    // Check threshold of reaction
-    int i_start = shell.threshold;
-    if (i_grid < i_start)
-      continue;
+  const auto& xs_upper = xt::row(cross_sections_, i_grid);
+  const auto& xs_lower = xt::row(cross_sections_, i_grid + 1);
 
-    // Evaluation subshell photoionization cross section
-    xs.photoelectric +=
-      std::exp(shell.cross_section(i_grid - i_start) +
-               f * (shell.cross_section(i_grid + 1 - i_start) -
-                     shell.cross_section(i_grid - i_start)));
-  }
+  xs.photoelectric = xt::sum(xt::exp(xs_upper + f * (xs_upper - xs_lower)))[0];
 
   // Calculate microscopic pair production cross section
   xs.pair_production = std::exp(
 
@@ -358,9 +358,9 @@ void sample_photon_reaction(Particle& p)
         continue;
 
       // Evaluation subshell photoionization cross section
-      double xs = std::exp(shell.cross_section(i_grid - i_start) +
-                           f * (shell.cross_section(i_grid + 1 - i_start) -
-                                 shell.cross_section(i_grid - i_start)));
+      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)));
 
       prob += xs;
       if (prob > cutoff) {