From b79d97d6ddecf22622abb75ce125cdc017c75afd Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 4 Dec 2019 13:14:34 -0600
Subject: [PATCH] Fix a number of problems related to photoatomic data

- OpenMC segfaulted when subshell cross sections were missing in HDF5 files
- IncidentPhoton.from_ace shouldn't work with mcplib84 and before
- Properly set threshold for subshell cross sections from eprdata12/14
- Make sure IncidentPhoton.from_endf works with only photoatomic data
---
 openmc/data/photon.py | 16 +++++++++++++---
 src/photon.cpp        |  5 +++++
 2 files changed, 18 insertions(+), 3 deletions(-)

diff --git a/openmc/data/photon.py b/openmc/data/photon.py
index 3766b697e..31514639a 100644
--- a/openmc/data/photon.py
+++ b/openmc/data/photon.py
@@ -616,17 +616,23 @@ class IncidentPhoton(EqualityMixin):
                 rx = PhotonReaction(mt)
                 data.reactions[mt] = rx
 
-                # Store cross section
+                # Store cross section, determining threshold
                 xs = ace.xss[idx : idx+n_energy].copy()
                 nonzero = (xs != 0.0)
                 xs[nonzero] = np.exp(xs[nonzero])
-                rx.xs = Tabulated1D(energy, xs, [n_energy], [5])
+                threshold = np.where(xs > 0.0)[0][0]
+                rx.xs = Tabulated1D(energy[threshold:], xs[threshold:],
+                                    [n_energy - threshold], [5])
                 idx += n_energy
 
                 # Copy binding energy
                 shell = _SUBSHELLS[d]
                 e = data.atomic_relaxation.binding_energy[shell]
                 rx.subshell_binding_energy = e
+        else:
+            raise ValueError("ACE table {} does not have subshell data. Only "
+                             "newer ACE photoatomic libraries are supported "
+                             "(e.g., eprdata14).".format(ace.name))
 
         # Add bremsstrahlung DCS data
         data._add_bremsstrahlung()
@@ -990,7 +996,11 @@ class IncidentPhoton(EqualityMixin):
             for mt, rx in self.reactions.items():
                 if mt >= 534 and mt <= 572:
                     shell = _REACTION_NAME[mt][1]
-                    e_f = self.atomic_relaxation.energy_fluorescence(shell)
+                    relax_data = self.atomic_relaxation
+                    if relax_data is not None:
+                        e_f = relax_data.energy_fluorescence(shell)
+                    else:
+                        e_f = 0.0
                     heating_xs += (energy - e_f) * rx.xs(energy)
 
         heat_rx = PhotonReaction(525)
diff --git a/src/photon.cpp b/src/photon.cpp
index ddde4bd6c..c2f41588a 100644
--- a/src/photon.cpp
+++ b/src/photon.cpp
@@ -110,6 +110,11 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
   std::vector<std::string> designators;
   read_attribute(rgroup, "designators", designators);
   auto n_shell = designators.size();
+  if (n_shell == 0) {
+    throw std::runtime_error{"Photoatomic data for " + name_ +
+      " does not have subshell data."};
+  }
+
   for (int i = 0; i < n_shell; ++i) {
     const auto& designator {designators[i]};