diff --git a/openmc/data/photon.py b/openmc/data/photon.py
index f8400b736d..4e266bc9d7 100644
--- a/openmc/data/photon.py
+++ b/openmc/data/photon.py
@@ -513,10 +513,6 @@ class IncidentPhoton(EqualityMixin):
         for mt in (502, 504, 515, 522, 525):
             data.reactions[mt] = PhotonReaction.from_ace(ace, mt)
 
-        # Get heating cross sections [eV*barn] from factors [MeV per collision]
-        data.reactions[525].xs.y *= sum([data.reactions[mt].xs.y for mt in
-                                         (502, 504, 515, 522)]) * EV_PER_MEV
-
         # Compton profiles
         n_shell = ace.nxs[5]
         if n_shell != 0:
@@ -997,7 +993,17 @@ class PhotonReaction(EqualityMixin):
 
         # Store cross section
         xs = ace.xss[idx : idx+n].copy()
-        if mt in (502, 504, 515, 522):
+        if mt == 525:
+            # Get heating cross sections [eV barn] from factors [MeV/collision]
+            # by multiplying the total xs (sum of (502, 504, 515, 522))
+            idx = ace.jxs[1] + n
+            total_xs = ace.xss[idx : idx + 4*n].copy()
+            total_xs = total_xs.reshape((4, n))
+            nonzero = (total_xs != 0.0)
+            total_xs[nonzero] = np.exp(total_xs[nonzero])
+            total_xs = total_xs.sum(axis=0)
+            xs *= total_xs * EV_PER_MEV
+        else:
             nonzero = (xs != 0.0)
             xs[nonzero] = np.exp(xs[nonzero])
         rx.xs = Tabulated1D(energy, xs, [n], [5])