Fix incorrect LCFA pair-creation spectrum

Error in positron energy spectrum visible for chi > 2, add test for T(chi) to guarantee this is sorted.
tgblackburn · Jan 23, 2023 · 307619f · 307619f
1 parent aa118c0
commit 307619f
Show file tree

Hide file tree

Showing 4 changed files with 93 additions and 8 deletions.
diff --git a/Cargo.lock b/Cargo.lock
diff --git a/Cargo.toml b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "ptarmigan"
-version = "1.2.0"
+version = "1.2.1"
 authors = ["tgblackburn"]
 edition = "2018"
 publish = false

diff --git a/docs/changelog.md b/docs/changelog.md
@@ -1,13 +1,17 @@
 # Changelog
 
-## v1.2.0
+## v1.2.1
 
 2023-01-23
 
 Added:
 
 * Looping over a range of a0s (PR #45).
 
+Fixed:
+
+* Incorrect sampling of LCFA pair-creation spectrum (evident only for chi > 2).
+
 ## v1.1.0
 
 2023-01-13

diff --git a/src/lcfa/pair_creation/mod.rs b/src/lcfa/pair_creation/mod.rs
@@ -46,7 +46,7 @@ fn spectrum(s: f64, chi: f64) -> f64 {
         .map(|(t, w)| {
             let xi = 2.0 / (3.0 * chi * s * (1.0 - s));
             let prefactor = (-xi * t.cosh() + t).exp();
-            w * prefactor * ((s / (1.0 - s) + (1.0 - s) / s) * (1.5 * t).cosh() + (t / 3.0).cosh() / t.cosh())
+            w * prefactor * ((s / (1.0 - s) + (1.0 - s) / s) * (2.0 * t / 3.0).cosh() + (t / 3.0).cosh() / t.cosh())
         })
         .sum()
 }
@@ -68,10 +68,10 @@ fn angular_spectrum(z: f64, s: f64, chi: f64) -> f64 {
 /// the cosine of the scattering angle, as well as the
 /// equivalent s and z for debugging purposes
 pub fn sample<R: Rng>(chi: f64, gamma: f64, rng: &mut R) -> (f64, f64, f64, f64) {
-    let max = if chi < 2.5 {
+    let max = if chi < 2.0 {
         spectrum(0.5, chi)
     } else {
-        spectrum(1.2 / chi, chi)
+        spectrum(0.9 * chi.powf(-0.875), chi)
     };
     let max = 1.2 * max;
 
@@ -151,16 +151,27 @@ mod tests {
 
         for (chi, target) in &pts {
             let result = auxiliary_t(*chi);
+
+            let prefactor = 1.0 / (3_f64.sqrt() * std::f64::consts::PI * chi);
+            let intgd: f64 = GAUSS_32_NODES.iter()
+                .zip(GAUSS_32_WEIGHTS.iter())
+                .map(|(t, w)| {
+                    let s = 0.5 * (1.0 + t);
+                    prefactor * 0.5 * w * spectrum(s, *chi)
+                })
+                .sum();
+
             let error = (result - target).abs() / target;
-            //println!("chi = {:.3e}, t(chi) = {:.6e}, error = {:.3e}", chi, result, error);
+            let intgd_error = (intgd - target).abs() / target;
+            println!("chi = {:>9.3e}, t(chi) = {:>12.6e} | {:>12.6e} [interp|intgd], error = {:.3e} | {:.3e}", chi, result, intgd, error, intgd_error);
             assert!(error < max_error);
         }
     }
 
     #[test]
     #[ignore]
     fn pair_spectrum_sampling() {
-        let chi = 1.0;
+        let chi = 2.0;
         let gamma = 1000.0;
         let mut rng = Xoshiro256StarStar::seed_from_u64(0);
         let path = format!("output/lcfa_pair_spectrum_{}.dat", chi);
@@ -172,4 +183,74 @@ mod tests {
             writeln!(file, "{:.6e} {:.6e}", s, z).unwrap();
         }
     }
+
+    static GAUSS_32_NODES: [f64; 32] = [
+        -9.972638618494816e-1,
+        -9.856115115452683e-1,
+        -9.647622555875064e-1,
+        -9.349060759377397e-1,
+        -8.963211557660521e-1,
+        -8.493676137325700e-1,
+        -7.944837959679424e-1,
+        -7.321821187402897e-1,
+        -6.630442669302152e-1,
+        -5.877157572407623e-1,
+        -5.068999089322294e-1,
+        -4.213512761306353e-1,
+        -3.318686022821276e-1,
+        -2.392873622521371e-1,
+        -1.444719615827965e-1,
+        -4.830766568773832e-2,
+        4.830766568773832e-2,
+        1.444719615827965e-1,
+        2.392873622521371e-1,
+        3.318686022821276e-1,
+        4.213512761306353e-1,
+        5.068999089322294e-1,
+        5.877157572407623e-1,
+        6.630442669302152e-1,
+        7.321821187402897e-1,
+        7.944837959679424e-1,
+        8.493676137325700e-1,
+        8.963211557660521e-1,
+        9.349060759377397e-1,
+        9.647622555875064e-1,
+        9.856115115452683e-1,
+        9.972638618494816e-1,
+    ];
+
+    static GAUSS_32_WEIGHTS: [f64; 32] = [
+        7.018610000000000e-3,
+        1.627439500000000e-2,
+        2.539206500000000e-2,
+        3.427386300000000e-2,
+        4.283589800000000e-2,
+        5.099805900000000e-2,
+        5.868409350000000e-2,
+        6.582222280000000e-2,
+        7.234579411000000e-2,
+        7.819389578700000e-2,
+        8.331192422690000e-2,
+        8.765209300440000e-2,
+        9.117387869576400e-2,
+        9.384439908080460e-2,
+        9.563872007927486e-2,
+        9.654008851472780e-2,
+        9.654008851472780e-2,
+        9.563872007927486e-2,
+        9.384439908080460e-2,
+        9.117387869576400e-2,
+        8.765209300440000e-2,
+        8.331192422690000e-2,
+        7.819389578700000e-2,
+        7.234579411000000e-2,
+        6.582222280000000e-2,
+        5.868409350000000e-2,
+        5.099805900000000e-2,
+        4.283589800000000e-2,
+        3.427386300000000e-2,
+        2.539206500000000e-2,
+        1.627439500000000e-2,
+        7.018610000000000e-3,
+    ];
 }