From f81f5b164f473a9ed64437e20c67feb98aa9a1b3 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 29 Jun 2020 09:42:15 -0500
Subject: [PATCH] Fix fission delayed photon scaling

---
 include/openmc/nuclide.h |  2 ++
 src/nuclide.cpp          |  8 +++-----
 src/physics.cpp          | 16 +++++++++++++++-
 3 files changed, 20 insertions(+), 6 deletions(-)

diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h
index 60bac81d2d..f3c811b516 100644
--- a/include/openmc/nuclide.h
+++ b/include/openmc/nuclide.h
@@ -82,6 +82,8 @@ public:
   std::unique_ptr<Function1D> total_nu_; //!< Total neutron yield
   std::unique_ptr<Function1D> fission_q_prompt_; //!< Prompt fission energy release
   std::unique_ptr<Function1D> fission_q_recov_; //!< Recoverable fission energy release
+  std::unique_ptr<Function1D> prompt_photons_; //!< Prompt photon energy release
+  std::unique_ptr<Function1D> delayed_photons_; //!< Delayed photon energy release
 
   // Resonance scattering information
   bool resonant_ {false};
diff --git a/src/nuclide.cpp b/src/nuclide.cpp
index 4e5d41b4be..8e199cb122 100644
--- a/src/nuclide.cpp
+++ b/src/nuclide.cpp
@@ -262,8 +262,6 @@ Nuclide::Nuclide(hid_t group, const std::vector<double>& temperature)
   }
 
   // Read fission energy release data if present
-  std::unique_ptr<Function1D> prompt_photons;
-  std::unique_ptr<Function1D> delayed_photons;
   if (object_exists(group, "fission_energy_release")) {
     hid_t fer_group = open_group(group, "fission_energy_release");
     fission_q_prompt_ = read_function(fer_group, "q_prompt");
@@ -271,12 +269,12 @@ Nuclide::Nuclide(hid_t group, const std::vector<double>& temperature)
 
     // We need prompt/delayed photon energy release for scaling fission photon
     // production
-    prompt_photons = read_function(fer_group, "prompt_photons");
-    delayed_photons = read_function(fer_group, "delayed_photons");
+    prompt_photons_ = read_function(fer_group, "prompt_photons");
+    delayed_photons_ = read_function(fer_group, "delayed_photons");
     close_group(fer_group);
   }
 
-  this->create_derived(prompt_photons.get(), delayed_photons.get());
+  this->create_derived(prompt_photons_.get(), delayed_photons_.get());
 }
 
 Nuclide::~Nuclide()
diff --git a/src/physics.cpp b/src/physics.cpp
index 1000c3c4b0..e3ede4d61b 100644
--- a/src/physics.cpp
+++ b/src/physics.cpp
@@ -4,6 +4,7 @@
 #include "openmc/bremsstrahlung.h"
 #include "openmc/constants.h"
 #include "openmc/eigenvalue.h"
+#include "openmc/endf.h"
 #include "openmc/error.h"
 #include "openmc/material.h"
 #include "openmc/math_functions.h"
@@ -567,8 +568,21 @@ void sample_photon_product(int i_nuclide, Particle& p, int* i_rx, int* i_product
 
     for (int j = 0; j < rx->products_.size(); ++j) {
       if (rx->products_[j].particle_ == Particle::Type::photon) {
+        // For fission, artificially increase the photon yield to account
+        // for delayed photons
+        double f = 1.0;
+        if (settings::delayed_photon_scaling) {
+          if (is_fission(rx->mt_)) {
+            if (nuc->prompt_photons_ && nuc->delayed_photons_) {
+              double energy_prompt = (*nuc->prompt_photons_)(p.E_);
+              double energy_delayed = (*nuc->delayed_photons_)(p.E_);
+              f = (energy_prompt + energy_delayed)/(energy_prompt);
+            }
+          }
+        }
+
         // add to cumulative probability
-        prob += (*rx->products_[j].yield_)(p.E_) * xs;
+        prob += f * (*rx->products_[j].yield_)(p.E_) * xs;
 
         *i_rx = i;
         *i_product = j;