Fix for erroneously non-zero tally results of photon threshold reacti…

…ons (#3242)

Co-authored-by: Paul Romano <[email protected]>

src/photon.cpp

@@ -328,15 +328,19 @@ PhotonInteraction::PhotonInteraction(hid_t group)
   }
 
   // Take logarithm of energies and cross sections since they are log-log
-  // interpolated
+  // interpolated. Note that cross section libraries converted from ACE files
+  // represent zero as exp(-500) to avoid log-log interpolation errors. For
+  // values below exp(-499) we store the log as -900, for which exp(-900)
+  // evaluates to zero.
+  double limit = std::exp(-499.0);
   energy_ = xt::log(energy_);
-  coherent_ = xt::where(coherent_ > 0.0, xt::log(coherent_), -500.0);
-  incoherent_ = xt::where(incoherent_ > 0.0, xt::log(incoherent_), -500.0);
+  coherent_ = xt::where(coherent_ > limit, xt::log(coherent_), -900.0);
+  incoherent_ = xt::where(incoherent_ > limit, xt::log(incoherent_), -900.0);
   photoelectric_total_ = xt::where(
-    photoelectric_total_ > 0.0, xt::log(photoelectric_total_), -500.0);
+    photoelectric_total_ > limit, xt::log(photoelectric_total_), -900.0);
   pair_production_total_ = xt::where(
-    pair_production_total_ > 0.0, xt::log(pair_production_total_), -500.0);
-  heating_ = xt::where(heating_ > 0.0, xt::log(heating_), -500.0);
+    pair_production_total_ > limit, xt::log(pair_production_total_), -900.0);
+  heating_ = xt::where(heating_ > limit, xt::log(heating_), -900.0);
 }
 
 PhotonInteraction::~PhotonInteraction()