Fix a bug in compton scattering shell selection

src/photon.cpp

@@ -466,31 +466,25 @@ void PhotonInteraction::compton_doppler(
       if (rn < c)
         break;
     }
+    double E = alpha * MASS_ELECTRON_EV;
 
     // Determine binding energy of shell
     double E_b = binding_energy_(shell);
 
-    // Determine p_z,max
-    double E = alpha * MASS_ELECTRON_EV;
-    if (E < E_b) {
-      *E_out = alpha / (1 + alpha * (1 - mu)) * MASS_ELECTRON_EV;
-      break;
-    }
+    // Resample if photon energy insufficient
+    if (E < E_b)
+      continue;
 
     double pz_max = -FINE_STRUCTURE * (E_b - (E - E_b) * alpha * (1.0 - mu)) /
                     std::sqrt(2.0 * E * (E - E_b) * (1.0 - mu) + E_b * E_b);
-    if (pz_max < 0.0) {
-      *E_out = alpha / (1 + alpha * (1 - mu)) * MASS_ELECTRON_EV;
-      break;
-    }
 
     // Determine profile cdf value corresponding to p_z,max
     double c_max;
-    if (pz_max > data::compton_profile_pz(n - 1)) {
+    if (std::abs(pz_max) > data::compton_profile_pz(n - 1)) {
       c_max = profile_cdf_(shell, n - 1);
     } else {
       int i = lower_bound_index(data::compton_profile_pz.cbegin(),
-        data::compton_profile_pz.cend(), pz_max);
+        data::compton_profile_pz.cend(), std::abs(pz_max));
       double pz_l = data::compton_profile_pz(i);
       double pz_r = data::compton_profile_pz(i + 1);
       double p_l = profile_pdf_(shell, i);
@@ -499,16 +493,21 @@ void PhotonInteraction::compton_doppler(
       if (pz_l == pz_r) {
         c_max = c_l;
       } else if (p_l == p_r) {
-        c_max = c_l + (pz_max - pz_l) * p_l;
+        c_max = c_l + (std::abs(pz_max) - pz_l) * p_l;
       } else {
         double m = (p_l - p_r) / (pz_l - pz_r);
-        c_max = c_l + (std::pow((m * (pz_max - pz_l) + p_l), 2) - p_l * p_l) /
+        c_max = c_l + (std::pow((m * (std::abs(pz_max) - pz_l) + p_l), 2) -
+                        p_l * p_l) /
                         (2.0 * m);
       }
     }
 
     // Sample value on bounded cdf
-    c = prn(seed) * c_max;
+    if (pz_max > 0.0) {
+      c = prn(seed) * c_max;
+    } else {
+      c = 1.0 + (c_max - 1.0) * prn(seed);
+    }
 
     // Determine pz corresponding to sampled cdf value
     tensor::View<const double> cdf_shell = profile_cdf_.slice(shell);
@@ -533,35 +532,23 @@ void PhotonInteraction::compton_doppler(
     double momentum_sq = std::pow((pz / FINE_STRUCTURE), 2);
     double f = 1.0 + alpha * (1.0 - mu);
     double a = momentum_sq - f * f;
-    double b = 2.0 * E * (f - momentum_sq * mu);
-    c = E * E * (momentum_sq - 1.0);
-
+    double b = 2.0 * (f - momentum_sq * mu);
+    c = (momentum_sq - 1.0);
     double quad = b * b - 4.0 * a * c;
-    if (quad < 0) {
-      *E_out = alpha / (1 + alpha * (1 - mu)) * MASS_ELECTRON_EV;
-      break;
-    }
+    if (quad < 0)
+      continue;
     quad = std::sqrt(quad);
-    double E_out1 = -(b + quad) / (2.0 * a);
-    double E_out2 = -(b - quad) / (2.0 * a);
-
-    // Determine solution to quadratic equation that is positive
-    if (E_out1 > 0.0) {
-      if (E_out2 > 0.0) {
-        // If both are positive, pick one at random
-        *E_out = prn(seed) < 0.5 ? E_out1 : E_out2;
-      } else {
-        *E_out = E_out1;
-      }
-    } else {
-      if (E_out2 > 0.0) {
-        *E_out = E_out2;
-      } else {
-        // No positive solution -- resample
-        continue;
-      }
-    }
-    if (*E_out < E - E_b)
+    double E_out1 = -(b + quad) / (2.0 * a) * E;
+    double E_out2 = -(b - quad) / (2.0 * a) * E;
+    // Determine solution to quadratic equation that is positive and minimal
+    if ((E_out1 < 0.0) && (E_out2 < 0.0))
+      continue;
+    if ((E_out1 > 0.0) && (E_out2 > 0.0))
+      *E_out = std::min(E_out1, E_out2);
+    else
+      *E_out = std::max(E_out1, E_out2);
+
+    if (prn(seed) * E <= *E_out)
       break;
   }
 

tests/regression_tests/atomic_relaxation/results_true.dat

@@ -1,9 +1,9 @@
 tally 1:
-1.956204E+00
-3.826732E+00
-7.918768E+04
-6.270688E+09
+1.920387E+00
+3.687888E+00
+7.854617E+04
+6.169501E+09
 0.000000E+00
 0.000000E+00
-9.208123E+05
-8.478953E+11
+9.214538E+05
+8.490772E+11

tests/regression_tests/photon_production/results_true.dat

@@ -1,8 +1,8 @@
 tally 1:
 8.610000E-01
 7.413210E-01
-9.491000E-01
-9.007908E-01
+9.493000E-01
+9.011705E-01
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -16,12 +16,12 @@ tally 2:
 1.573004E+00
 4.296434E-04
 1.845934E-07
-2.350047E-01
-5.522722E-02
+2.362831E-01
+5.582969E-02
 0.000000E+00
 0.000000E+00
-2.350047E-01
-5.522722E-02
+2.362831E-01
+5.582969E-02
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -53,40 +53,40 @@ tally 3:
 4.104374E+12
 4.296582E-04
 1.846062E-07
-2.297000E-01
-5.276209E-02
-4.196651E+00
-1.761188E+01
+2.344000E-01
+5.494336E-02
+4.593008E+00
+2.109572E+01
 0.000000E+00
 0.000000E+00
-2.297000E-01
-5.276209E-02
-4.196651E+00
-1.761188E+01
+2.344000E-01
+5.494336E-02
+4.593008E+00
+2.109572E+01
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.774484E+05
-3.148794E+10
+1.792949E+05
+3.214667E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.774484E+05
-3.148794E+10
+1.792949E+05
+3.214667E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.474427E+04
-2.173936E+08
+1.452085E+04
+2.108550E+08
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.474427E+04
-2.173936E+08
+1.452085E+04
+2.108550E+08
 0.000000E+00
 0.000000E+00
 tally 4:
@@ -102,16 +102,16 @@ tally 4:
 4.104374E+12
 0.000000E+00
 0.000000E+00
-2.297000E-01
-5.276209E-02
-4.196651E+00
-1.761188E+01
+2.344000E-01
+5.494336E-02
+4.593008E+00
+2.109572E+01
 0.000000E+00
 0.000000E+00
-2.297000E-01
-5.276209E-02
-4.196651E+00
-1.761188E+01
+2.344000E-01
+5.494336E-02
+4.593008E+00
+2.109572E+01
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -122,8 +122,8 @@ tally 4:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.774484E+05
-3.148794E+10
+1.792949E+05
+3.214667E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -134,8 +134,8 @@ tally 4:
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.474427E+04
-2.173936E+08
+1.452085E+04
+2.108550E+08
 0.000000E+00
 0.000000E+00
 tally 5:

tests/regression_tests/photon_production_fission/results_true.dat

@@ -1,12 +1,12 @@
 k-combined:
-2.297165E+00 1.955494E-02
+2.255088E+00 7.363219E-02
 tally 1:
-2.696393E+00
-2.423937E+00
+2.640439E+00
+2.329186E+00
 0.000000E+00
 0.000000E+00
-2.696393E+00
-2.423937E+00
+2.640439E+00
+2.329186E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
@@ -18,52 +18,52 @@ tally 1:
 0.000000E+00
 0.000000E+00
 tally 2:
-2.672029E+00
-2.380251E+00
-4.288244E+08
-6.130569E+16
+2.625521E+00
+2.303289E+00
+4.210498E+08
+5.923082E+16
 0.000000E+00
 0.000000E+00
-2.672029E+00
-2.380251E+00
-4.288244E+08
-6.130569E+16
+2.625521E+00
+2.303289E+00
+4.210498E+08
+5.923082E+16
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.711908E+05
-9.769096E+09
+1.826906E+05
+1.114925E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.711908E+05
-9.769096E+09
+1.826906E+05
+1.114925E+10
 0.000000E+00
 0.000000E+00
 tally 3:
-2.649127E+00
-2.339294E+00
-4.288244E+08
-6.130569E+16
+2.658439E+00
+2.355784E+00
+4.210498E+08
+5.923082E+16
 0.000000E+00
 0.000000E+00
-2.649127E+00
-2.339294E+00
-4.288244E+08
-6.130569E+16
+2.658439E+00
+2.355784E+00
+4.210498E+08
+5.923082E+16
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.711908E+05
-9.769096E+09
+1.826906E+05
+1.114925E+10
 0.000000E+00
 0.000000E+00
 0.000000E+00
 0.000000E+00
-1.711908E+05
-9.769096E+09
+1.826906E+05
+1.114925E+10
 0.000000E+00
 0.000000E+00

tests/regression_tests/photon_source/results_true.dat

@@ -1,5 +1,5 @@
 tally 1:
-2.263761E+02
-5.124615E+04
+2.193647E+02
+4.812086E+04
 0.000000E+00
 0.000000E+00

tests/regression_tests/pulse_height/inputs_true.dat

@@ -15,8 +15,8 @@
   </geometry>
   <settings>
     <run_mode>fixed source</run_mode>
-    <particles>100</particles>
-    <batches>5</batches>
+    <particles>10000</particles>
+    <batches>1</batches>
     <source type="independent" strength="1.0" particle="photon">
       <space type="point">
         <parameters>0.0 0.0 0.0</parameters>